A proppant transport simulation method based on a physically encoded neural network

Abstract

The present application relates to a kind of physical encoding neural network-based proppant migration simulation method, it is related to artificial intelligence and the cross technical field of oil and gas hydraulic fracturing engineering, including: the concentration field of proppant in hydraulic fracture and the physical field of fluid in current time are input pressure reasoning network, the pressure field of current time fluid is obtained;Proppant velocity field is calculated based on pressure field in current time;Concentration field and velocity field in current time are input concentration reasoning network, and the concentration field and physical field of next time are obtained;Next time is regarded as new current time, and the steps of above-mentioned pressure field solution, velocity field calculation and concentration field and physical field update are iteratively executed, until reaching preset simulation length, and the concentration field simulation result of proppant in each time is obtained.

Description

Technical Field

[0001] The embodiments in this specification relate to the interdisciplinary technical field of artificial intelligence and oil and gas hydraulic fracturing engineering, specifically to a proppant transport simulation method based on a physical coding neural network. Background Technology

[0002] In the process of hydraulic fracturing of shale gas and other low-permeability oil and gas reservoirs, the distribution pattern of proppant in the fracture space directly determines the oil and gas extraction efficiency and economic benefits. Therefore, studying the migration and diffusion law of proppant has important engineering value.

[0003] Existing proppant transport simulation methods cannot simultaneously achieve both computational efficiency and physical consistency. Traditional numerical simulations offer strong physical consistency but are computationally expensive, making it difficult to support rapid iterations; pure data-driven models are fast but do not incorporate physical laws, resulting in severe error accumulation and insufficient generalization ability in long-term simulations; PINN treats physical constraints as soft penalties, leading to gradient imbalances in strongly convection-dominated problems and difficulty in capturing high-frequency abrupt changes, resulting in insufficient accuracy in long-term simulations. Summary of the Invention

[0004] The purpose of the embodiments in this specification is to provide a proppant transport simulation method based on a physical coding neural network, so as to overcome the problem that existing methods cannot simultaneously achieve both computational efficiency and physical consistency.

[0005] To address the aforementioned technical problems, this specification provides, in one aspect, a proppant transport simulation method based on a physical coding neural network, comprising: The concentration field of the proppant and the physical property field of the fluid in the hydraulic fracture at the current moment are input into the pressure inference network to obtain the pressure field of the fluid at the current moment. The pressure inference network includes a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the degree of deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy the preset isobaric boundary conditions and / or zero flux boundary conditions. Calculate the velocity field of the proppant at the current moment based on the pressure field; The current concentration field and velocity field are input into the concentration inference network to obtain the concentration field and physical property field at the next moment. The concentration inference network is coupled with multiple Π-type convolutional modules. The Π-type convolutional modules include differential convolutional layers and Π-type convolutional layers. The differential convolutional layers are used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolutional layers are used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution. The rheological term is used to characterize the rheological feedback effect of proppant viscosity changes and particle interactions on the concentration distribution. The next moment is taken as the new current moment, and the above steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field are iteratively executed until the preset simulation time is reached, so as to obtain the simulation results of the proppant concentration field at each moment.

[0006] Furthermore, embodiments of this specification provide a proppant transport simulation device, comprising: The pressure inference module is used to input the concentration field of the proppant and the physical property field of the fluid in the hydraulic fracture at the current moment into the pressure inference network to obtain the pressure field of the fluid at the current moment. The pressure inference network includes a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the degree of deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy the preset isobaric boundary conditions and / or zero flux boundary conditions. The calculation module is used to calculate the velocity field of the proppant at the current moment based on the pressure field; The concentration inference module is used to input the current concentration field and velocity field into the concentration inference network to obtain the concentration field and physical property field at the next moment. The concentration inference network is coupled with multiple Π-type convolution modules. The Π-type convolution modules include differential convolution layers and Π-type convolution layers. The differential convolution layers are used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolution layers are used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution. The rheological term is used to characterize the rheological feedback effect of proppant viscosity changes and particle interactions on the concentration distribution. The iteration module is used to take the next moment as the new current moment and iteratively execute the above steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field until the preset simulation time is reached, so as to obtain the simulation results of the proppant concentration field at each moment.

[0007] Furthermore, embodiments of this specification provide a computer device, including: Memory, used to store computer programs; A processor is used to execute the computer program to implement the above-described proppant transport simulation method.

[0008] In another aspect, embodiments of this specification provide a computer storage medium storing computer program instructions that, when executed, implement the aforementioned proppant transport simulation method.

[0009] In another aspect, embodiments of this specification provide a computer program product, including a computer program that, when executed by a processor, implements the above-described proppant transport simulation method.

[0010] As can be seen from the technical solutions provided in the embodiments of this specification above, the embodiments of this specification, through the combined effect of the Poisson regularization term and the boundary mask term in the pressure inference network, ensure that the solution of the pressure field is simultaneously constrained by the elliptic pressure Poisson equation and the hard constraint of the boundary conditions. During long-term time-series simulations, the pressure field always conforms to physical laws and the actual engineering boundary conditions. Furthermore, in the concentration inference network, the differential convolutional layer performs deterministic calculations of the convection term with fixed weights, and the Π-type convolutional layer approximates the rheological term in the form of a high-order polynomial product. These two layers respectively handle the explicitly expressible linear differential terms and the complex nonlinear terms in the physical equations, achieving accurate recursion of the concentration field with a clear division of labor. In addition, through closed-loop iteration of pressure, velocity, concentration, and physical properties, the simulation results of the concentration field at each time step are output, achieving a balance between computational efficiency and physical consistency, and ensuring the accuracy and stability of long-term time-series simulations. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments or prior art of this specification, the accompanying drawings used in the description of the embodiments or prior art will be briefly introduced below.

[0012] Figure 1 This is a flowchart of a proppant transport simulation method based on a physical coding neural network provided in the embodiments of this specification; Figure 2 A complete implementation flowchart of a specific embodiment provided in this specification; Figure 3 This is a schematic diagram of the internal architecture of the Π-type convolution module in a specific embodiment provided in this specification. Figure 4 A schematic diagram of the stress inference network structure and boundary mask layer in a specific embodiment provided in this specification; Figure 5 This is a diagram of the time-flow loop architecture of a coupled system provided in a specific embodiment of the embodiments of this specification; Figure 6 Forward difference in one specific embodiment provided in this specification (·) / x of 3 3. Schematic diagram of the differential convolution operator; Figure 7 Forward difference in one specific embodiment provided in this specification (·) / y's 3 3. Schematic diagram of the differential convolution operator; Figure 8 A schematic diagram of the reference concentration field of FDM in a specific embodiment provided in this specification; Figure 9A schematic diagram of the predicted concentration field of the prior art ConvLSTM in a specific embodiment provided for the embodiments of this specification; Figure 10 A schematic diagram of the prediction absolute error of the prior art ConvLSTM in a specific embodiment provided for the embodiments of this specification; Figure 11 A schematic diagram of the prediction relative error of the prior art ConvLSTM in a specific embodiment provided in this specification; Figure 12 A schematic diagram of the predicted concentration field of prior art PINN in a specific embodiment provided for the embodiments of this specification; Figure 13 A schematic diagram of the prediction absolute error of prior art PINN in a specific embodiment provided for the embodiments of this specification; Figure 14 A schematic diagram illustrating the prediction relative error of prior art PINN in a specific embodiment provided for the present specification. Figure 15 A schematic diagram of the predicted concentration field of the proppant transport simulation method of the present invention in a specific embodiment provided in this specification; Figure 16 A schematic diagram of the prediction absolute error of the proppant transport simulation method of the present invention in a specific embodiment provided in this specification; Figure 17 A schematic diagram of the relative prediction error of the proppant migration simulation method of the present invention in a specific embodiment provided in this specification; Figure 18 This is a schematic diagram of the structural composition of a proppant transport simulation device provided in the embodiments of this specification; Figure 19 This is a schematic diagram of the structural composition of a computer device provided in the embodiments of this specification. Detailed Implementation

[0013] The technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this specification, and not all embodiments. Based on the embodiments in this specification, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this specification.

[0014] It should be noted that the terms "first," "second," etc., used in this specification, claims, and the foregoing drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, apparatus, product, or device that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.

[0015] In some embodiments, hydraulic fractures can be tensile fracture surfaces formed in oil and gas reservoirs by high-pressure fluids during fracturing operations, filled with proppant to maintain their conductivity. The geometry of hydraulic fractures is typically described by three spatial dimensions: length, height, and aperture (width), where the width is much smaller than the length and height, thus simplifying the fluid motion within the fracture into a two-dimensional flow problem. The fracture walls are solid-liquid interfaces, impenetrable to fluid, and the fracture tip is a pressure boundary.

[0016] In some embodiments, the pressure field of the fluid within a hydraulic fracture can be the pressure distribution of the fluid at various spatial locations within the fracture. The pressure field is determined by the fluid injection, fracture geometry, and physical properties, and its spatial distribution affects the direction and rate of fluid transport. Pressure propagates much faster within the fracture than the fluid flow velocity, and changes in the pressure field are instantly transmitted throughout the entire fracture region, thereby driving fluid flow.

[0017] In some embodiments, the physical property field of the fluid within a hydraulic fracture can be the spatial distribution of the physical properties of the mixed fluid within the fracture, including the density field and viscosity field. The mixed fluid consists of pure fracturing fluid and proppant particles suspended therein; therefore, the physical property field depends on the concentration distribution of the proppant. Spatial non-uniformity of density and viscosity directly affects the solution of the pressure field and the reconstruction of the velocity field.

[0018] In some embodiments, the concentration field of the proppant within the hydraulic fracture can be the volume fraction distribution of proppant particles at various spatial locations within the fracture. After the proppant is injected into the fracture with fracturing fluid, it migrates and settles under the influence of flow and gravity, gradually accumulating at the bottom of the fracture to form a sandbar. The spatial distribution of the concentration field directly determines the conductivity distribution after the fracture is closed.

[0019] In some embodiments, the velocity field of the proppant within a hydraulic fracture can be the vector distribution of the transport velocity of proppant particles at various spatial locations within the hydraulic fracture. The velocity field mainly includes the horizontal velocity along the fracture length and the vertical velocity along the fracture height, wherein the vertical velocity is affected by gravitational settling and is related to particle size and fluid viscosity. The velocity field is an important driving factor for the evolution of the concentration field.

[0020] In some embodiments, a field can be a discretized numerical distribution of a physical quantity at various spatial locations within the computational domain. For hydraulic fractures, the fracture computational domain can be divided into several grid cells, with each grid cell or grid node storing a numerical value of the physical quantity. The values ​​at all grid points together constitute the field of that physical quantity. For example, the concentration field is the set of proppant volume fraction values ​​at each grid point, and the pressure field is the set of fluid pressure values ​​at each grid point. When using neural network methods, the aforementioned physical quantity field can be organized into a matrix or tensor form, serving as the input or output of the network.

[0021] This specification provides an embodiment of a proppant transport simulation method based on a physical coding neural network, referring to... Figure 1 As shown, the specific implementation may include the following steps: S101: Input the concentration field of the proppant and the physical property field of the fluid in the hydraulic fracture at the current moment into the pressure inference network to obtain the pressure field of the fluid at the current moment; the pressure inference network includes a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the degree of deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy the preset isobaric boundary conditions and / or zero flux boundary conditions.

[0022] In some embodiments, the pressure inference network may include a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the degree of deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy preset isobaric boundary conditions and / or zero flux boundary conditions.

[0023] In existing methods for solving pressure fields using neural networks, physical constraints are typically incorporated into the loss function as soft constraints. This means that only a penalty term is used to encourage the network output to satisfy the physical equations and boundary conditions, without directly imposing constraints on the network's forward propagation process. This can lead to the pressure field deviating from physical laws during the inference stage. Therefore, in this embodiment, the pressure inference network incorporates both a Poisson regularization term and a boundary mask term.

[0024] The Poisson regularization term characterizes the degree to which the pressure field output by the pressure inference network deviates from the elliptic pressure Poisson equation. During the training phase, this deviation guides the update direction of the network weights, gradually bringing the network output closer to the true pressure field that satisfies the equation. The elliptic pressure Poisson equation describes the pressure distribution of incompressible fluids in a non-uniform viscosity field. The smaller the deviation of the pressure field from this equation, the more physically plausible the pressure field is.

[0025] Boundary mask terms are used to characterize the enforced constraints on the values ​​of grid points in the boundary regions of the pressure field output by the pressure inference network. The boundary conditions for hydraulic fractures include a constant pressure boundary condition where the pressure at the fracture tip equals the preset formation pressure, and a zero flux boundary condition where the fluid normal velocity at the fracture wall is zero. The boundary mask terms do not rely on backpropagation during network training; instead, they directly replace the values ​​of the corresponding boundary grid points in the network output, ensuring that the pressure field strictly satisfies the aforementioned boundary conditions at the output.

[0026] This pressure inference network can employ an encoder-decoder architecture. The encoder is used for stepwise downsampling to expand the receptive field and extract the long-range influence of the spatial distribution of the local field on the global pressure balance. The decoder is used for stepwise upsampling to restore spatial resolution and, in conjunction with skip connections, transmits shallow local detail information to deeper layers to preserve local pressure features near the wall region and injection end. Information is transmitted between the encoder and decoder through feature maps. The feature map output from each layer of the encoder is concatenated with the feature map output from the corresponding layer of the decoder via skip connections in the channel dimension before being input into the next layer of the decoder.

[0027] Through the above mechanism, this pressure inference network can output a pressure field that satisfies the physical equations and boundary constraints in a single forward propagation without iteratively solving a large sparse linear system of equations during the inference phase, thereby improving computational efficiency and maintaining physical consistency.

[0028] In some embodiments, the concentration field of the proppant within the hydraulic fracture at the current moment can be obtained. This concentration field can be input into the constitutive equations to obtain the fluid property field at the current moment. The property field includes the mixing density field and the equivalent viscosity field.

[0029] In existing methods, neural networks typically map directly from input to output, treating physical property parameters either as constants or implicitly within the network's black-box mapping. This lack of explicit constitutive mapping between the concentration field and the physical property field makes it difficult to guarantee physical consistency between them. Therefore, in this embodiment, after obtaining the concentration field of the proppant within the hydraulic fracture at the current moment, it does not directly input the concentration field into the pressure inference network. Instead, it first inputs it into the physical property constitutive equation for explicit calculation.

[0030] Constitutive equations characterize the mapping relationship between proppant concentration and the physical properties of the mixed fluid, specifically including the mixing density equation and the equivalent viscosity equation. The mixing density increases linearly with concentration, reflecting the physical law that proppant enrichment leads to an increase in the mass of local fluid particles, thereby enhancing the effect of gravitational volume forces. The equivalent viscosity increases non-linearly with increasing concentration, approaching infinity when the local concentration approaches the theoretical saturation volume fraction of the proppant. This property characterizes the physical phenomenon of a sharp increase in flow resistance when proppant becomes trapped in cracks.

[0031] By substituting the concentration field at the current moment into the constitutive equations of the physical properties point by point, the physical property field of the fluid at the current moment can be obtained. This physical property field includes the mixing density field and the equivalent viscosity field, which are represented by spatial distribution matrices to characterize the density and viscosity values ​​of the mixed fluid at various locations throughout the entire domain.

[0032] The advantage of this explicit calculation is that the physical field strictly follows the physical laws described by the constitutive equation under any working condition, and will not produce non-physical density or viscosity prediction values ​​due to insufficient training data coverage, and the calculation process does not require additional training resources.

[0033] In some embodiments, the mixed density field and equivalent viscosity field at the current moment can be input into the encoder of the pressure inference network. The encoder extracts multi-scale global features of the input physical property field through a stepwise downsampling operation. The multi-scale global features are used to characterize the long-range influence of the spatial distribution of the physical property field on the global pressure balance.

[0034] The pressure Poisson equation is an elliptic equation, and its solution exhibits global coupling characteristics—a change in physical properties at any point within the crack will instantaneously affect the entire pressure field distribution. If the encoder's receptive field is insufficient to cover the entire area, the network cannot accurately capture the long-range propagation characteristics of pressure. Therefore, the encoder in this embodiment extracts multi-scale global features of the input physical property field through a step-by-step downsampling operation.

[0035] The encoder consists of multiple downsampling modules connected in series. Each downsampling module includes a convolutional layer for extracting local features and a pooling layer for compressing spatial dimensions. The mixed density field and equivalent viscosity field at the current time step are input to the encoder as a dual-channel input tensor, with spatial dimensions equal to the number of grid rows and columns in the crack computation domain. Each downsampling module performs convolution and pooling operations on this tensor sequentially: the convolution operation extracts local spatial features of the physical property field at this scale using multiple trainable convolutional kernels and outputs a multi-channel feature map; the pooling operation reduces the spatial dimension of the feature map by half by taking the maximum or average value within a local window.

[0036] After multi-layer downsampling, the shallow layer output retains high-resolution details of local variations in the physical property field, while the deep layer output encodes low-resolution global features of the global physical property distribution. The deepest layer's feature map has the smallest spatial dimension and its receptive field covers the entire computational domain. The value of each grid node in this layer incorporates the global physical property field distribution information, characterizing the long-range influence of the spatial distribution of the physical property field on the global pressure balance. This long-range influence feature is then passed to the decoder to reconstruct the complete pressure field.

[0037] In some embodiments, the deepest feature map output by the encoder can be input into the decoder of the stress inference network. The decoder recovers the spatial resolution through stepwise upsampling operations and skip connections to obtain the candidate stress field. Skip connections are used to concatenate the shallow spatial features output by each layer of the encoder with the deep semantic features output by the corresponding layer of the decoder in the channel dimension to preserve the local details of the physical field.

[0038] The deepest feature map output by the encoder has low spatial resolution. If it is directly mapped to a pressure field, detailed pressure information in local areas such as the injection end and the wall will be lost, leading to inaccurate calculation of pressure gradients near the wall and injection end. To address this, the decoder can restore spatial resolution by performing stepwise upsampling operations combined with skip connections.

[0039] The decoder may include multiple upsampling modules, the number of which is equal to the number of downsampling modules in the encoder. Each upsampling module receives shallow feature maps passed from skip connections of the corresponding layer in the encoder and deep semantic feature maps output by the previous layer's upsampling module. The shallow feature maps output by each layer of the encoder preserve local details of the physical field, including the location of viscosity abrupt changes near the wall and density gradients near the injection end.

[0040] At each stage, the decoder first upsamples the feature map output from the previous module (e.g., through transposed convolution or interpolation) to double the spatial dimension. Then, through skip connections, it concatenates the corresponding shallow feature map from the encoder with the upsampled deep feature map along the channel dimension. Finally, it performs convolution on the concatenated tensor to fuse shallow detail information and deep semantic information. After successive upsampling and skip connections, the decoder outputs a feature map with the same spatial dimension as the input physical field. This feature map contains both global pressure distribution trends and local boundary detail features, serving as the candidate pressure field output.

[0041] In some embodiments, candidate pressure fields can be input into the boundary mask layer of the pressure inference network. The boundary mask layer directly assigns values ​​to the candidate pressure fields according to a preset boundary mask matrix, replaces the pressure values ​​of the grid points in the candidate pressure fields that are the boundaries of the boundary mask matrix with preset boundary condition values, and outputs the replaced candidate pressure fields as the pressure fields of the fluid at the current moment.

[0042] Although the candidate pressure field output by the decoder satisfies the Poisson equation constraint in the internal region, its boundary values ​​may still deviate from the preset engineering boundary conditions. If only the soft constraints in the loss function are relied upon to approximate the boundary conditions, the boundary accuracy during the inference stage is difficult to guarantee. Therefore, this embodiment sets a boundary mask layer to directly assign and correct the candidate pressure field.

[0043] The boundary mask layer pre-stores a boundary mask matrix with the same dimensions as the computational domain space. This matrix takes values ​​of 0 or 1 at each grid location, where a value of 1 corresponds to the boundary region of the hydraulic fracture (including the tip boundary and wall boundary), and a value of 0 corresponds to the interior region. The boundary mask layer reads this matrix and iterates through all grid points of the candidate pressure field: when the mask value corresponding to a grid point is 1, the pressure value of that grid point is replaced with a preset boundary condition value; when the mask value is 0, the pressure value of that grid point remains unchanged.

[0044] The preset boundary condition values ​​include a constant pressure boundary value at the hydraulic fracture tip and a zero flux boundary value at the fracture wall. The specific value of the constant pressure boundary value at the tip is equal to the preset formation pressure. The implementation of the zero flux boundary condition at the wall includes, but is not limited to, setting the pressure gradient along the wall normal direction at the wall to zero, which is equivalent to setting the pressure value of the grid point at the wall to be equal to the pressure value of its adjacent internal grid point, so as to ensure that the fluid normal velocity at the wall is zero.

[0045] After the direct assignment described above, the values ​​of all boundary grid points in the candidate pressure field are forcibly corrected to satisfy the engineering boundary conditions. This corrected pressure field is then used as the pressure field output of the fluid at the current moment. This boundary mask layer participates in the forward computation during both the training and inference phases of the pressure inference network. During the training phase, the trainable parameters of the encoder and decoder are updated through gradient backpropagation, while during the inference phase, the aforementioned assignment operation is directly executed.

[0046] S102: Calculate the velocity field of the proppant at the current moment based on the pressure field.

[0047] In some embodiments, the physical property field may include a mixed density field and an equivalent viscosity field.

[0048] In existing methods, velocity field calculations typically assume uniform and constant fluid properties, neglecting the spatial differences in density and viscosity caused by uneven proppant concentration distribution, leading to deviations from actual physical conditions. Therefore, in this embodiment, the property field specifically includes two independent spatial distribution matrices: a mixed density field and an equivalent viscosity field.

[0049] The mixing density field characterizes the density distribution of the mixed fluid at all grid locations across the entire domain. It is calculated point-by-point from the current concentration field using the mixing density equation, and its physical meaning is the volume-weighted average density value of the proppant and fracturing fluid within each grid cell. The equivalent viscosity field characterizes the apparent viscosity distribution of the mixed fluid at all grid locations across the entire domain. It is calculated point-by-point from the current concentration field using the equivalent viscosity equation, and its physical meaning is the equivalent viscosity value of the mixed proppant and fracturing fluid within each grid cell. Both property fields are two-dimensional matrices with the same spatial dimension as the concentration field. The matrix elements correspond one-to-one with each grid location in the fracture computational domain, jointly providing property input for the velocity field calculation.

[0050] In some embodiments, the horizontal pressure gradient and the vertical pressure gradient can be calculated based on the pressure field.

[0051] The velocity field is driven by the pressure gradient, and its calculation accuracy directly depends on the accuracy of the pressure gradient. If the gradient is directly calculated using numerical difference on the discrete pressure field, numerical truncation errors may lead to inaccurate gradient calculations in regions where pressure changes are gradual. Therefore, this embodiment calculates the horizontal and vertical pressure gradients based on the pressure field at the current moment.

[0052] The pressure field is stored in matrix form, with row indices corresponding to vertical spatial positions and column indices corresponding to horizontal spatial positions. The horizontal pressure gradient is calculated, but is not limited to: for each row of the pressure field matrix, dividing the pressure difference between adjacent column elements by the corresponding horizontal spatial step size yields the horizontal pressure gradient value at each column position in that row. Similarly, the vertical pressure gradient is calculated, but is not limited to: for each column of the pressure field matrix, dividing the pressure difference between adjacent row elements by the corresponding vertical spatial step size yields the vertical pressure gradient value at each row position in that column.

[0053] In terms of calculation methods, the horizontal pressure gradient can be implemented using one of three discretization schemes, including but not limited to forward differencing, backward differencing, or central differencing. The vertical pressure gradient is calculated similarly and will not be elaborated further here.

[0054] In some embodiments, the horizontal velocity field can be calculated based on the equivalent viscosity field and horizontal pressure gradient at the current moment; the horizontal velocity field is used to characterize the proppant migration velocity along the length of the hydraulic fracture.

[0055] Horizontal velocity is driven by the horizontal pressure gradient and is inversely proportional to the fluid's equivalent viscosity—the higher the viscosity, the lower the velocity produced by the same pressure gradient. To accurately calculate the horizontal velocity field, this embodiment calculates the horizontal velocity field based on the current equivalent viscosity field and the horizontal pressure gradient.

[0056] The calculation of the horizontal velocity field is based on the parallel plate laminar flow assumption. Under this assumption, the crack width is much smaller than the length, and the fluid flow in the crack can be regarded as laminar motion between two parallel plates. The horizontal velocity is parabolic along the crack width direction, and the average horizontal velocity along the width direction is proportional to the pressure gradient and inversely proportional to the fluid viscosity.

[0057] In one implementation, the horizontal pressure gradient matrix, equivalent viscosity field matrix, and crack width constant can be substituted point by point into a preset formula to obtain the horizontal velocity field matrix. In another implementation, a differential convolution operator can be used for calculation. Specifically, a spatial partial differential operation is performed on the pressure field using a fixed-weight differential convolution kernel to obtain the horizontal pressure gradient matrix, which is then divided point by point by the equivalent viscosity field matrix and multiplied by a constant coefficient. The former is suitable for scenarios with high computational speed requirements, while the latter is suitable for scenarios requiring unified processing of multiple differential operations within the same computational framework.

[0058] The horizontal velocity field is stored in the form of a two-dimensional matrix to characterize the distribution of proppant migration velocity along the length of the hydraulic fracture (i.e., the horizontal direction) at all grid locations in the global domain.

[0059] In some embodiments, the base vertical velocity field can be calculated based on the mixed density field and vertical pressure gradient at the current moment.

[0060] The vertical velocity is driven not only by the vertical pressure gradient but also by gravity—a denser fluid mixture tends to move downwards under gravity. Ignoring gravity would underestimate the driving factors of proppant settling. Therefore, this embodiment calculates the foundation's vertical velocity field based on the current mixing density field and vertical pressure gradient.

[0061] The calculation of the basic vertical velocity field is also based on the assumption of laminar flow on a parallel flat plate, but the contribution of the gravity term needs to be superimposed on the pressure gradient term.

[0062] In one embodiment, the vertical pressure gradient matrix, the mixed density field matrix, and the equivalent viscosity field matrix can be substituted point by point into a preset formula to obtain the basic vertical velocity field. In another embodiment, the vertical pressure gradient can be calculated using a differential convolution operator, and then multiplied point by point by the mixed density field matrix with a gravity term, followed by division by the equivalent viscosity field matrix.

[0063] The basic vertical velocity field is stored in the form of a two-dimensional matrix, as one of the components of the vertical velocity field.

[0064] In some embodiments, the sedimentation correction velocity field can be calculated based on the concentration field at the current moment.

[0065] Proppant particles exhibit slip-sedimentation relative to the mixed fluid. This sedimentation behavior conforms to classical free sedimentation laws in pure fracturing fluid zones, but it is significantly attenuated in high-concentration sandbar accumulation zones due to interparticle interactions. Using a constant sedimentation velocity would overestimate the sedimentation rate in the sandbar zone. Therefore, this embodiment calculates a sedimentation correction velocity field based on the concentration field at the current moment.

[0066] A concentration masking factor can be added to the classic Stokes free settling formula. This factor causes the settling correction rate to decrease monotonically with increasing concentration, and approaches zero when the local concentration approaches 1, characterizing the physical phenomenon of impeded particle settling in high-concentration sand embankment areas. When the local concentration approaches zero, the concentration masking factor approaches 1, degenerating into the classic single-particle free settling formula.

[0067] In one implementation, the concentration field, equivalent viscosity field, and various constants can be substituted point by point into a preset formula to obtain the sedimentation correction velocity field. In another implementation, a Π-type convolutional layer can be used to perform multi-order polynomial product operations on the concentration field and equivalent viscosity field to approximate the value of the sedimentation correction velocity field. This is suitable for scenarios where sedimentation behavior has additional complex influencing factors and is difficult to describe accurately with a single analytical formula.

[0068] The settlement correction velocity field is stored in the form of a two-dimensional matrix and serves as a correction term for the vertical velocity field.

[0069] In some embodiments, the vertical velocity field can be determined based on the basic vertical velocity field and the settlement-corrected velocity field; the vertical velocity field is used to characterize the proppant migration velocity in the crack height direction.

[0070] The true vertical velocity is not solely determined by the pressure gradient and gravity-driven base flow; it also requires the superposition of the slip and settling velocities of the proppant particles relative to the fluid. If the base vertical velocity field is directly used as the final estimate of the vertical velocity field, the impact of impeded settlement in the sand embankment area on concentration redistribution cannot be accurately depicted, leading to distorted sand embankment morphology predictions. Therefore, this embodiment determines the vertical velocity field based on the base vertical velocity field and the settlement correction velocity field.

[0071] The methods for determining the vertical velocity field include, but are not limited to: adding the basic vertical velocity field and the settlement correction velocity field point by point to obtain the vertical velocity value at each grid position in the vertical velocity field. The meaning of this superposition operation is that the total velocity of the fluid in the vertical direction is composed of two parts: the basic flow velocity driven by the pressure gradient and gravity (calculated by the laminar flow formula of parallel plate) and the slip settlement velocity of the proppant particles relative to the fluid (calculated by the hindered settlement formula).

[0072] At the operational level, this point-by-point addition is an element-by-element addition operation of two matrices of the same dimension, with the output matrix having the same spatial dimension as the input matrix. The sign of each element in the vertical velocity field indicates the velocity direction—positive values ​​are upward, and negative values ​​are downward. In the gravity-dominated region (lower part of the crack), the vertical velocity field is dominated by downward settling, consistent with the physical process of sand embankment accumulation from bottom to top in the concentration field.

[0073] The vertical velocity field is stored in the form of a two-dimensional matrix to characterize the distribution of proppant migration velocity in the hydraulic fracture height direction (i.e., the vertical direction) at all grid locations across the entire domain.

[0074] In some embodiments, the velocity field at the current moment can be determined based on the horizontal velocity field and the vertical velocity field.

[0075] The horizontal and vertical velocity fields characterize the proppant's migration velocities in two orthogonal directions. If these are treated as independent variables, the vector characteristics of the velocity field—i.e., directional information—cannot be fully represented. Therefore, this embodiment determines the velocity field at the current moment based on the horizontal and vertical velocity fields.

[0076] The velocity field can be determined in ways including but not limited to: combining the horizontal and vertical velocity fields along the channel dimension to form a two-dimensional velocity field tensor with two channels, where the first channel is the horizontal velocity field matrix and the second channel is the vertical velocity field matrix. At each grid position in the spatial dimension, the two elements of this two-channel tensor correspond to the horizontal and vertical velocity components at that position, respectively, and together they constitute the velocity vector at that position.

[0077] In one implementation, the two single-channel matrices can be merged into a single two-channel tensor. In another implementation, the horizontal and vertical velocity fields can be stored as two independent two-dimensional matrices and retrieved separately using indices in subsequent calculations.

[0078] The velocity field is stored in the form of a two-channel tensor, with each spatial location containing both horizontal and vertical velocity components. This is used to calculate the concentration flux in the concentration convection conservation equation and drive the update of the concentration field at the next moment.

[0079] S103: Input the current concentration field and velocity field into the concentration inference network to obtain the concentration field and property field at the next moment; the concentration inference network is coupled with multiple Π-type convolutional modules, the Π-type convolutional modules include differential convolutional layers and Π-type convolutional layers. The differential convolutional layers are used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolutional layers are used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution. The rheological term is used to characterize the rheological feedback effect of proppant viscosity change and particle interaction on the concentration distribution.

[0080] In some embodiments, the concentration inference network is coupled with multiple Π-type convolutional modules. Each Π-type convolutional module includes a differential convolutional layer and a Π-type convolutional layer. The differential convolutional layer is used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolutional layer is used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution, and the rheological term is used to characterize the rheological feedback effect of proppant viscosity changes and particle interactions on the concentration distribution.

[0081] Existing physical information neural networks, when dealing with convection-dominated problems, typically assign all terms of the physical equations to the network for fitting, without distinguishing between precisely expressible linear differential terms and complex nonlinear terms. This forces the network to learn two different mapping relationships simultaneously within the same set of parameters, resulting in significant training difficulties and a challenge in balancing the accuracy of both. Therefore, the concentration inference network in this embodiment couples multiple Π-shaped convolutional modules.

[0082] The Π-type convolutional module is used to characterize the global estimation of the spatial term in the concentration convection conservation equation. A single Π-type convolutional module consists of two parallel branches: a differential convolutional layer and a Π-type convolutional layer. The two branches process the same input separately and then sum the outputs. The differential convolutional layer is used to characterize the deterministic computational branch that performs spatial partial differential operations based on fixed differential weights to solve for the estimation of the convection term in the concentration convection conservation equation. The Π-type convolutional layer is used to characterize the trainable approximation branch that performs multi-order polynomial product operations to fit the estimation of the rheological term in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution, i.e., the concentration flux change caused by the velocity field carrying proppant particles through space. The rheological term is used to characterize the rheological feedback effect of proppant viscosity changes and particle interactions on the concentration distribution, i.e., the nonlinear coupling effect of the concentration field on its own evolution by changing the viscosity of the mixed fluid and the interparticle forces.

[0083] The sum of the outputs of the two branches yields the overall estimate of the spatial term in the concentration convection conservation equation by this Π-shaped convolutional module. When multiple Π-shaped convolutional modules are connected in series, the output of the previous module serves as the input of the next module, progressively extracting higher-order nonlinear features.

[0084] In some embodiments, the concentration field and velocity field at the current moment can be input into the differential convolutional layer in the concentration inference network to obtain a first estimate of the convection term in the concentration convection conservation equation.

[0085] The convection term in the concentration convection conservation equation has a clear mathematical form: the dot product of the velocity field and the concentration gradient. The concentration gradient is the spatial partial derivative of the concentration field, and its calculation is essentially performing a spatial differentiation operation on the concentration field matrix, rather than approximating an unknown function. If this calculation process is delegated to a trainable network for data-driven learning, unnecessary fitting errors will be introduced, and it may gradually deviate from the true physical values ​​in long-term time-series extrapolations. Therefore, this embodiment inputs the current concentration field and velocity field into the differential convolutional layer of the concentration inference network to obtain the first estimate of the convection term in the concentration convection conservation equation.

[0086] The differential convolutional layer contains a differential convolution kernel with fixed weights. These fixed weights are the discretization coefficient matrix of the spatial partial differential operators in the convection term of the concentration convection conservation equation. The values ​​of each element in this discretization coefficient matrix are determined by the corresponding spatial partial differential operators (such as...). / x or / The mathematical form of y is uniquely determined and remains constant after the differential convolutional layer is constructed, and does not participate in network training. The differential convolutional layer uses the fixed-weight differential convolutional kernel as a sliding window template, sliding it row by row and column by column on the input concentration field matrix, performing a weighted summation operation at each window position, and outputting the spatial partial derivative value of the concentration field at that position.

[0087] After performing spatial partial differential operations on the concentration field to obtain the concentration gradient, the concentration gradient matrix is ​​multiplied point-by-point with the velocity field matrix to obtain the first estimate of the convection term in the concentration convection conservation equation. This first estimate has the same spatial dimension as the input concentration field.

[0088] In some embodiments, the concentration field and velocity field at the current moment can be input into the Π-type convolutional layer in the concentration inference network to obtain a second estimate of the rheological term in the concentration convection conservation equation.

[0089] The rheological term in the concentration convection conservation equation is a nonlinear function of the concentration and velocity fields, encompassing complex physical processes involving nonlinear changes in equivalent viscosity with concentration, multi-particle interactions, and obstructed sedimentation correction. Its specific functional form is difficult to derive precisely from pure theory. Attempting to calculate this term accurately using a fixed differential operator would be impossible due to the lack of an explicit mathematical expression. Therefore, this embodiment inputs the current concentration and velocity fields into a Π-shaped convolutional layer in the concentration inference network to obtain a second estimate of the rheological term in the concentration convection conservation equation.

[0090] A Π-type convolutional layer can include multiple parallel trainable convolutional layers, element-wise multiplication layers, and combined convolutional layers. Multiple parallel trainable convolutional layers extract features from the input tensor at different scales, outputting multiple feature maps. The element-wise multiplication layer multiplies the values ​​of these multiple feature maps at the same spatial location point by point, obtaining a product feature map. The combined convolutional layer performs a convolution operation on this product feature map, fusing multi-channel information into a single-channel output. This sequence of parallel extraction, element-wise multiplication, and output fusion is mathematically equivalent to a high-order polynomial combination of the input features, giving it the ability to approximate any continuous nonlinear function. It can be used to fit the nonlinear part of the rheological term that cannot be precisely expressed by a fixed differential operator.

[0091] During the training process of the Π-type convolutional layer, its trainable parameters are updated by minimizing the physical residuals of the concentration convection conservation equation, so that the output of the Π-type convolutional layer gradually approaches the true rheological term value.

[0092] In some embodiments, the time rate of change of the concentration field at the current moment can be calculated based on the first estimate and the second estimate.

[0093] The concentration convection conservation equation relates the temporal rate of change of concentration to the spatial term (the sum of the convection and rheological terms)—the temporal rate of change equals the overall estimate of the spatial term. If the convection and rheological terms are calculated separately and then not fused, a complete estimate of the spatial term cannot be obtained, thus hindering the understanding of the temporal evolution of the concentration field. Therefore, this embodiment calculates the temporal rate of change of the concentration field at the current moment based on the first and second estimates.

[0094] The calculation of the rate of change over time includes, but is not limited to: adding the first estimate of the convection term output by the differential convolutional layer and the second estimate of the rheological term output by the Π-type convolutional layer point by point along the channel dimension to obtain the overall estimate of the spatial term in the concentration convection conservation equation. The addition operation means: element-wise addition of two matrices of the same dimension, where each element of the output matrix is ​​equal to the sum of the element values ​​of the first estimate matrix and the second estimate matrix at the same spatial location.

[0095] The overall estimate of the spatial term is used as the rate of change of the concentration field at the current moment. This rate of change has the same spatial dimension as the concentration field at the current moment, and the value at each grid location represents the instantaneous rate of change of the proppant concentration at that location at the current moment—positive values ​​indicate that the concentration increases with time, and negative values ​​indicate that the concentration decreases with time.

[0096] In some embodiments, the concentration field increment from the current moment to the next moment can be calculated based on the rate of change over time and a preset time step to determine the concentration field at the next moment.

[0097] After obtaining the rate of change of the concentration field at the current moment, the evolution of the concentration field from the current moment to the next moment needs to be realized through time discretization. The key to this discretization process is to determine the concentration field increment. To this end, this embodiment calculates the concentration field increment from the current moment to the next moment based on the rate of change of time and a preset time step, so as to determine the concentration field at the next moment.

[0098] The calculation of concentration field increments includes, but is not limited to, multiplying the time rate of change matrix by a preset time step to obtain the concentration field increment matrix. This multiplication is an element-wise multiplication of a scalar and a matrix—the preset time step is a scalar constant, which is multiplied by each element in the time rate of change matrix to obtain the corresponding element in the concentration field increment matrix. This concentration field increment matrix represents the change in concentration at each grid location during the time interval from the current time to the next time step. The value of the concentration field increment can be positive or negative—a positive increment indicates that the concentration at that location will increase in the next time step, and a negative increment indicates that the concentration will decrease in the next time step.

[0099] The determination of the concentration field at the next time step includes, but is not limited to, adding the current concentration field matrix to the concentration field increment matrix point by point to obtain the concentration field matrix at the next time step. This point-by-point addition is an element-by-element addition of matrices of the same dimension, and each element of the output matrix is ​​equal to the sum of the element value at the corresponding position of the current concentration field and the element value at the corresponding position of the concentration field increment. The value at each grid position in the resulting concentration field matrix at the next time step is the proppant volume concentration at that position in the next time step.

[0100] In some embodiments, the concentration field at the next moment can be saturated and clamped based on a preset concentration threshold.

[0101] In long-term iterative simulations, accumulated numerical errors may cause concentration values ​​at some grid points to exceed physically permissible ranges—for example, below 0 (negative concentration) or above the theoretical maximum accumulation concentration (above 1). If these non-physical concentration values ​​are not corrected in time, they will spread throughout the computational domain via convective transport in subsequent iterations, causing the entire simulation results to fail. Therefore, this embodiment uses a preset concentration threshold to saturate and clamp the concentration field at the next time step.

[0102] The preset concentration thresholds include a lower concentration limit and an upper concentration limit. The lower concentration limit is usually set to 0, representing the lowest physical limit at which no proppant particles exist at that grid location. The upper concentration limit is usually set to the theoretical saturation volume fraction of the proppant, c_max, representing the highest physical limit at which no more proppant can be accommodated in a tightly packed state at that grid location.

[0103] The specific operation of saturation clamping includes, but is not limited to: traversing all grid positions of the concentration field matrix at the next time step, judging and correcting the concentration value at each position. When the concentration value at a certain position is lower than the lower concentration threshold, the concentration value at that position is replaced with the lower concentration threshold. When the concentration value at a certain position is higher than the upper concentration threshold, the concentration value at that position is replaced with the upper concentration threshold. When the concentration value at a certain position is within the range between the lower and upper concentration thresholds, the concentration value at that position remains unchanged.

[0104] In one implementation, values ​​exceeding the boundary can be directly set as boundary values. In another implementation, values ​​exceeding the boundary can be compressed to the vicinity of the boundary by a certain proportion instead of being directly truncated, to avoid numerical abrupt changes introduced by the clamping operation affecting the physical evolution of adjacent regions. After saturation clamping, the values ​​at each grid position of the concentration field are all within the physically permissible range.

[0105] In some embodiments, the concentration field at the next moment after clamping can be input into the fluid's property constitutive equation to obtain the property field at the next moment; the property constitutive equation is used to characterize the mapping relationship between the fluid's mixing density and equivalent viscosity and the proppant concentration.

[0106] The concentration field itself is a physical quantity describing the distribution of the proppant. Changes in the proppant distribution ultimately affect the pressure and velocity fields by altering the properties of the mixed fluid, thus forming a complete physical feedback loop. If only the concentration field is updated without updating the corresponding property field, the property data used for pressure calculation at the next time step will not reflect the current state of the concentration field, resulting in a lag in the pressure field solution. Therefore, this embodiment inputs the clamped concentration field for the next time step into the fluid's constitutive equations to obtain the property field for the next time step.

[0107] Constitutive equations of physical properties are used to characterize the mapping relationship between the mixing density and equivalent viscosity of a fluid and the proppant concentration.

[0108] The methods for obtaining the physical property field include, but are not limited to: applying the mapping relationship point by point to all grid positions of the concentration field matrix at the next time step, calculating the mixing density value and equivalent viscosity value at each position, and forming the mixing density field matrix and the equivalent viscosity field matrix, respectively. These two matrices together constitute the physical property field at the next time step.

[0109] In one implementation, numerical calculations can be performed independently for each grid location based on the property constitutive equation. The advantage of this approach is that the calculation process is completely transparent and does not introduce any additional approximation errors. In another implementation, the mapping relationship between concentration and property can be learned through a small neural network (such as a single-hidden-layer fully connected network). The advantage of this approach is that it has stronger expressive power when the property constitutive equation has a complex form or when multiphase interface effects exist.

[0110] This property field was then used to solve for the pressure in the next iteration, ensuring that the property field remained physically consistent with the current concentration field.

[0111] S104: Take the next moment as the new current moment, iteratively execute the above steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field until the preset simulation time is reached, and obtain the simulation results of the proppant concentration field at each moment.

[0112] In some embodiments, the next moment can be taken as the new current moment, and the steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field can be executed iteratively until the preset simulation time is reached to obtain the simulation results of the proppant concentration field at each moment.

[0113] A single time step only yields the concentration field at the next moment, while engineering practice requires spatiotemporal evolution data of the proppant from injection to final settling, and the simulation process must have an automatic termination condition. Without setting an iteration termination check, the simulation will continue indefinitely; without saving the results after each iteration, only the field data at the last moment can be obtained, failing to present the complete migration process. Therefore, this embodiment uses the next moment as the new current moment, iteratively executing the steps of pressure field solution, velocity field calculation, and concentration and property field updates until the preset simulation duration is reached, obtaining the proppant concentration field simulation results at each moment.

[0114] The specific implementation of the iterative deduction includes, but is not limited to: after completing one update of the concentration field and property field from the current time to the next time, taking the next time as the new current time, and repeating the complete sequence of steps: updating property field, solving the pressure field with the pressure inference network, reconstructing the velocity field, and updating the concentration field with the concentration inference network. After each iteration, the current time is advanced forward by a preset time step.

[0115] The preset simulation duration is a pre-defined total simulation time, used to characterize the physical time range that the proppant migration simulation needs to cover. After each iteration, the currently accumulated time steps are compared with the preset simulation duration. If the physical time corresponding to the current accumulated time steps multiplied by the preset time step is less than the preset simulation duration, the next iteration continues. If the physical time corresponding to the current accumulated time steps multiplied by the preset time step reaches or exceeds the preset simulation duration, the iterative simulation terminates.

[0116] During the iteration process, the concentration field, pressure field, velocity field, and physical property field at each moment are stored, forming a multiphysics sequence arranged in chronological order. The simulation results of the proppant concentration field at each moment are output in time series form to characterize the spatiotemporal evolution data of the proppant in the hydraulic fracture from injection to final settling.

[0117] In one implementation, the total number of iterations is calculated by dividing a preset simulation duration by a preset time step before the simulation begins. The simulation terminates when the number of iterations reaches the total number of iterations. In another implementation, the accumulated physical time is calculated after each iteration, and the simulation terminates when the accumulated physical time reaches the preset simulation duration.

[0118] In some embodiments, the concentration field, pressure field, and velocity field at each time point can be inversely normalized to obtain concentration field cloud maps, pressure field cloud maps, and velocity field cloud maps at each time point, which are used to characterize the sandbank accumulation morphology and concentration front migration position, pressure gradient distribution characteristics, and fluid flow direction of the proppant in the hydraulic fracture space.

[0119] Concentration inference networks and pressure inference networks typically normalize input and output data to accelerate training convergence. Therefore, the numerical ranges of the physical quantities output by the network are compressed to specific intervals (e.g., [0,1] or [-1,1]), not directly corresponding to the actual physical dimensions. Directly using these normalized values ​​for analysis or visualization will not yield concentration, pressure, and velocity values ​​with real physical meaning. Therefore, this embodiment performs inverse normalization on the concentration, pressure, and velocity fields at each time point to obtain concentration field contour maps, pressure field contour maps, and velocity field contour maps at each time point.

[0120] The specific implementation of the inverse normalization process includes, but is not limited to: obtaining scaling parameters used to normalize each physical quantity during the data preprocessing stage. These scaling parameters include the minimum, maximum (for min-max normalization), or mean and standard deviation (for z-score normalization) of each physical quantity. Substituting each element value in the concentration field matrix at each time step into the corresponding inverse normalization transformation formula, the normalized concentration values ​​are restored to the true proppant volume fraction values. The same inverse normalization operation is performed on the pressure field and velocity field matrices, respectively restoring them to the true physical values ​​with pressure dimensions (Pascals) and velocity dimensions (m / s).

[0121] After inverse normalization, the concentration, pressure, and velocity fields at each time step are all numerical matrices with actual physical dimensions. When generating a concentration field contour map based on the concentration field matrix at each time step, the concentration value at each grid location is mapped to grayscale or color value. When generating a pressure field contour map based on the pressure field matrix at each time step, the pressure value at each grid location is mapped to grayscale or color value. When generating a velocity field contour map based on the velocity field matrix at each time step, the velocity value at each grid location is mapped to grayscale or color value, or the direction is indicated by arrows.

[0122] Concentration field contour maps are used to characterize the sandbank accumulation morphology and concentration front migration location of proppant within the hydraulic fracture space—the sandbank accumulation morphology is represented by the shape and extent of high-concentration areas, and the concentration front migration location is represented by the advancing front of the concentration isolines. Pressure field contour maps are used to characterize the pressure gradient distribution characteristics—the magnitude and direction of the pressure gradient are represented by the density and orientation of the pressure isolines. Velocity field contour maps are used to characterize the direction of fluid flow—the direction and magnitude of the velocity vector are represented by the direction and length of the arrow.

[0123] In some embodiments, the concentration field at each time point can be compared with a reference concentration field point by point to determine the simulation accuracy.

[0124] The baseline concentration field is a reference value of the concentration field calculated under the same working conditions using a high-precision numerical simulation method (such as the fine-grid finite difference method). Its grid density is higher than the grid density used when generating the simulation samples in this method, and it is used as a reference standard to measure the accuracy of the simulation results of this method. Specific implementation methods for point-by-point comparison include, but are not limited to: for the concentration field matrix at each time step, traversing all grid positions, comparing the concentration value simulated by this method at that position with the concentration value at the corresponding position in the baseline concentration field, and calculating the difference between the two. This difference can be calculated in various forms, including absolute difference and squared difference.

[0125] Based on the point-by-point comparison results, multiple error indices are calculated across the entire domain. The maximum absolute error across the entire domain is the maximum absolute difference between the simulated value and the reference value across all grid locations, used to characterize the degree of extreme deviation in the simulation results. The average absolute error across the entire domain is the average of the absolute differences across all grid locations, used to characterize the overall deviation of the simulation results. The average relative error across the entire domain is the average ratio of the absolute differences across all grid locations to the absolute values ​​of the reference values, used to characterize the relative deviation of the simulation results at different concentration levels. These error indices together constitute a quantitative evaluation system for the accuracy of the simulation results.

[0126] In one implementation, one or more key time steps in the simulation process (such as early injection, mid-migration, and late sedimentation) can be selected, and error calculations can be performed only on the concentration fields at these times. In another implementation, the error of the concentration fields at all output times can be calculated time-by-time, and the time average can be obtained to comprehensively reflect the accuracy performance throughout the entire simulation process.

[0127] In some embodiments, the simulation time can be compared with the time taken by a benchmark numerical simulation method to determine simulation efficiency.

[0128] Simulation time consists of total time and average single-step time. Total time is the total computation time from the start to the end of the simulation, including forward propagation and post-processing time for all iterations, measured in seconds. Average single-step time is the total time divided by the total number of iterations, representing the average computation time for a single time step. Specific methods for recording simulation time include, but are not limited to: recording a system timestamp before the first iteration begins and recording a system timestamp after the last iteration, using the difference between the two timestamps as the total time.

[0129] The benchmark numerical simulation method may include, but is not limited to, the high-precision finite difference method. Its simulation time is obtained by recording the time taken to run the same simulation under the same conditions on a high-performance computing node. Dividing the total time of this method by the total time of the benchmark numerical simulation method yields the speedup ratio—the larger the speedup ratio, the more significant the efficiency improvement of this method.

[0130] In one implementation, the total time of this method is directly compared with that of the benchmark method within a preset simulation duration to verify the single-run efficiency advantage of this method under specific conditions. In another implementation, within the same total simulation time, the number of multi-parameter scenarios that this method and the benchmark method can complete are counted respectively. The cumulative efficiency advantage of this method in high-throughput simulation is demonstrated by comparing the number of completed scenarios—because this method has a shorter single-run simulation time, it can cover more parameter combinations within the same total time, thus making it more suitable for large-scale uncertainty quantification analysis or parameter optimization sweep scenarios.

[0131] In some embodiments, the elliptic pressure Poisson equation can be constructed using the following formula: ; In the formula, For vector differential operators, Fluid pressure, For fluid mixing density, The equivalent viscosity of the fluid. This is the acceleration due to gravity.

[0132] The elliptic pressure Poisson equation describes the pressure distribution of incompressible fluids in non-uniform physical fields. The divergence operator on the left-hand side of the equation... Gradient operator is used to characterize the net flux of the vector field within brackets at various spatial locations, and its result represents the degree of convergence or divergence of the vector field at that location. This equation characterizes the rate and direction of pressure change in various spatial directions. Dividing by the equivalent viscosity indicates that the pressure gradient is less efficient in regions with higher viscosity and more efficient in regions with lower viscosity; the equivalent viscosity acts as a modulation coefficient to regulate the spatial propagation capability of the pressure gradient. The outer divergence sums the components of the viscosity-modulated pressure gradient in each direction, yielding the net driving force at that location. The gradient operator on the right-hand side of the equation characterizes the spatial variation of the product field of mixing density and gravitational acceleration, where the mixing density varies with the proppant concentration distribution, and the gravitational acceleration direction is vertically downward. Dividing by the equivalent viscosity μ and summing the results using the outer divergence yields the effective buoyancy driving force generated by the coupling of the non-uniform density field and gravity. Essentially, this equation represents the balance between the diffusion intensity of pressure under viscosity modulation and the gravity-buoyancy source term driven by density differences; the pressure field is determined by the combined effect of both.

[0133] In some embodiments, the concentration convection conservation equation can be constructed using the following formula: ; ; ; ; In the formula, , and These are the time term, convection term, and rheological term, respectively. Fluid transport velocity, including horizontal transport velocity. and vertical transport speed , The sliding settlement velocity, For vector differential operators, Fluid pressure, For fluid mixing density, The equivalent viscosity of the fluid. It is the acceleration due to gravity. The average width of the hydraulic fracture. This refers to the volume concentration of the proppant. The true density of the proppant particles. The density of pure fracturing fluid. The equivalent spherical particle size of the proppant, For horizontal step size, The vertical step size. For time.

[0134] The concentration convection conservation equation is used to describe the evolution of proppant concentration over time. (Time term) This term characterizes the rate of change of proppant concentration over time at various spatial locations; positive values ​​indicate an increase in concentration, and negative values ​​indicate a decrease. (Convection term) The concentration flux change generated by fluid flow carrying proppant particles through space is characterized by calculations that rely on the dot product of horizontal and vertical transport velocities with the concentration gradient, where the horizontal velocity is determined by the horizontal pressure gradient. The ratio to the equivalent viscosity determines the vertical velocity, which is determined by the vertical pressure gradient. The ratio of the superposition term of the mixture density and gravitational acceleration to the equivalent viscosity is determined, and the slip settling velocity is superimposed. Rheological term. This is used to characterize the feedback effect of complex factors on concentration distribution, such as the nonlinear change of equivalent viscosity with concentration, interparticle interactions, and hindered sedimentation correction, which cannot be explicitly expressed by a fixed differential operator. Among them, the hindered sedimentation velocity includes a concentration masking factor. This causes the settling velocity to decrease as the concentration increases. The overall physical meaning of this equation is that the instantaneous rate of change of proppant concentration at a certain location is jointly determined by the convective contribution of fluid flow carrying particles into and out of the control volume, the nonlinear viscosity change, and the contribution of particle interaction.

[0135] In some embodiments, a first simulation sample including multiple time-varying material property fields and corresponding pressure fields can be obtained. An initial pressure inference network is constructed; the initial pressure inference network includes an encoder, a decoder, and a skip connection layer. The skip connection layer is used to concatenate the shallow spatial features output from each layer of the encoder with the deep semantic features output from the corresponding layer of the decoder. A boundary mask layer is set at the output of the decoder, which is used to constrain the pressure field output by the decoder according to a preset boundary mask matrix. The initial pressure inference network is trained using the sum of the pressure field prediction loss of the first simulation sample and the discrete residual loss of the elliptic pressure Poisson equation as the first loss, until the first loss satisfies the convergence condition, thus obtaining the pressure inference network.

[0136] Existing physical information neural networks typically employ a simple weighted sum of data loss and physical loss during training, followed by backpropagation to update all network parameters. This lack of hard constraints on the network output results in the physical constraints being approximated only statistically, failing to guarantee that the network output accurately meets engineering boundary conditions in the boundary region. Furthermore, the weights of the physical loss are difficult to uniformly determine under different physical property distributions, easily leading to gradient imbalance. Therefore, this embodiment trains the initial stress inference network through the following steps.

[0137] First simulation samples, including material property fields and corresponding pressure fields at multiple time points, are obtained. These first simulation samples provide input data and supervision for training the initial pressure inference network. The material property fields in the first simulation samples include a mixed density field and an equivalent viscosity field, while the pressure field is a numerical solution obtained by solving the elliptic pressure Poisson equation under the same material property field conditions using the high-precision finite difference method. The paired material property field-pressure field data at multiple time points are derived from simulation results at different time steps, covering various possible patterns of spatial distribution of the material property fields, thus ensuring that the training samples cover as wide a range of operating conditions as possible. Each sample in the first simulation samples contains an input tensor and a target tensor. The input tensor is the material property field at a certain time point (a two-channel matrix, corresponding to the mixed density field and the equivalent viscosity field, respectively), and the target tensor is the pressure field at that time point (a single-channel matrix). The material property fields and pressure fields satisfy the physical mapping relationship of the elliptic pressure Poisson equation, and this mapping relationship is learned by the network during training.

[0138] An initial stress inference network is constructed. This network maps the input physical property field to a corresponding stress field, and its architecture is an encoder-decoder structure. The encoder extracts multi-scale global features of the input physical property field through progressive downsampling operations. It consists of multiple downsampling modules cascaded together. Each downsampling module contains a convolutional layer and a pooling layer. The convolutional layer extracts local spatial features, and the pooling layer compresses the spatial dimension to expand the receptive field. The decoder restores spatial resolution through progressive upsampling operations. It consists of multiple upsampling modules cascaded together. Each upsampling module contains an upsampling layer and a convolutional layer. The upsampling layer expands the spatial dimension, and the convolutional layer fuses features and generates the output. Skip connection layers concatenate the shallow spatial features output from each layer of the encoder with the deep semantic features output from the corresponding layer of the decoder along the channel dimension. The shallow output feature map of the encoder has a high spatial resolution, preserving local details of the input physical field (such as the precise spatial location of viscosity abrupt change and the specific distribution of density gradient near the injection end). This information is directly passed to the corresponding layer of the decoder through the skip connection layer. It is then concatenated with the deep semantic feature map of the decoder layer after upsampling in the channel dimension and input into the convolutional layer of the same layer for feature fusion. This allows the decoder to utilize both the shallow local details and the deep global distribution trend information when reconstructing the pressure field.

[0139] A boundary mask layer is set at the output of the decoder. This layer constrains the pressure field output by the decoder according to a preset boundary mask matrix. This boundary mask matrix is ​​a binary matrix with the same spatial dimension as the pressure field, where each element is either 0 or 1. A value of 1 corresponds to a boundary region in the computational domain (including the boundaries of hydraulic crack tips and walls), while a value of 0 corresponds to an interior region. At the end of the forward propagation of the initial pressure inference network, the boundary mask layer performs an assignment operation: for grid points in the pressure field matrix output by the decoder that correspond to a value of 1 in the boundary mask matrix, the pressure value of the grid point is replaced with a preset boundary condition value; for grid points that correspond to a value of 0, the pressure value remains unchanged. This assignment operation, along with other parts of the initial pressure inference network, participates in the forward computation, and the assigned pressure field becomes the final output of the network in the current forward propagation. During training, the boundary mask layer allows gradients to propagate from the boundary region to the interior region. The network learns through backpropagation how to align the pressure field of the interior region with the fixed boundary values, rather than learning the boundary values ​​themselves.

[0140] The initial pressure inference network is trained using the sum of the pressure field prediction loss and the discrete residual loss of the elliptic pressure Poisson equation for the first simulation sample as the first loss, until the first loss meets the convergence condition. The pressure field prediction loss characterizes the deviation between the pressure field output by the initial pressure inference network and the true pressure field value in the first simulation sample. Specifically, it is expressed as mean square error, calculated by averaging the squares of the differences between the predicted and true values ​​at all grid locations. The discrete residual loss of the elliptic pressure Poisson equation characterizes the degree of deviation of the output pressure field from the elliptic pressure Poisson equation. Specifically, it is the sum of the squares of the residuals obtained by substituting the output pressure field and the input physical property field into the elliptic pressure Poisson equation using a conserved finite volume discretization scheme. This scheme is obtained by integrating over each grid control volume, ensuring strict conservation of flow rate between adjacent grids. The first loss is obtained by weighted summing of the pressure field prediction loss and the discrete residual loss of the elliptic pressure Poisson equation. During the gradient backpropagation of the first loss, all trainable parameters of the encoder and decoder in the initial stress inference network are updated. The assignment operation of the boundary mask layer only participates in the forward computation and does not introduce trainable parameters. During training, the numerical change of the first loss on the validation set is monitored in real time. When the first loss no longer continues to decrease (i.e., the validation set loss no longer decreases within several consecutive training cycles) or the decrease is lower than a preset threshold, it is determined that the first loss meets the convergence condition, training is stopped, and the current network parameters are fixed as the parameters of the stress inference network after training.

[0141] In some embodiments, the first loss can be determined using the following formula: ; ; ; In the formula, The first loss, The predicted loss for the pressure field of the first simulation sample. The discrete residual loss of the elliptic pressure Poisson equation for the first simulation sample. The number of samples in the first simulation. This is a predicted value for the pressure field. This represents the true value of the pressure field. Let be the number of grid cells in the pressure field. This represents the discrete residuals of the elliptic pressure Poisson equation.

[0142] In some embodiments, a second simulation sample including a concentration field at multiple time points can be obtained; an initial concentration inference network is constructed; the initial concentration inference network includes multiple cascaded Π-type convolutional modules, each Π-type convolutional module including a differential convolutional layer and a Π-type convolutional layer; the differential convolutional layer has a differential convolutional kernel with fixed weights, the fixed weights being determined by the discretization coefficient matrix of the spatial partial differential operator of the convection term; the Π-type convolutional layer includes multiple parallel trainable convolutional layers, element-wise multiplication layers, and combined convolutional layers, the multiple parallel trainable convolutional layers being used to process the same input physics layer respectively. Feature extraction is performed on the field to generate multiple feature maps. Element-wise multiplication layers are used to multiply the multiple feature maps element by element to obtain a multi-order polynomial combination result. Combination convolutional layers are used to fuse the multi-order polynomial combination result into a rheological term estimate. The mean square error of the difference between the time term and the convection and rheological terms in the concentration convection conservation equation corresponding to the second simulation sample is used as the second loss to train the initial concentration inference network. During the training process, the weight parameters of the differential convolutional layer are kept constant, and the trainable parameters of the Π-type convolutional layer are updated until the second loss meets the convergence condition, thus obtaining the concentration inference network.

[0143] Existing physical information neural network training typically requires a large amount of labeled data as supervision signals. However, in the application scenario of this application, the true distribution of the concentration field is difficult to obtain directly through experimental means. If labeled data is generated entirely through finite difference simulation, the generation cost is high and the amount of data is limited. Furthermore, if the soft-constrained PINN training method is used, and the physical residual is used as a loss function term and weighted and summed with the data loss, the physical constraint only indirectly affects the network output in the form of a penalty term. This cannot guarantee that the concentration field strictly satisfies the conservation equation in the forward calculation, and the difference in gradient magnitude between different loss terms in the convection-dominated problem can easily lead to training imbalance. Therefore, this embodiment uses a pure physical residual-driven unsupervised training method to train the initial concentration inference network.

[0144] A second set of simulated samples containing concentration fields at multiple time points is obtained. These samples provide initial conditions and a physical evolution process for training the initial concentration inference network. The concentration fields in the second set of simulated samples are numerical solutions obtained by solving the concentration convection conservation equation using the high-precision finite difference method. Each sample contains a sequence of concentration fields from the initial time point to multiple subsequent time points. During training, the second set of simulated samples does not provide labeled ground truth values ​​of the concentration fields as a supervisory target. Instead, it only provides the initial concentration field as the starting point for network inference, and the concentration field sequences at subsequent time points for calculating the discrete approximation of the time term in the concentration convection conservation equation. In other words, the concentration fields in the second set of simulated samples are not the labeled data that the network output should approximate, but rather intermediate reference values ​​in the physical residual calculation process. The concentration field output by the network is guided to an evolutionary trajectory that satisfies the concentration convection conservation equation under the constraint of the physical residual loss function, rather than being forced to fit a specific simulation result. Therefore, this training method can optimize and update network parameters without relying on labeled data.

[0145] An initial concentration inference network is constructed. This network receives the current concentration and velocity fields and processes them step-by-step through multiple cascaded Π-shaped convolutional modules, outputting the concentration field for the next time step. Each Π-shaped convolutional module consists of multiple cascaded Π-shaped convolutional modules. Each module includes two parallel branches: a differential convolutional layer and a Π-shaped convolutional layer. These two branches receive the same input, and their outputs are summed to form the module's output. When multiple Π-shaped convolutional modules are cascaded, the output of the preceding module serves as the input to the following module, progressively extracting higher-order physical features and nonlinear mappings.

[0146] The differential convolutional layer contains a differential convolution kernel with fixed weights, which are determined by the discretization coefficient matrix of the spatial partial differential operator for the convection term. The spatial partial differential operator includes... / x and / y represents the first-order spatial partial derivatives of the concentration field along the horizontal and vertical directions, respectively. The values ​​of each element in the discretization coefficient matrix are uniquely determined by the discretization form of the operator (such as forward difference, backward difference, or central difference), and remain constant after the differential convolutional layer is constructed, without changing during the training process. In the forward propagation of the initial concentration inference network, the differential convolutional layer performs a sliding window weighted summation operation on the input concentration field matrix through the fixed-weight differential convolution kernel to obtain the spatial partial derivative values ​​of the concentration field (i.e., the concentration gradient). Then, the concentration gradient is multiplied point-by-point with the input velocity field to output a numerical estimate of the convection term in the concentration convection conservation equation. The function of the differential convolutional layer is to solve the linear differential terms in the physical equations that can be precisely expressed through deterministic calculations, without introducing any trainable parameters.

[0147] Π-type convolutional layers are used to numerically estimate the rheological term in the concentration-convection conservation equation in the form of a high-order polynomial product. A Π-type convolutional layer includes multiple parallel trainable convolutional layers, element-wise multiplication layers, and combined convolutional layers. Multiple parallel trainable convolutional layers perform convolution operations on the same input physical field (a tensor formed by concatenating the concentration field and velocity field along the channel dimension) to extract spatial features at different scales, generating multiple feature maps. The number of channels in each feature map can be the same or different. The element-wise multiplication layer multiplies the values ​​of the multiple feature maps at the same spatial location point by point to obtain a product feature map. This operation is mathematically equivalent to constructing a high-order polynomial combination of the input features. When the number of parallel convolutional layers is N, the output of the element-wise multiplication layer corresponds to the Nth-order polynomial cross term of the input features—each feature map can be regarded as a certain nonlinear transformation of the input features, and their product is a high-order combination of the transformation results. The combined convolutional layer is used to perform convolution operations on the product feature map, fusing multi-channel information into a single-channel output, which is the numerical estimate of the rheological term. All trainable parameters of the Π-type convolutional layer (including the kernel weights and biases of the parallel convolutional layers and the kernel weights and biases of the combined convolutional layers) are updated during training through gradient backpropagation.

[0148] The mean squared error of the difference between the time term and the convection and rheological terms in the concentration convection conservation equation corresponding to the second simulation sample is used as the second loss to train the initial concentration inference network. During training, the weight parameters of the differential convolutional layer remain constant, and the trainable parameters of the Π-type convolutional layer are updated until the second loss meets the convergence condition, thus obtaining the concentration inference network. Specifically, in each training batch, the concentration field at a certain moment in the second simulation sample is input into the initial concentration inference network, and the network outputs the predicted value of the concentration field at the next moment. The values ​​of the three terms in the concentration convection conservation equation are calculated from the current concentration field, the current velocity field (calculated from the pressure field), and the predicted value of the concentration field at the next moment: the time term is obtained by subtracting the current concentration field from the predicted value of the concentration field at the next moment and dividing by a preset time step; the convection term is estimated by the convection term output by the differential convolutional layer; and the rheological term is estimated by the rheological term output by the Π-type convolutional layer. The difference between the time term and the sum of the convection and rheological terms is squared and averaged over all grid positions and all training samples to obtain the second loss.

[0149] During the gradient backpropagation of the second loss, the fixed weights of the differential convolutional layer remain constant and do not participate in parameter updates. Only the trainable parameters of the Π-type convolutional layer (including the weights and biases of parallel and combined convolutional layers) receive gradients and are updated. The physical significance of this freezing strategy is that the differential convolutional layer is responsible for handling the explicitly known linear differential terms in the physical equations, whose mathematical form is known and does not require training. Freezing ensures that this part of the calculation is not affected by the bias of the training data and always strictly follows the physical conservation laws. The Π-type convolutional layer is responsible for fitting the unknown nonlinear rheological terms in the physical equations that are difficult to express explicitly. Its parameters need to be learned through a data-driven approach to approximate this nonlinear mapping relationship. The different parameter update strategies of the two branches—the physically known part is not updated, and the physically unknown part is updated—reflect the different roles they play at the physical modeling level. During training, the numerical change of the second loss on the validation set is monitored in real time. When the second loss no longer continues to decrease or the decrease is lower than a preset threshold, it is determined that the second loss meets the convergence condition, training is stopped, and the current network parameters are fixed as the parameters of the concentration inference network after training.

[0150] In some embodiments, the second loss can be determined using the following formula: ; In the formula, This is the second loss. , and These represent the time term, convection term, and rheological term of the concentration-convection conservation equation. For fluid transport velocity, For vector differential operators, This refers to the volume concentration of the proppant. For time, This represents the number of samples in the second simulation.

[0151] In some embodiments, the two-phase flow of fluid within a hydraulic fracture can be simplified to a two-dimensional parallel flat plate laminar flow to construct a multiphysics coupled equation system including the concentration convection conservation equation, the pressure Poisson equation, the property constitutive equation, and the velocity equation; boundary conditions are set; the boundary conditions include: the left boundary of the hydraulic fracture is an injection boundary with constant proppant concentration, the right boundary of the hydraulic fracture is a tip boundary with constant pressure, the upper and lower boundaries of the hydraulic fracture are wall boundaries with no slip and no penetration, and the proppant concentration within the hydraulic fracture is zero at the initial moment; the multiphysics coupled equation system is spatially discretized on a discretized spatial grid, and solved according to the boundary conditions at a preset time step to obtain simulation samples including property fields, pressure fields, velocity fields, and corresponding concentration fields at multiple moments.

[0152] In existing technologies, numerical simulation of solid-liquid two-phase flow within three-dimensional hydraulic fractures is extremely costly. Directly solving the three-dimensional flow equations requires dividing the data into fine meshes in the length, height, and width directions, with the computational cost increasing cubically with the number of meshes. In reality, the length and height of a hydraulic fracture are much greater than its width, and the flow gradient along the fracture width is much smaller than that along the fracture length and height. Therefore, high-resolution solutions along the fracture width are unnecessary. To address this, this embodiment simplifies the solid-liquid two-phase flow within the hydraulic fracture as a two-dimensional parallel flat plate laminar flow, constructing a multiphysics coupled equation system including the concentration convection conservation equation, the pressure Poisson equation, the property constitutive equation, and the velocity equation.

[0153] The simplification of two-dimensional parallel flat plate laminar flow specifically means: neglecting the flow variation along the fracture width direction, reducing the three-dimensional flow problem to a two-dimensional flow problem along the fracture length direction (horizontal direction) and the fracture height direction (vertical direction). Under this simplification, the fracture width appears as the average width in the velocity equation, representing the average effect after integration along the fracture width direction. Solid-liquid two-phase flow is based on the Euler-Euler two-fluid continuous medium hypothesis—treating the proppant particle phase and the fracturing fluid-based liquid phase as two interpenetrating continuous media, each satisfying its own conservation equations. The two phases are coupled through concentration (concentration determines the density and viscosity of the mixture, and the mixture as a whole is treated as a quasi-single-phase fluid). The multiphysics coupling equation set includes the concentration convection conservation equation describing the spatiotemporal evolution of proppant concentration, the pressure Poisson equation describing the spatial distribution of the pressure field, the constitutive equation describing the changes in density and viscosity of the mixed fluid with concentration, and the horizontal and vertical velocity equations describing the relationship between the velocity field and the pressure field. The above set of equations couples the four physical fields—concentration field, pressure field, property field, and velocity field—together. The property field is determined by the concentration field through the constitutive equation, the pressure field is determined by the property field through the Poisson equation, the velocity field is determined by the pressure field through the velocity equation, and the concentration field is determined by the velocity field through the convection conservation equation, forming a complete closed-loop coupling relationship.

[0154] Define boundary conditions. Boundary conditions provide the boundary constraints and initial states required to solve the above partial differential equations, ensuring the uniqueness of the solution. Boundary conditions include boundary conditions and initial conditions. Boundary conditions include: setting the left boundary of the hydraulic fracture as an injection boundary with a constant proppant concentration, where the volumetric proppant concentration remains at a preset constant injection concentration value, representing the continuous supply of proppant carried into the fracture by the injected fluid during fracturing. Setting the right boundary of the hydraulic fracture as a tip boundary with constant pressure, where the fluid pressure remains at a preset formation pressure value, representing the pressure balance between the fracture tip and the reservoir matrix. Setting the upper and lower boundaries of the hydraulic fracture as wall boundaries with no slip and no penetration, where the normal velocity component of the fluid is zero (the wall is impenetrable), and the tangential velocity component has a large gradient near the wall due to viscosity but no slip, representing the rigid constraint of the fracture wall on the fluid flow. Initial conditions include: the proppant concentration within the hydraulic fracture is zero everywhere at the initial moment, representing the initial state before proppant injection into the fracture.

[0155] Based on the multiphysics coupling equations and boundary conditions, the multiphysics coupling equations are spatially discretized on a discretized spatial grid, and solved hourly according to a preset time step, yielding simulation samples including the physical property field, pressure field, velocity field, and corresponding concentration field at multiple time points. The discretized spatial grid divides the crack computational domain into several regular grid cells, and the values ​​of each physical quantity (concentration, pressure, density, viscosity, velocity) are defined at the grid nodes or grid centers. Specific methods of spatial discretization include, but are not limited to, using the finite difference method to approximate the partial derivative terms in the multiphysics coupling equations with the differences of the physical quantities at the grid nodes, transforming the partial differential equations into an algebraic equation system. Specific methods of time stepping include, but are not limited to, using explicit or implicit time discretization schemes, starting from the initial time, solving the algebraic equations step by step according to a preset time step, obtaining the numerical solutions of each physical quantity at each time step.

[0156] In the finite difference method implementation, the second-order spatial derivative term is discretized using a central difference scheme, and the first-order spatial derivative term is discretized using a forward or backward difference scheme. An algebraic equation is established for each grid node, and the distribution of physical quantities across the entire field at that moment is obtained by assembling and solving a large sparse linear system of equations. In the finite volume method implementation, the computational domain is divided into several control volumes. For each control volume, the partial differential equation is integrated and discretized into an algebraic equation. This method naturally satisfies the conservation of physical quantities and is more accurate in flux calculation. There are several implementation methods for time-progression. In the explicit time-progression method, the concentration value at the next moment is explicitly calculated from the concentration value at the current moment and the spatial term, without iterative solution, but the time step size is limited by stability conditions. In the implicit time-progression method, the concentration value at the next moment, together with the physical quantities at the current moment and the next moment, constitutes a coupled system of equations, which needs to be solved iteratively within each time step. The time step size is not limited by stability conditions, but the computational cost per step is high.

[0157] In the uniform mesh implementation, the grid step size is constant in both the horizontal and vertical directions, which is computationally simple but cannot automatically refine the mesh in high-gradient regions. In the adaptive mesh implementation, the mesh is automatically refined in concentration fronts or regions with abrupt viscosity changes to improve local resolution, while maintaining a coarse mesh in other regions to control computational load. The application of boundary conditions depends on the aforementioned boundary conditions and is achieved by setting fixed values ​​for boundary mesh nodes (concentration injection boundary) or fixed pressure values ​​(tip pressure boundary) and a constraint that the normal velocity at the wall is zero. After solving, the physical property field, pressure field, velocity field, and concentration field at each time step are stored in chronological order to form a simulation sample database. This simulation sample database is used to provide input-output paired samples for training the pressure inference network and the concentration inference network. For the pressure inference network, the physical property field in the simulation samples serves as the input, and the pressure field as the output target; for the concentration inference network, the concentration field sequence in the simulation samples is used to approximate the time term in the physical residual loss function.

[0158] In some embodiments, the mixed density equation and the equivalent viscosity equation can be constructed using the following formulas: ; ; In the formula, For mixed density, This refers to the volume concentration of the proppant. The true density of the proppant particles. The density of pure fracturing fluid. The equivalent viscosity of the mixture, The initial dynamic viscosity of the base fluid without proppant. This is the theoretical saturation volume concentration of the proppant. It is a negative rheological index.

[0159] The mixing density equation is used to characterize the relationship between the density of a fluid mixture and the proppant concentration. This equation is achieved through a volume fraction weighted average: when... When the density of the mixture is 0, the density of the mixture is equal to the density of the pure fracturing fluid. When the density approaches 1, the mixed density approaches the true density of the proppant.

[0160] The equivalent viscosity equation is used to characterize the nonlinear relationship between the equivalent viscosity of a fluid mixture and the proppant concentration. This equation is realized in power-law form: when... When the viscosity is 0, the mixing viscosity equals the base fluid viscosity. As it gradually increases from 0, the ratio ( ) increases, (1- As the effective viscosity decreases, the negative exponent β causes the reciprocal of this decrease to be amplified sharply, resulting in a rapid increase in the equivalent viscosity; when When the equivalent viscosity approaches the theoretical saturation volume concentration, it tends to infinity, which mathematically characterizes the physical phenomenon of a sharp increase in flow resistance during sand blockage.

[0161] Existing pressure inference networks typically treat the mixed density field and equivalent viscosity field as two independent input channels to the encoder when processing physical property fields, failing to differentiate the varying importance of different regions within the physical property field for pressure field solution. Regions with abrupt viscosity changes in the equivalent viscosity field (such as sand-blocking fronts) are often the main sources of pressure field error, yet these regions only occupy a small portion of the computational domain. If the encoder applies the same feature extraction weights to all regions, the network's computational resources cannot be concentrated on high-gradient regions, leading to insufficient accuracy in pressure field solution near viscosity interfaces. Furthermore, when abrupt changes in the concentration field cause a significant shift in the equivalent viscosity field compared to the previous moment, the pressure field needs to adjust accordingly within a short timeframe. However, in the encoder-decoder structure, physical property information requires multiple layers of downsampling and upsampling to reach the output. During deep feature transmission, physical property information may attenuate, causing the pressure field response to lag behind changes in the physical property field, potentially leading to iterative instability in long-term time-series inferences. Therefore, this embodiment incorporates a physical property field gradient-aware attention module, physical property residual connection paths, and adaptive gating units into the pressure inference network.

[0162] In some embodiments, the pressure inference network includes a property field gradient sensing attention module. The above-described inputting the current concentration field and property field into the pressure inference network to obtain the current pressure field may further include: inputting the current equivalent viscosity field into the module to obtain the spatial gradient magnitude distribution of the equivalent viscosity field, and generating an attention weight map based on this distribution, where the attention weights are positively correlated with the gradient magnitudes. The attention weight map is then multiplied element-wise with the feature maps output from each layer of the encoder to obtain a weighted feature map. The weighted feature map is then input into the decoder to reconstruct the pressure field.

[0163] The gradient-aware attention module for physical properties generates an attention weight map based on the drastic spatial changes in the equivalent viscosity field, automatically focusing the computational resources of the pressure inference network on the interface region with rapidly changing viscosity. The residual connection path directly transmits physical property information from the input to the decoder output, shortening the propagation path of physical property information from input to output and mitigating the attenuation of physical property information in deep networks. The adaptive gating unit dynamically adjusts the contribution weights of the residual connection path based on the global rate of change of the equivalent viscosity field, enhancing the response speed of the residual path when physical properties change abruptly and maintaining the stability of the main path when physical property changes are gradual, achieving a dynamic balance between transient response speed and steady-state computational accuracy in pressure field solving.

[0164] The equivalent viscosity field at the current moment is input into the property field gradient sensing attention module. This module calculates the spatial gradient magnitude distribution of the equivalent viscosity field and generates an attention weight map accordingly. The calculation of the spatial gradient magnitude distribution includes, but is not limited to: spatially differencing the equivalent viscosity field matrix along both the horizontal and vertical directions to obtain the horizontal and vertical gradient matrices; taking the square root of the sum of the squares of these two gradient matrices to obtain the gradient magnitude at each grid location. This gradient magnitude characterizes the degree of spatial variation of the equivalent viscosity at that location. The generation of the attention weight map includes, but is not limited to: using the gradient magnitude at each grid location as a weight value; locations with larger gradient magnitudes have higher weight values, and locations with smaller gradient magnitudes have lower weight values. The attention weight map has the same spatial dimension as the equivalent viscosity field.

[0165] In one implementation, the gradient magnitudes of all grid locations can be normalized to the [0,1] interval. The normalized values ​​are then directly used as attention weights, with the location having the largest gradient magnitude having a weight of 1 and the location having the smallest gradient magnitude having a weight of 0. In another implementation, a gradient magnitude threshold can be set. Grid locations with gradient magnitudes higher than this threshold have a weight of 1 (activation), while those lower have a weight of 0 (suppression). This method can completely eliminate interference from low-gradient regions, allowing the network's computational resources to be fully concentrated in high-gradient regions. However, discontinuities may exist on both sides of the boundary.

[0166] In some embodiments, the pressure inference network further includes a property residual connection path and an adaptive gating unit. The above-mentioned input of the weighted feature map into the decoder to reconstruct the pressure field may further include: during the decoding process, convolving the property field through the property residual connection path to obtain the property residual, and adding it point-by-point with the pressure field output by the decoder to obtain the candidate pressure field. The global rate of change of the equivalent viscosity field relative to the previous time step is calculated and input into the adaptive gating unit to obtain the gating weights, which are positively correlated with the global rate of change. Based on the gating weights, the candidate pressure field and the pressure field output by the decoder are combined to obtain the pressure field at the current time step.

[0167] The attention weight map is element-wise multiplied with the feature maps output by each layer of the encoder. During downsampling, the encoder outputs multiple feature maps, each with a fixed scaling factor relative to the input physical property field. Before element-wise multiplication, the attention weight map undergoes spatial dimension matching—it is scaled to the same spatial dimension as the current encoder layer's feature map using interpolation. Specific operations for element-wise multiplication include, but are not limited to: multiplying the weight value at each grid position in the attention weight map by the feature values ​​of all channels at the same spatial position in the encoder feature map (i.e., all channels at the same position share the same weight value), resulting in a weighted feature map. This weighting operation amplifies the responses in the encoder feature map corresponding to high gradient regions of equivalent viscosity and suppresses responses corresponding to low gradient regions, thus guiding the decoder to focus more on interface regions with drastic viscosity changes when reconstructing the pressure field.

[0168] The weighted feature map is input into the decoder, which recovers the spatial resolution through stepwise upsampling operations and skip connections, outputting the stress field. During the decoding process, the stress inference network also performs the following steps.

[0169] The current physical property field is input into the physical property residual connection path. This path consists of a set of combined convolutional layers used to extract features from the input physical property field (a two-channel tensor composed of a mixed density field and an equivalent viscosity field), outputting a tensor with the same spatial dimension as the decoder's output pressure field. This tensor serves as the physical property residual, used to correct the decoder's output pressure field point by point. The physical property residual is then added point by point to the decoder's output pressure field to obtain the candidate pressure field. This point-by-point addition operation includes, but is not limited to, element-wise addition of matrices of the same dimension—the value at each grid position in the candidate pressure field is equal to the sum of the corresponding value in the decoder's output pressure field and the corresponding value in the physical property residual. The role of the physical property residual connection path is to provide a pathway for physical property information that bypasses the encoder, allowing it to directly participate in the reconstruction of the pressure field based on deep features, mitigating the attenuation of physical property information after passing through deep networks.

[0170] The candidate pressure field retains the contribution of physical property information through the residual path, but the contribution weight of the residual path should not remain constant under all operating conditions. When the equivalent viscosity field changes significantly compared to the previous moment (such as a sharp increase in viscosity at the start of sand plugging), the contribution of the residual path needs to be increased to accelerate the transient response of the pressure field; when the equivalent viscosity field changes gradually, the contribution of the residual path should be decreased to maintain the stability of the pressure field and avoid introducing unnecessary disturbances. To this end, the global rate of change of the equivalent viscosity field is calculated and input into the adaptive gating unit.

[0171] The global rate of change characterizes the overall degree of change in the equivalent viscosity field at the current moment compared to the equivalent viscosity field at the previous moment. The calculation of the global rate of change includes, but is not limited to: traversing all grid positions of the equivalent viscosity field, calculating the relative difference between the equivalent viscosity value at the current moment and the equivalent viscosity value at the previous moment at each grid position, and averaging this relative difference across all grid positions to obtain the global rate of change. The mathematical form of the relative difference is the absolute value of the difference divided by the equivalent viscosity value at the previous moment (with a small constant added to prevent division by zero). A larger global rate of change indicates a more drastic change in the viscosity field.

[0172] The adaptive gating unit outputs the gating weights for the connection paths of the physical property residuals based on the global rate of change. The gating weights range from [0,1] and are used to control the contribution ratio of the physical property residuals in the pressure field reconstruction. The specific value of the gating weights is positively correlated with the global rate of change—the larger the global rate of change, the closer the gating weights are to 1, and the greater the contribution of the physical property residuals; the smaller the global rate of change, the closer the gating weights are to 0, and the smaller the contribution of the physical property residuals.

[0173] In one implementation, the global rate of change can be normalized and directly used as the gating weight, which is suitable for situations where the global rate of change is relatively uniformly distributed. In another implementation, the global rate of change can be mapped to a continuous value in the interval [0,1] using a logistic function. When the global rate of change is below a certain threshold, the gating weight approaches 0 (the residual path is almost completely closed), and when the global rate of change is above the threshold, the gating weight rapidly rises to approach 1 (the residual path is fully open).

[0174] The pressure field at the current moment is obtained by weighting the candidate pressure field and the pressure field output by the decoder based on the gating weights. The specific formula for the weighted combination includes, but is not limited to: P_final = g·P_candidate + (1-g)·P_decoder, where P_candidate is the candidate pressure field (including the contribution of physical property residuals), P_decoder is the original pressure field output by the decoder, and g is the gating weight output by the adaptive gating unit. When g=1, the pressure field is completely determined by the candidate pressure field (residual paths are fully open); when g=0, the pressure field is completely determined by the decoder output (residual paths are closed); when g is in the interval (0,1), the pressure field is the weighted average of the two. The pressure field after weighted combination is the final pressure field output at the current moment. This adaptive gating mechanism enables the pressure inference network to achieve fast transient response by relying on residual paths when viscosity changes abruptly, and to return to the main path to ensure steady-state accuracy when viscosity is flat, thus balancing the dual requirements of transient response speed and steady-state calculation accuracy.

[0175] Existing technologies typically assume a constant fracture width or use it as a simple preset parameter when dealing with fluid flow within hydraulic fractures, failing to consider the impact of fracture geometric deformation caused by pressure field changes during proppant transport on the flow and concentration fields. In actual hydraulic fracturing projects, fracture width is not fixed—it increases when the pressure within the fracture increases, causing the fracture wall to elastically open and the width to increase; conversely, it decreases when the pressure decreases due to compression under closure stress. Changes in fracture width directly alter the flow cross-sectional area, thus affecting the spatial distribution of the velocity field and the transport rate of the concentration field. Without considering this fluid-structure interaction effect, simulation accuracy will significantly decrease under conditions of high net pressure or geometrically non-uniform fractures. Therefore, this embodiment performs a fracture width update step after each iteration's property field update and before the pressure field solution. Meanwhile, due to the distortion of the physical mesh caused by the change in crack width, the standard fixed-weight differential convolution kernel can only be accurately equivalent to the finite difference operator on equidistant regular meshes. When used directly on non-uniform meshes, it will produce calculation errors. Therefore, the differential convolution kernel can be adaptively corrected to maintain its differential calculation accuracy.

[0176] In some embodiments, after each iteration of the physical property field update and before the pressure field solution, a crack width update is also performed: the crack width change is calculated based on the pressure field and the crack width spatial distribution field is determined; a coordinate transformation relationship is established based on the local rate of change of the crack width spatial distribution field; an equidistant physical grid is mapped to a non-uniform body-fitted computational grid according to the coordinate transformation relationship, with the grid density being positively correlated with the crack width change rate; and a fixed-weight differential convolution kernel is applied to the body-fitted computational grid to perform spatial partial differential operations.

[0177] The fracture width variation is calculated based on the pressure field to determine the spatial distribution field of the fracture width. The fracture width variation characterizes the degree of change in the width of the hydraulic fracture at each spatial location relative to the reference width. The calculation methods for the fracture width variation include, but are not limited to: based on the analytical or numerical relationship between fracture aperture and net pressure in elasticity, substituting the difference between the pressure value at each grid location in the pressure field at the current moment and the reference pressure value (such as the initial formation pressure) into this relationship to obtain the fracture width variation at each grid location. The reference width is the initial fracture width when not subjected to pressure disturbance. The reference width at each grid location at the current moment is added to the corresponding fracture width variation to obtain the actual fracture width value at each grid location at the current moment. The actual fracture width values ​​at all grid locations constitute the fracture width spatial distribution field. The spatial non-uniformity of this fracture width spatial distribution field reflects the influence of the spatial distribution differences of the pressure field on the fracture geometry—regions with high pressure have larger fracture widths (greater fracture opening), while regions with low pressure have smaller fracture widths (greater fracture closure).

[0178] In one implementation, the crack width variation can be independently calculated for each grid location using analytical formulas from linear elastic fracture mechanics and the relationship between one-dimensional crack aperture and net pressure (such as the crack width-pressure relationship in the PKN or KGD models). This implementation has the advantages of low computational cost and high speed, and is suitable for conditions with relatively regular crack morphologies. In another implementation, a simplified numerical model based on the finite element or boundary element method can be used to establish a displacement discontinuity equation on the crack surface and solve for the coupling relationship between the pressure field and the crack width field. This implementation can more accurately characterize the crack deformation response under non-uniform pressure distribution, but the computational cost is relatively high.

[0179] Coordinate transformation relationships are established based on the local rate of change of the slot width spatial distribution field. These relationships characterize the mapping from physical coordinate space to computational coordinate space. The local rate of change of the slot width spatial distribution field is the spatial gradient of the slot width value along the horizontal and vertical directions. Methods for establishing coordinate transformation relationships include, but are not limited to: associating the offset of the grid node position in the computational coordinate space relative to the physical coordinate space with the local rate of change of the slot width spatial distribution field—regions with a large rate of change of slot width have their distances in physical space magnified when mapped to computational space (mesh refinement); regions with a small rate of change of slot width have their distances in physical space reduced when mapped to computational space (mesh sparseness). The specific mathematical form of the coordinate transformation relationship is as follows: ξ / x = f( W / x) η / y = f( W / y), where (ξ,η) are the computational coordinates, (x,y) are the physical coordinates, W is the spatial distribution field of the seam width, and f is a monotonically increasing function. The specific value of this coordinate transformation relationship at each grid location is uniquely determined by the local rate of change of the seam width at that location.

[0180] In one implementation, physical coordinates can be directly mapped to computational coordinates using an analytical function based on algebraic transformation. The partial derivatives of the mapping function are determined by the local rate of change of the slot width spatial distribution field. In another implementation, the coordinate transformation relationship can be transformed into an elliptic partial differential equation based on differential equation transformation. Solving this equation yields a smooth mapping from physical space to computational space. This method ensures the continuity and smoothness of the mapping and is suitable for applications with drastic slot width variations and requiring fine control of mesh quality.

[0181] Based on coordinate transformation relationships, equidistant physical grids are mapped to non-uniform body-fitted computational grids, with grid density positively correlated with the crack width variation rate. An equidistant physical grid is a regular grid divided with a constant step size in physical coordinate space, where each grid cell has a consistent area and shape. A non-uniform body-fitted computational grid is an irregular grid formed in computational coordinate space by adjusting the grid node positions according to the crack width spatial distribution field. Specific mapping operations include, but are not limited to: for each grid node in physical coordinate space, calculating its corresponding position in computational coordinate space based on the coordinate transformation relationship at that node location. After mapping all grid nodes to computational coordinate space, the area and shape of each grid cell are no longer uniform—regions with a large crack width variation rate have densely distributed grid nodes, small grid areas, and high resolution; regions with a small crack width variation rate have sparsely distributed grid nodes, large grid areas, and low computational cost. The body-fitting meaning of the computational grid is that its grid boundary fits the geometric boundary of the hydraulic crack; grid nodes on the physical boundary always remain on the boundary and do not deviate from the boundary position due to coordinate transformation.

[0182] In some embodiments, the fixed-weight differential convolution kernel is an adaptive differential convolution operator with Jacobian matrix correction. When performing spatial partial differential operations on the physics field: the Jacobian matrix correction coefficients of the coordinate transformation relationship can be calculated. These correction coefficients are multiplied into the weight coefficients of the differential convolution kernel to obtain the adaptive convolution kernel weights at the current time step. A sliding-window weighted summation is then performed on the input physics field on the body-fitted computational grid using the differential convolution kernel with these adaptive weights, outputting the spatial partial derivative values ​​of the physics field.

[0183] A fixed-weight differential convolution kernel is applied to a body-fitted computational grid to perform spatial partial differential operations. When performing sliding-window weighted summation on a body-fitted computational grid, the differential accuracy of the standard fixed-weight convolution kernel is affected because the grid is non-uniform (it is equivalent to the weights of a finite difference operator on a uniform physical grid, but no longer precisely equivalent on a non-uniform grid). To ensure that the differential convolution kernel can still accurately perform spatial partial differential operations on a non-uniform body-fitted computational grid, this embodiment further constructs the fixed-weight differential convolution kernel as an adaptive differential convolution operator with Jacobian matrix correction.

[0184] When performing spatial partial differential operations on a physical field, the Jacobian matrix correction coefficients for the coordinate transformation relationship are calculated. The Jacobian matrix characterizes the degree of local scaling and deformation at each grid location during the coordinate transformation from physical coordinate space to computational coordinate space. The elements of the inverse Jacobian matrix of the two-dimensional coordinate transformation are determined by the partial derivatives of the slit width spatial distribution field. The Jacobian correction coefficients are the element values ​​in the inverse Jacobian matrix corresponding to the direction of that partial derivative.

[0185] The correction coefficients are multiplied into the weight coefficients of the differential convolution kernel to obtain the adaptive convolution kernel weights at the current time step. The original fixed weight coefficients of the differential convolution kernel correspond to the finite difference coefficients on the equidistant physical grid, which are suitable for calculating partial derivatives in physical coordinate space. On the body-fitted computational grid, the computational coordinate distance corresponding to the same physical distance is no longer constant, and the original weight coefficients need to be corrected using Jacobian correction coefficients. The specific correction operations include, but are not limited to: multiplying the weight coefficients in the differential convolution kernel corresponding to the partial derivatives in the horizontal direction by the element in the inverse of the Jacobian matrix (which represents the coordinate transformation from physical coordinate space to computational coordinate space) corresponding to the horizontal transformation ratio between physical and computational coordinates; and multiplying the weight coefficients in the differential convolution kernel corresponding to the partial derivatives in the vertical direction by the element in the inverse of the Jacobian matrix corresponding to the vertical transformation ratio between physical and computational coordinates. The corrected weight coefficients have different values ​​at different spatial locations on the body-fitted computational grid—the more severe the grid deformation, the greater the correction amount; and the gentler the grid deformation, the closer the correction amount is to the original fixed value. This modification essentially applies spatially dependent positional weighting to the weight coefficients of the convolution kernel, ensuring that when the same physical convolution template performs sliding window weighted summation at any position on a non-uniform grid, its output accurately corresponds to the partial derivative value in physical space, thereby maintaining the consistency of the differential computation capability of the fixed differential convolution kernel under different grid shapes.

[0186] In one implementation, the Jacobian correction coefficients can be directly multiplied into the fixed weight coefficients to obtain new convolution kernel weights, making the calculation simple and direct. In another implementation, a complete differential homeomorphism can be established between the physical coordinate space and the computational coordinate space. This transforms all partial differential operators in the physical space into differential operators in the computational space through coordinate transformation. This approach is not limited to simple multiplication and correction of the fixed weight coefficients but involves a complete coordinate transformation of the convolution kernel. This method offers higher accuracy in regions with extreme mesh deformation and drastic changes in the Jacobian matrix.

[0187] By applying a differential convolution kernel with adaptive weights to the input physics field on a body-fitted computational mesh using a sliding window weighted summation, the spatial partial derivative of the physics field is output. At each mesh node of the body-fitted computational mesh, the adaptively weighted differential convolution kernel uses the local adaptive weight coefficients at that node as a sliding window template to perform a weighted summation operation on the values ​​of the input physics field in the neighborhood of that node. The result of the weighted summation is the spatial partial derivative of the physics field in physical coordinate space at that mesh node location. This spatial partial derivative is then used for velocity field reconstruction calculations or to solve the convection term in the concentration convection conservation equation. This adaptive differential convolution mechanism enables the fixed-weight differential convolution kernel to accurately calculate spatial partial derivatives even on non-uniform meshes, achieving the embedding of the geometric nonlinear characteristics of fluid-structure interaction at the differential operator level.

[0188] In some embodiments, the Π-type convolutional layer is equipped with a spatially adaptive activation gating. The method may further include: inputting the current concentration field into the spatial activation gating unit to calculate the concentration value and concentration gradient magnitude at each spatial location, and comparing them with preset concentration thresholds and gradient thresholds, respectively. When the concentration value is lower than the concentration threshold and the gradient magnitude is lower than the gradient threshold, the location is determined to be a non-significant region, and a mask value of 0 is generated; otherwise, it is determined to be a significant region, and a mask value of 1 is generated. The first estimate output by the Π-type convolutional layer is multiplied element-wise with the mask to obtain a third estimate. The temporal rate of change of the concentration field at the current time can be calculated based on the third and second estimates.

[0189] When calculating rheological terms, the Π-type convolutional layer performs the same high-order polynomial product operation on all spatial locations within the computational domain, without distinguishing the significant differences in the degree of nonlinear rheological effects at different locations. In regions where the proppant concentration is close to zero (pure fracturing fluid zone), there are almost no nonlinear rheological effects, and in regions with gentle concentration gradients, the nonlinear effects are also very weak. However, the high-order polynomial product operation of the Π-type convolutional layer still performs the same nonlinear transformation in these regions as in high-concentration regions. This increases unnecessary computation and, because these regions lack effective physical residual monitoring signals, the Π-type convolutional layer may produce overfitting outputs in these regions, introducing unnecessary numerical perturbations in long-term time-series simulations. Therefore, the Π-type convolutional layer in this embodiment is equipped with a spatially adaptive activation gating.

[0190] Spatial adaptive activation gating is used to determine whether the high-order nonlinear fitting operation of the Π-type convolutional layer needs to be enabled at each position based on the spatial distribution characteristics of the concentration field. It generates a binary mask matrix with the same spatial dimension as the concentration field. The value of each position of the mask matrix is ​​used to control whether the output of the Π-type convolutional layer participates in the subsequent calculation at that position.

[0191] The concentration field at the current moment is input into the spatial activation gating unit. The spatial activation gating unit first calculates the concentration value and concentration gradient magnitude at each spatial location. The concentration value is directly taken from the element value of the concentration field matrix at that location. The concentration gradient magnitude is calculated as follows: the concentration field matrix is ​​spatially differencing along the horizontal and vertical directions to obtain the horizontal and vertical concentration gradient matrices. The square root of the sum of the squares of the two gradient matrices is used to obtain the concentration gradient magnitude at each grid location. The concentration gradient magnitude is used to characterize the degree of spatial variation of the concentration field at that location. This parameter, together with the concentration value, constitutes two independent dimensions for judging whether the nonlinear rheological effect is significant—regions with high concentration but gentle concentration gradients (such as the core area of ​​a fully deposited sand embankment) and regions with low concentration but steep concentration gradients (such as the location of a concentration front) may both have significant nonlinear rheological effects. Relying solely on the concentration value or solely on the concentration gradient is insufficient to fully cover the above two cases; both must be combined for judgment.

[0192] The concentration values ​​at each spatial location are compared with preset concentration thresholds, and the concentration gradient amplitude at each spatial location is compared with a preset gradient threshold. The preset concentration threshold is used to distinguish whether the proppant concentration is sufficient to produce a significant nonlinear rheological effect. This threshold is determined comprehensively based on the physical properties of the proppant particles (such as particle size, density, and theoretical maximum packing concentration) and the basic rheological properties of the fracturing fluid. The preset gradient threshold is used to distinguish whether the concentration gradient is sufficient to allow spatial activation gating to determine that the location is a significant region of nonlinear rheological effect. When the concentration value at a spatial location is lower than the preset concentration threshold and the concentration gradient amplitude is lower than the preset gradient threshold, the location is determined to be a non-significant region (i.e., the contribution of the nonlinear rheological effect at this location is negligible), and a mask value of 0 is generated. Otherwise, the location is determined to be a significant region (i.e., the contribution of the nonlinear rheological effect at this location is not negligible), and a mask value of 1 is generated. The judgment logic is executed sequentially for all spatial locations, eventually forming a binary matrix (i.e., a spatial activation mask) with the same spatial dimension as the concentration field. In this matrix, the set of grid locations with a value of 1 represents the salient region of the nonlinear rheological effect, while the set of grid locations with a value of 0 represents the non-salient region.

[0193] The rheological term estimate output by the Π-type convolutional layer is multiplied element-wise with the spatial activation mask to obtain the gated rheological term estimate. The specific operation of this element-wise multiplication is as follows: traverse all grid positions of the rheological term estimate matrix, and multiply the rheological term estimate at each position by the mask value (0 or 1) of the spatial activation mask at the same position. When the mask value is 0, the product result at that position is 0, meaning the rheological term estimate of the Π-type convolutional layer at that position is completely masked, and the rheological term contribution at that position is zero. When the mask value is 1, the product result at that position is equal to the original value of the rheological term estimate of the Π-type convolutional layer at that position, meaning the rheological term estimate of the Π-type convolutional layer at that position is fully preserved. Therefore, the gated rheological term estimate obtained in this way is non-zero only at positions where nonlinear rheological effects are significant, and zero at other positions. Based on the gated rheological and convection term estimates, the temporal rate of change of the concentration field at the current moment is calculated. At this time, the concentration field change in insignificant regions is driven only by the contribution of the convection term, while the concentration field change in significant regions is driven by both the convection and rheological terms.

[0194] In one implementation, a binary mask with values ​​of 0 or 1 can be generated as described above. In another implementation, the spatial activation gating unit outputs a soft mask with continuous values ​​in the interval [0,1]. The mask value is obtained by mapping the concentration value and the concentration gradient magnitude through a sigmoid function—when both the concentration value and the concentration gradient magnitude are very low, the mask value is close to 0; as either value increases, the mask value gradually and smoothly rises to near 1. This method creates a smooth transition between salient and insignificant regions, avoiding the numerical discontinuities that may be introduced by hard gating.

[0195] In some embodiments, the Π-type convolutional layer is a multi-timescale Π-type convolutional layer, including a fast-channel sublayer and a slow-channel sublayer. The above-described method of inputting the current concentration field and velocity field into the differential convolutional layer in the concentration inference network to obtain a first estimate of the convection term in the concentration convection conservation equation may further include: inputting the current concentration field and velocity field into the fast-channel sublayer, and outputting the transient response component in the rheological term through a multi-order polynomial product operation, used to characterize the immediate contribution of the nonlinear coupling between physical quantities at the current moment to the evolution of the concentration field. Obtaining a concentration field sequence from multiple consecutive historical moments, including the current moment, inputting it into the slow-channel sublayer, and outputting the gradually changing trend component in the rheological term through temporal feature extraction and a multi-order polynomial product operation, used to characterize the long-term cumulative contribution of inter-particle interactions during the continuous deposition of the proppant. Superimposing the transient response component and the gradually changing trend component yields the first estimate.

[0196] In actual physical processes, rheological effects comprise two components with distinct properties: the first is the transient response component, which includes the immediate feedback of viscosity to concentration changes and the real-time effect of velocity gradients on the concentration field. This component is directly related to the current physical state, changing rapidly and over a short timescale. The second is the gradual trend component, which represents the long-term evolutionary trend of gradually increasing interparticle contact and effective stress during proppant accumulation. This component depends on the accumulated history of the concentration field, changing slowly over a long timescale. A single Π-type convolutional layer can only learn the input-output mapping at a fixed time step and cannot distinguish between the two different timescales of physical processes, resulting in insufficient ability to capture the gradual trend component. In long-term time series extrapolation, gradual errors may gradually accumulate. Therefore, the Π-type convolutional layer in this embodiment is a multi-timescale Π-type convolutional layer, including fast-channel Π-type convolutional sublayers and slow-channel Π-type convolutional sublayers.

[0197] The fast-channel Π-type convolutional sublayer is used to output the transient response component of the rheological term, taking the current concentration and velocity fields as inputs. The fast-channel Π-type convolutional sublayer has the same structure as the basic Π-type convolutional layer, consisting of a combination of multiple parallel trainable convolutional layers, element-wise multiplication layers, and combined convolutional layers. The input to the fast-channel sublayer is a tensor formed by concatenating the current concentration and velocity fields along the channel dimension. After multiplication by multiple polynomials, the transient response component is output. This transient response component characterizes the instantaneous contribution of the nonlinear coupling between physical quantities at the current moment to the evolution of the concentration field—that is, the rheological effect directly determined by the current concentration and velocity states, such as the adjustment of flow resistance caused by the instantaneous change in viscosity with concentration, and the instantaneous feedback of the velocity gradient to the concentration distribution.

[0198] The slow-channel Π-type convolutional sublayer is used to take a concentration field sequence from multiple consecutive historical moments, including the current moment, as input and output the slowly varying trend component in the rheological term. The slow-channel Π-type convolutional sublayer adds temporal feature extraction capabilities to the basic Π-type convolutional layer. Its structure includes a temporal feature extraction module and a spatial polynomial product module. The temporal feature extraction module extracts a temporal feature map reflecting the long-term evolution trend of the concentration field from the input multi-moment concentration field sequence. The spatial dimension of this temporal feature map is the same as the concentration field at the current moment, and the temporal feature values ​​at each grid position characterize the cumulative change trend of the concentration in the time dimension at that position (such as the rate of continuous increase in concentration, the acceleration of concentration change, etc.). The spatial polynomial product module performs multi-order polynomial product operations on the temporal feature map and the velocity field at the current moment, outputting the slowly varying trend component. This slowly varying trend component characterizes the long-term cumulative contribution of interparticle interactions during the continuous proppant deposition process—such as the increase in effective stress due to the gradual increase in sand embankment height as the proppant injection amount increases, and the gradual strengthening of the interparticle contact force network, etc., which evolve slowly over time.

[0199] In one implementation, a one-dimensional temporal convolutional layer can be used to perform convolution operations on the input multi-time-step concentration field sequence along the time dimension, extracting the local trend features of concentration changes over time at each spatial location. In another implementation, a simplified version of a gated recurrent unit or long short-term memory network can be used to successively update the hidden state along the time dimension, extracting the long-term evolution features of concentration at each spatial location over time. This approach can model temporal dependencies over a longer range and is suitable for characterizing the precursor features of sand blockage that appear after proppant has been continuously accumulating for hundreds of time steps.

[0200] The transient response component is superimposed with the slowly varying trend component to obtain the final output of the multi-timescale Π-shaped convolutional layer, which is the complete estimate of the rheological term. This superposition operation is an element-wise addition of two matrices of the same dimension, and each element of the output matrix is ​​equal to the sum of the corresponding elements of the transient response component matrix and the slowly varying trend component matrix. This complete rheological term estimate is then multiplied element-wise by a spatially adaptive activation-gated mask and superimposed with the convection term estimate to calculate the temporal rate of change of the concentration field.

[0201] The parameter update strategies for the fast and slow channel sublayers are the same—all trainable parameters of both sublayers are updated through backpropagation of the physical residual loss during the training of the concentration inference network. However, the sources of the gradient signals received by the parameters during training differ for the two channels: the parameters of the fast channel sublayer are mainly affected by the physical residual at the current time step, while the parameters of the slow channel sublayer are constrained by temporal information provided by concentration field sequences from multiple historical time steps, in addition to being affected by the physical residual at the current time step. Therefore, the supervision signals received by the two sublayers during training differ in time scale. This multi-timescale structure enables the Π-type convolutional layer to simultaneously handle nonlinear rheological effects at two different time scales—transient response and long-term evolution—within the same framework.

[0202] The following is a specific embodiment of this specification: Figure 2 The complete implementation process of this specific embodiment is shown. Figure 3 This diagram illustrates the internal architecture of the Π-type convolution module in this specific embodiment. Figure 4 This diagram illustrates the stress inference network structure and boundary mask layer in this specific embodiment. Figure 5 The diagram shows the time-flow loop architecture of the coupled system in this specific embodiment. Figure 6 The forward difference is shown in this specific embodiment. (·) / x of 3 3. Schematic diagram of the differential convolution operator; Figure 7 The forward difference is shown in this specific embodiment. (·) / y's 3 3. Schematic diagram of the differential convolution operator; Figures 8-17 Table 1 shows a comparison of the simulation cloud maps and errors of existing technologies ConvLSTM, PINN, and the above-mentioned proppant transport simulation method of the present invention with the reference concentration field of FDM. Table 2 shows a comparison of the calculation time of different simulation methods and a comparison of the quantitative error of the concentration field at typical times. It can be seen that the above-mentioned proppant transport simulation method of the present invention achieves a balance between computational efficiency and physical consistency, and ensures the accuracy and stability of long-term extrapolation.

[0203] Table 1

[0204] Table 2

[0205] Based on the above-described proppant transport simulation method, this specification also provides embodiments of a proppant transport simulation device. For example... Figure 18 As shown, the proppant transport simulation device 1800 may specifically include the following modules: The pressure inference module 1801 is used to input the concentration field of the proppant and the physical property field of the fluid in the hydraulic fracture at the current moment into the pressure inference network to obtain the pressure field of the fluid at the current moment. The pressure inference network includes a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy the preset isobaric boundary conditions and / or zero flux boundary conditions. Calculation module 1802 is used to calculate the velocity field of the proppant at the current moment based on the pressure field; The concentration inference module 1803 is used to input the current concentration field and velocity field into the concentration inference network to obtain the concentration field and physical property field at the next moment. The concentration inference network is coupled with multiple Π-type convolutional modules. The Π-type convolutional modules include differential convolutional layers and Π-type convolutional layers. The differential convolutional layers are used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolutional layers are used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution. The rheological term is used to characterize the rheological feedback effect of proppant viscosity change and particle interaction on the concentration distribution. The iteration module 1804 is used to take the next moment as the new current moment and iteratively execute the above steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field until the preset simulation time is reached, so as to obtain the simulation results of the proppant concentration field at each moment.

[0206] This specification also provides a computer device for a proppant transport simulation method based on a physical coding neural network, including a processor and a memory for storing processor-executable instructions. In specific implementations, the processor can execute the above-described proppant transport simulation method according to the instructions.

[0207] To execute the above instructions more accurately, please refer to... Figure 19 As shown in the embodiments of this specification, another specific computer device 1900 is also provided, wherein the computer device 1900 includes a network communication port 1901, a processor 1902 and a memory 1903, and the above structures are connected by internal cables so that the various structures can perform specific data interaction.

[0208] The processor 1902 can be specifically used to execute the above-mentioned proppant transport simulation method.

[0209] The memory 1903 can be used to store the corresponding instruction program.

[0210] In this embodiment, the network communication port 1901 can be a virtual port bound to different communication protocols, thereby enabling the sending or receiving of different data. For example, the network communication port can be a port responsible for web data communication, a port responsible for FTP data communication, or a port responsible for email data communication. Furthermore, the network communication port can also be a physical communication interface or communication chip. For example, it can be a wireless mobile network communication chip, such as GSM or CDMA; it can also be a Wi-Fi chip; or it can be a Bluetooth chip.

[0211] In this embodiment, the processor 1902 can be implemented in any suitable manner. For example, the processor can take the form of a microprocessor or processor and a computer-readable medium storing computer-readable program code (e.g., software or firmware) executable by the (micro)processor, logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers, and embedded microcontrollers, etc. This specification is not limiting.

[0212] In this embodiment, the memory 1903 includes volatile memory and non-volatile memory. The memory 1903 can include multiple layers. In digital systems, anything that can store binary data can be a memory; in integrated circuits, a circuit with storage function but no physical form is also called a memory, such as RAM, FIFO, etc.; in a system, a storage device with a physical form is also called a memory, such as a memory stick, TF card, etc.

[0213] Furthermore, embodiments of this specification also provide a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described... Figure 3 Instructions for the proppant transport simulation method shown.

[0214] Furthermore, embodiments of this specification provide a computer program product, including a computer program that, when executed by a processor, implements the above-described... Figure 1 The method for simulating proppant transport is shown.

[0215] It should be understood that in the various embodiments of this specification, the sequence number of each process does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this specification.

[0216] It should also be understood that, in the embodiments of this specification, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Additionally, the character " / " in this specification generally indicates that the preceding and following related objects have an "or" relationship.

[0217] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) that include computer-usable program code.

[0218] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0219] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0220] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational tasks to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The task is a function specified in one or more boxes.

[0221] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for simulating proppant transport based on a physical coding neural network, characterized in that, include: Input the concentration field of the proppant and the physical property field of the fluid in the hydraulic fracture at the current moment into the pressure inference network to obtain the pressure field of the fluid at the current moment. The pressure inference network includes a Poisson regularization term and a boundary mask term. The Poisson regularization term is used to constrain the deviation of the fluid pressure field from the elliptic pressure Poisson equation, and the boundary mask term is used to constrain the fluid pressure field to satisfy preset isobaric boundary conditions and / or zero flux boundary conditions. Calculate the velocity field of the proppant at the current moment based on the pressure field; The current concentration field and velocity field are input into the concentration inference network to obtain the concentration field and physical property field at the next moment. The concentration inference network is coupled with multiple Π-type convolutional modules. The Π-type convolutional modules include differential convolutional layers and Π-type convolutional layers. The differential convolutional layers are used to perform spatial partial differential operations based on fixed differential weights to fit the convection term estimate in the concentration convection conservation equation. The Π-type convolutional layers are used to perform multi-order polynomial product operations to fit the rheological term estimate in the concentration convection conservation equation. The convection term is used to characterize the convective transport effect of fluid flow on the concentration distribution. The rheological term is used to characterize the rheological feedback effect of proppant viscosity changes and particle interactions on the concentration distribution. The next moment is taken as the new current moment, and the above steps of solving the pressure field, calculating the velocity field, and updating the concentration field and physical property field are iteratively executed until the preset simulation time is reached, so as to obtain the simulation results of the proppant concentration field at each moment.

2. The method according to claim 1, characterized in that, The method further includes: Construct an elliptic pressure Poisson equation using the following formula: ； In the formula, For vector differential operators, Fluid pressure, For fluid mixing density, The equivalent viscosity of the fluid. This is the acceleration due to gravity.

3. The method according to claim 1, characterized in that, The method further includes: Concentration convection conservation equations are constructed using the following formula: ； ； ； ； In the formula, , and These are the time term, convection term, and rheological term, respectively. Fluid transport velocity, including horizontal transport velocity. and vertical transport speed , The sliding settlement velocity, For vector differential operators, Fluid pressure, For fluid mixing density, The equivalent viscosity of the fluid. It is the acceleration due to gravity. The average width of the hydraulic fracture. This refers to the volume concentration of the proppant. The true density of the proppant particles. The density of pure fracturing fluid. The equivalent spherical particle size of the proppant, For horizontal step size, The vertical step size, For time.

4. The method according to claim 1, characterized in that, The method further includes: Obtain the first simulation sample, which includes the physical property field and the corresponding pressure field at multiple time points; Construct an initial stress inference network; the initial stress inference network includes an encoder, a decoder and a skip connection layer, the skip connection layer is used to concatenate the shallow spatial features output by each layer of the encoder with the deep semantic features output by the corresponding layer of the decoder; A boundary mask layer is set at the output of the decoder. The boundary mask layer is used to constrain the pressure field output by the decoder according to the preset boundary mask matrix. The initial pressure inference network is trained by using the sum of the pressure field prediction loss of the first simulation sample and the discrete residual loss of the elliptic pressure Poisson equation as the first loss until the first loss satisfies the convergence condition, thus obtaining the pressure inference network.

5. The method according to claim 1, characterized in that, The method further includes: Obtain a second simulation sample that includes concentration fields at multiple time points; An initial concentration inference network is constructed. The initial concentration inference network includes multiple cascaded Π-type convolutional modules. Each Π-type convolutional module includes a differential convolutional layer and a Π-type convolutional layer. The differential convolutional layer has a differential convolutional kernel with fixed weights, which are determined by the discretization coefficient matrix of the spatial partial differential operator of the convection term. The Π-type convolutional layer includes multiple parallel trainable convolutional layers, element-wise multiplication layers, and combined convolutional layers. The multiple parallel trainable convolutional layers are used to extract features from the same input physical field to generate multiple feature maps. The element-wise multiplication layers are used to multiply the multiple feature maps element-wise to obtain a multi-order polynomial combination result. The combined convolutional layer is used to fuse the multi-order polynomial combination result into a rheological term estimate. The mean square error of the difference between the time term and the convection and rheological terms in the concentration convection conservation equation corresponding to the second simulation sample is used as the second loss to train the initial concentration inference network. During the training process, the weight parameters of the differential convolutional layer are kept constant, and the trainable parameters of the Π-type convolutional layer are updated until the second loss meets the convergence condition, thus obtaining the concentration inference network.

6. The method according to claim 4 or 5, characterized in that, The method further includes: The two-phase flow of fluid within a hydraulic fracture is simplified to a two-dimensional parallel flat plate laminar flow, in order to construct a multi-physics field coupled equation set including the concentration convection conservation equation, the pressure Poisson equation, the property constitutive equation, and the velocity equation. Define boundary conditions; the boundary conditions include: the left boundary of the hydraulic fracture is an injection boundary with constant proppant concentration, the right boundary of the hydraulic fracture is a tip boundary with constant pressure, the upper and lower boundaries of the hydraulic fracture are wall boundaries with no slippage and no penetration, and the proppant concentration in the hydraulic fracture is zero at the initial moment. The multiphysics coupling equations are spatially discretized on a discretized spatial grid and solved according to the predetermined time step based on the boundary conditions, resulting in simulation samples including the physical property field, pressure field, velocity field and corresponding concentration field at multiple time points.

7. The method according to claim 1, characterized in that, The physical property field includes a mixed density field and an equivalent viscosity field; the calculation of the proppant velocity field at the current moment based on the pressure field includes: Calculate the horizontal and vertical pressure gradients based on the pressure field; The horizontal velocity field is calculated based on the equivalent viscosity field and horizontal pressure gradient at the current moment; the horizontal velocity field is used to characterize the proppant migration velocity along the length of the hydraulic fracture. Calculate the vertical velocity field of the base based on the current mixed density field and vertical pressure gradient; Calculate the sedimentation correction velocity field based on the current concentration field; The vertical velocity field is determined based on the basic vertical velocity field and the settlement-corrected velocity field; the vertical velocity field is used to characterize the proppant migration velocity in the crack height direction. Determine the velocity field at the current moment based on the horizontal and vertical velocity fields.

8. The method according to claim 1, characterized in that, The step of inputting the current concentration field and velocity field into the concentration inference network to obtain the concentration field and property field at the next moment includes: The concentration field and velocity field at the current moment are input into the differential convolutional layer in the concentration inference network to obtain the first estimate of the convection term in the concentration convection conservation equation; The concentration field and velocity field at the current moment are input into the Π-type convolutional layer in the concentration inference network to obtain the second estimate of the rheological term in the concentration convection conservation equation; Based on the first and second estimates, calculate the time rate of change of the concentration field at the current moment; The concentration field increment from the current moment to the next moment is calculated based on the rate of change over time and a preset time step, so as to determine the concentration field at the next moment. Based on a preset concentration threshold, the concentration field at the next moment is saturated and clamped. The concentration field at the next moment after clamping is input into the fluid's constitutive equation to obtain the fluid's property field at the next moment; the constitutive equation is used to characterize the mapping relationship between the fluid's mixing density and equivalent viscosity and the proppant concentration.

9. The method according to claim 8, characterized in that, The constitutive equations of the physical properties include a mixed density equation and an equivalent viscosity equation; the method further includes: The mixed density equation and equivalent viscosity equation are constructed using the following formulas: ； ； In the formula, For mixed density, This refers to the volume concentration of the proppant. The true density of the proppant particles. The density of pure fracturing fluid. The equivalent viscosity of the mixture, The initial dynamic viscosity of the base fluid without proppant. This is the theoretical saturation volume concentration of the proppant. is the rheological index.

10. The method according to claim 1, characterized in that, The method further includes: The concentration field, pressure field and velocity field at each time point are inversely normalized to obtain the concentration field cloud map, pressure field cloud map and velocity field cloud map at each time point, which are used to characterize the sandbank accumulation morphology and concentration front migration position, pressure gradient distribution characteristics and fluid flow direction of the proppant in the hydraulic fracture space. The concentration field at each time point is compared with the reference concentration field to determine the simulation accuracy. The simulation time is compared with that of the benchmark numerical simulation method to determine the simulation efficiency.

Images