JP2023099406A - Technique for digitalized multiphase fluid-structural coupled (fsi) penetration numerical value simulation in interior rock core experiment - Google Patents

Abstract

To provide a technique for a multiphase fluid-structural coupled (FSI) penetration numerical value simulation in an interior rock core experiment.SOLUTION: A technique includes four steps of: measuring a parameter of a sample after a rock core sample is dried, performing vacuum air exhaustion, and applying pressure on saturated water or oil (step 1); performing a nuclear magnet resonance scan and acquiring a two-dimensional image (step 2); acquiring a three-dimensional data volume of the rock core, constructing a three-dimensional irregular void network model, substituting the three-dimensional data volume into the node of the model, calculating a hole-slot diameter adjacent to adjacent nodes, and performing simulation calculation of a penetration rate of the model (step 3); securing the penetration rate of the rock core close to an actual value, constructing a void network model of the rock core, and performing a penetration numerical value simulation by combining the model and an unsteady fluid-structural coupled (FSI) multiphase penetration numerical value simulation technique (Step 4).SELECTED DRAWING: Figure 1

Description

Field of the Invention The present invention relates to the field of oil and gas field development, especially to the digitalized multi-phase fluid-structure interaction (FSI) infiltration numerical simulation method in indoor rock core experiments, which is applied to general sandstone oil and gas minerals, In addition, it is also applied to non-general oil and gas minerals such as dense oil and gas, shale oil and gas and natural gas hydrates, and is also applied to CO2 geological sealed storage technology.

Oil and natural gas are one of the important energies for maintaining the rapid development of the national economy. Whether to do so has long been an important issue during the oil and gas field development process. Since the structure of pits and throats inside rocks in underground deposits is actually complicated, it is difficult to clarify the rules of fluid infiltration in them by experimental means. Many researchers look for ways to help improve oil extraction rates by using porous media models to simulate the flow of different types of fluids inside rocks. Digital core (digital rock core) technology, as a branch of the porous media model, can be used in various fields such as geology, seismic, geological logging, exploitation and extraction enhancement in the oil and natural gas industry. The digital core technology can preserve the microscopic physical characteristics of the rock core, so it can ensure that the rock can be used infinitely, and is an important rock physical experiment numerical simulation platform. The effects of various microscopic factors (e.g., pore connectivity and wettability) on the reservoir can be quantitatively studied, in addition to physical properties that cannot be directly measured by traditional physical experiments (e.g., The relative permeability of the three phases of oil, gas, and water) can be calculated. In view of its wide range of applications, the research on digital cores has important implications for improving the extraction rate of oil and natural gas.

In the application of digital cores to conduct various rock physics experiments, the most basic task is to construct a three-dimensional digital core model that is accurate and suitable for actual rocks. In the past 15 years, with the innovation of experimental instruments and the breakthrough of new theories, domestic and foreign research teams have continuously put forward new construction methods of digital core models. After years of research, digital rock core construction methods are divided into three types: numerical reconstruction method, physical experiment method and mixing method. Numerical reconstruction methods are based on a small amount of 2D thin-section images, exploit the information contained in the 2D images, and additionally reconstruct a 3D digital core through random simulation methods or sedimentary rock process simulation methods. It is a way to The accuracy and modeling efficiency of the model constructed by this method are relatively low, and the selection of constraints in the modeling method causes the simulation results to be random, so the true rock core characteristics of the reserve stratum cannot be obtained. Hard to reproduce. The physical experiment method uses laboratory instruments to photograph or scan rock core samples to obtain a large number of two-dimensional photographs of rock cores, and then through modeling programs or software to obtain two-dimensional images. Refers to superimposing or reconstructing a dimensional photograph onto a three-dimensional digital core. However, the method is limited by the resolution and accuracy of the laboratory equipment (e.g., CT scanner), so the scale of the model to be constructed is relatively small (generally millimeter scale), and the representativeness of the model and its industrial application are limited. is greatly restricted. Extraction and analysis of microscopic parameters of rock cores with stalactite fracture features are difficult, costly physical experiments, and long cycles. The mixed method can build a relatively accurate three-dimensional model by combining various modeling methods and taking advantage of each, but the scale of the model to be built is still different from the actual rock core. At the same time, limited by the incompleteness of the pore-scale infiltration theory and the mathematical model, the simulation results in the model constructed by the method differ from the actual experimental results to varying degrees. In recent years, with the rapid progress of magnetic resonance imaging (MRI) and the remarkable increase in computing power (hash rate) of computer GPU chips, digitized three-dimensional rock samples corresponding to real rock samples in laboratory experiments have been developed. Build a pore network model of the rock core, develop a pore network model multiphase infiltration theory, a mathematical model and numerical simulation method, and digitize the pore network model of the rock core and the pore network model multiphase infiltration numerical value Combining the simulation methods, an important physical foundation was laid for carrying out the rock core seepage simulation study. At the same time, it has created the technical method of digitized rock core pore network model and pore network model multi-phase infiltration numerical simulation method, which has great significance for the development of oil and natural gas related industries. .

The purpose of the present invention is to overcome the shortcomings of existing technology, specifically to provide a digitized multi-phase fluid-structure interaction (FSI) infiltration numerical simulation method for laboratory rock core experiments. Then, by combining the digitized pore network model of the rock core with the fluid-structure interaction (FSI) seepage numerical simulation method, the process of the laboratory rock core flow experiment can be simulated and reproduced. Digitalized single-phase and multi-phase (oil-water, gas-water and oil-gas) fluid-structure interaction (FSI) seepage numerical simulation analysis and testing in core experiments can be realized.

The following technical solutions achieve the objectives of the present invention.

A digitized multiphase fluid-structure interaction (FSI) infiltration numerical simulation method for laboratory rock core experiments, which includes the following steps.

Step 1: First select the rock core sample, after drying, measure the hole/throat length, porosity and permeability of the rock core, additionally evacuate the rock core sample, simulating formation water or formation crude oil for saturation;
Step two: performing nuclear magnetic resonance MRI/T ₂ scan measurements on the rock core sample to obtain two-dimensional images of the rock core at different cross-sections;
Step 3: Performing interpolation on the 2D images to obtain a 3D data volume of the rock core MRI/ _T2 ,
Step 4: constructing a three-dimensional random void network model and substituting the MRI/T ₂ three-dimensional data volume into the nodes of the three-dimensional random void network model;
Step 5: Calculate the hole throat radius of each neighboring node in the three-dimensional irregular pore network model through the conversion factor α, and simulate the permeability of the digitized pore network model of the rock core; In addition, by adjusting the magnitude of the conversion factor α, we ensured that the permeability of the pore network was close to the measured permeability of the rock core, and in addition, the pore network of the digitized rock core building a model;
Step 6: Combining the digitized rock core pore network model with the unsteady fluid-structure interaction (FSI) multiphase infiltration numerical simulation method to perform the rock core fluid-structure interaction (FSI) infiltration numerical simulation. Contain six steps.

The second step specifically includes the following contents. 1. Using a nuclear magnetic resonance apparatus, perform a nuclear magnetic resonance scan on the rock core sample after the evacuation and saturation numerical simulation, and additionally generate a _T2 distribution using the SIRT inversion algorithm. 2. Section processing is performed for each position of the rock core. After selecting a number of nuclear magnetic resonance scans for the rock core to acquire rock core nuclear magnetic resonance MRI two-dimensional images at different section positions of the rock core, the section position coordinates and the pixel data of the two-dimensional images are divided into two Containing two contents of saving in dimensional image TXT text.

The third step specifically includes the following contents. Using an interpolation algorithm and additionally combining the 2D image TXT text acquired in step 2 and the rock core hole/throat length in step 1, the data volume in nuclear magnetic resonance imaging, i.e. the nucleus of the rock core Interpolation is performed on the magnetic resonance MRI 2D images so that the data volume size is adapted to the requirements of the microscopic hole and throat features of the reaction rock, which is then used for rock core MRI/T. Including obtaining _two 3D tensor data volumes.

The step 4 specifically includes the following substeps.

S401: Set the size of the model, the number of ligands and the average hole/throat length according to the scale of the rock core,
S402: Adopting the three-dimensional regular cubic network structure of the computer programming language and matrix computing library Eigen structure to generate an X×Y×Z three-dimensional regular cubic network;
S403: The total number of nodes in the cubic network is (X - 1 x (Y - 1) x (Z - 1), each node represents one gap, and the nodes are connected through the throat , installed so that the other part is a rock framework,
S404: 1, in the constructed cubic network model, six throats are connected to each other around each node representing the void, the throat length is the rock average hole-throat length <l>, x of the cubic network model, Let the side lengths in the y and z directions be L _x =(X-1)<l>, _Ly =(Y-1)<l> and L _z =(Z-1)<l>, respectively, and the cube Between all grid nodes in the network model, realization of perfect connections through round tubes, and the ratio of the radii of the _voids to the _throats of 1; and removing all points in the yoz plane of each layer whose distance from the center point is greater than _Ly , so that the cube network model becomes a plunger-like model that actually matches the rock core shape things and
S405: 1, Substituting the numerical values in the three-dimensional tensor data volume of nuclear magnetic resonance MRI/ _T2 to each node in the three-dimensional regular cubic network model (the connecting line value between two nodes is the hole/throat ), 2. For every hole/throat radius R in the network, take the mean of the MRI/ _T binaries at the two neighboring nodes; 3. Coordinate each node in the model to the spherical surface by randomly moving within the inner space, creating an irregular network space structure and causing random variations in the hole/throat length;4; Obtaining pore network models with different connectivity features by randomly removing them.

Step 5 is specifically the following steps: 1. Assuming that the initial value of the conversion factor is α, calculate the hole-throat radius R _i according to the method of step 4 and through the conversion factor α. After obtaining the digitized pore network model of the rock core by and then verifying that the simulated permeability agrees with the measured permeability of the rock core; recalculate the hole-throat radius R _i of the pore network, adjusting the conversion factor α, until it essentially agrees with the actual measured permeability of the rock core, and the digitized rock core and obtaining a digitized rock core pore network model corresponding to the actual rock core sample by calculating its permeability value.

Said step 6 specifically includes the following substeps.

S601: Analyzing the unsteady flow of fluid in the hole-throat channel based on the rock core-sized pore network coupled unsteady flow model to analyze the fluid axial velocity v _r distribution in the circular tube bundle;
S602: Constructing a non-stationary single-phase liquid fluid-structure interaction (FSI) infiltration mathematical model employing the finite volume method and the Renault transport equation, based on the hypothetical conditions to which the single-phase liquid infiltration process should fit;
S603: Unsteady single- _phase gas-fluid- _structure interaction ( FSI) building a mathematical model of penetration,
S604: Utilizing unsteady single-phase liquid fluid-structure interaction (FSI) seepage mathematical model and unsteady single-phase gas fluid-structure interaction (FSI) seepage mathematical model, and characterizing mixed fluid in void network model Constructing a fluid-structure interaction (FSI) multiphase seepage mathematical model under a process of unmixed exclusion, binding parameters;
S605: Combining the digitized rock core pore network model with the unsteady fluid-structure interaction (FSI) multi-phase infiltration numerical simulation method to perform the rock core fluid-structure interaction (FSI) infiltration numerical simulation. including one step.

で示す通りである（式の中で、C_t=C_ρ+C_pであり、C_tが総合圧縮係数であり、C_pが空隙圧縮係数であり（単位：Pa^-1）、C_ρが液体圧縮係数であり（単位：Pa^-1）、Δtが時間ステップであり、q_ijが管束内流体の体積流量であり、g_ijが管束内液体の水力伝導率であり、<p>_ijが隣のノードiとjとの間の平均圧力であり、Δp_ijが隣のノードiとjとの間の管束内流体圧力差であり、V_pi0が初期時刻におけるノードiの空隙体積であり、R_0ijが初期時刻におけるノードiとjとの間の穴・スロート管束の半径であり、l_ijが隣のノードiとjとの間の穴・スロート長さであり、nが制御ボリュームの中心ノードiにつながって続いているノード数であり、μが流体粘度である）。 Specifically, the unsteady single-phase liquid fluid-structure interaction (FSI) seepage mathematical model is expressed as follows:

(In the formula, C _t = C _ρ + C _p , C _t is the total compression coefficient, C _p is the pore compression coefficient (unit: Pa ⁻¹ ), and C _ρ is is the liquid compression coefficient (unit: Pa ⁻¹ ), Δt is the time step, q _ij is the volumetric flow rate of the fluid in the tube bundle, g _ij is the hydraulic conductivity of the liquid in the tube bundle, and <p> _ij is is the average pressure between neighboring nodes i and j, _Δpij is the intra-bundle fluid pressure difference between neighboring nodes i and j, _Vpi0 is the void volume of node i at the initial time, R _0ij is the radius of the hole-throat bundle between nodes i and j at the initial time, l _ij is the hole-throat length between adjacent nodes i and j, and n is the center of the control volume is the number of nodes connected to and following node i, and μ is the fluid viscosity).

で示す通りである（式の中で、g_gijが管束内ガス流動の水力伝導率であり、μ_gがガス粘度であり、p/Zμ_gが非線形項であり、pが制御ボリューム中心ノード位置の空隙率流体圧力であり、R_gがガス常数である）。 Specifically, the unsteady single-phase gas-fluid-structure interaction (FSI) seepage mathematical model is expressed as follows:

(where g _gij is the hydraulic conductivity of the gas flow in the tube bundle, μ _g is the gas viscosity, p/Z μ _g is the nonlinear term, and p is the control volume center node position is the porosity fluid pressure and R _g is the gas constant).

で示す通りである（式の中で、g_ijが管束内流体の水力伝導率であり、p_cijがノードiとjとの間の穴・スロートチャンネルの毛管力であり、μ_effがノードiとjとの間の混合流体の有効粘度であり、C_IとC_Dがそれぞれ注入流体と被排除流体の圧縮率因子であり、S_IとS_Dがそれぞれiを中心とする制御ボリューム内の注入流体と被排除流体の飽和度である）。 The fluid-structure interaction (FSI) seepage mathematical model is specifically expressed as follows:

(where g _ij is the hydraulic conductivity of the fluid in the tube bundle, p _cij is the capillary force in the hole/throat channel between nodes i and j, and μ _eff is the node i is the effective viscosity of the mixed fluid between and j, C _I and C _D are the compressibility factors of the injected and ejected fluids, respectively, and S _I and S _D are, respectively, the values in the control volume centered at i. is the saturation of the injected and rejected fluids).

The unsteady fluid-structure interaction (FSI) multiphase seepage numerical simulation method specifically includes the following steps.

線形化後の流体－構造連成（FSI）多相浸透数学モデルを取得すること、

S702：陰的数値シミュレーション手法を採用して、ソース・シンク項Q_iを導入して線形化後の流体－構造連成（FSI）多相浸透数学モデルに対して離散化処理を行った後、

を得ること（式の中で、上付き文字tが現在時刻の下での流動状態を表し、上付き文字t+Δtが次の時刻の流動状態を表す）、
S703：ステップ7023公式における異なる時刻の流動状態について分離と合併を行うことによって、

を得ること、
S704：空隙ネットワークモデルにおけるすべてのネットワーク・ノードを走査した後、空隙ネットワークモデルにおけるすべての制御ボリュームが、皆現在の時刻と次の時刻の流動状態を、ステップS703の公式に代入して整理してから、下記のマトリックス、即ち、 [A]^t+Δt[X]^t+Δt=[B]^tを形成した後（式の中で、[A]^t+Δtが、N×Nサイズの流体水力伝導率関連スパースマトリックスであり、Nが、空隙ネットワークモデルのノード数であり、[X]^t+Δtと[B]^tが、長さがNである二つのベクトルであり、[X]^t+Δtが次の時刻の圧力場ベクトルであり、[B]^tが前の時刻の圧力場と境界条件関連ベクトルである）、GPU代数多重グリッド一般化最小残差アルゴリズムを利用して、以上のマトリックスの解を求めることによって、現在時刻の空隙ネットワークモデルにおける流体圧力場分布を得ること、
S705：非定常流体－構造連成（FSI）多相浸透数値シミュレーションのプロセス中において、固定時間ステップ又は可変時間ステップの方式を採用することによって、当該ステップ幅内の流体界面に変位を発生させ、各時間ステップの下で、界面移動後の新しい流体界面位置を計算し、且つ、全体的な空隙ネットワークモデルにおけるすべての穴・スロートチャンネルの水力伝導率とモデルにおける各相流体の飽和度をアップデートした後、全体的な空隙ネットワーク・スペースに浸入流体がいっぱいになるか、又はある飽和度数値に達するまで、流体圧力場分布の解を求めることという五つのステップを含むこと。 S701: First, by linearizing the nonlinear capillary pressure term in the fluid-structure interaction (FSI) multiphase seepage mathematical model,

obtaining a linearized fluid-structure interaction (FSI) multiphase infiltration mathematical model;

S702: Adopting an implicit numerical simulation method to introduce the source-sink term Q _i to perform discretization on the linearized fluid-structure interaction (FSI) multiphase seepage mathematical model,

to obtain (in the formula, the superscript t represents the flow state under the current time, and the superscript t+Δt represents the flow state at the next time),
S703: By separating and merging the flow states at different times in the step 7023 formula,

to obtain
S704: After scanning all the network nodes in the pore network model, all the control volumes in the pore network model are arranged by substituting the current and next flow states into the formula of step S703. , after forming the following matrix: [A] ^t+Δt [X] ^t+Δt =[B] ^t (where [A] ^t+Δt is the N×N size fluid hydraulic is the conductivity-related sparse matrix, where N is the number of nodes in the void network model, [X] ^{t + Δt} and [B] ^t are two vectors of length N, and [X] ^{t +} where ^Δt is the pressure field vector at the next time and [B] ^t is the pressure field and boundary condition related vector at the previous time), we utilize the GPU algebraic multigrid generalized minimum residual algorithm to obtain the above matrix Obtaining the fluid pressure field distribution in the void network model at the current time by solving for
S705: During the process of unsteady fluid-structure interaction (FSI) multi-phase infiltration numerical simulation, by adopting the method of fixed time step or variable time step to generate displacement at the fluid interface within the step width, Under each time step, we calculated the new fluid interface position after interface movement and updated the hydraulic conductivity of all holes and throat channels in the global void network model and the saturation of each phase fluid in the model. Then, including five steps of solving the fluid pressure field distribution until the entire void network space is filled with the infiltrating fluid or some saturation value is reached.

The present invention relates to a digitized rock core pore network model and unsteady fluid-structure interaction (FSI) infiltration mathematical model simulation method based on nuclear magnetic resonance MRI/ _T2 , and thus provided in the present invention. The detailed realization process of the method is specifically divided into the following two flows.

1. To construct a pore network model of the digitized rock core based on nuclear magnetic resonance MRI/ _T2 .

(1) After selecting a rock core sample and drying it, measure the length, diameter, porosity and permeability of the rock core sample. A vacuum is applied to the rock core sample to simulate formation water or formation oil for saturation. Carry out nuclear magnetic resonance scanning to obtain nuclear magnetic resonance imaging MRI data and nuclear magnetic resonance _T2 spectrum data of corresponding rock core samples (for the scanning process, it can be realized by referring to nuclear magnetic resonance imaging , "Nuclear Magnetic Resonance Imaging", 2004, Higher Education Publishing House).Perform MRI measurements on different positions of the rock core and acquire MRI two-dimensional images of the rock core at the corresponding positions.Scanning of experimental instruments Depending on the accuracy and within the processing range of the experimental equipment, select the appropriate section position and number of sections to acquire nuclear magnetic resonance MRI two-dimensional images at different section positions of the rock core.The section position coordinates and pixel data are Save to TXT text.

(2) Acquire a three-dimensional tensor data volume of nuclear magnetic resonance MRI/ _T2 of the rock core. By interpolating the data volume in nuclear magnetic resonance imaging with the two-dimensional image TXT text and the rough range of the rock core hole and throat lengths acquired in step (1), the data volume scale is determined by the reaction rock meet the microscopic hole and throat feature requirements of Algorithms such as trilinear interpolation and Kriging interpolation can be used as the interpolation algorithm. The basic parameters of the interpolation are determined by the actual length, diameter and spatial position of the two-dimensional image slices of the rock core. Acquire the MRI/T ₂ three-dimensional tensor data volume through interpolation.

(3) Build a space irregular void network model. Set the size of the model, the number of ligands and the average hole/throat length according to the scale of the rock core. First, adopt the regular cube network structure of C++ language and matrix computing library Eigen structure to generate an X×Y×Z three-dimensional regular cube network (the values of X, Y and Z are , determined by the scale of the MRI/T2 three-dimensional tensor data volume obtained in the step). We set the total number of nodes in the network model to be (X-1×(Y-1)×(Z-1), each node represents one gap, and there is a throat channel ( In the constructed network model, six throats are connected to each other around each node representing a void, and the length of the throat is set to the average hole length of the rock. The throat length <l>, and the x-, y-, and z-direction side lengths of the model are L _x =(X-1)<l>, _Ly =(Y-1)<l>, and L _z =( Z-1) <l> Realize complete connections through round tubes between all grid nodes in the model (at this time, the number of arbitrary node ligands z = 6), and make gaps and throats by setting L _y and L _z to the actual rock core diameters, and removing all points in the yoz plane of each layer whose distance from the center point is greater than L _y . is a plunger-like model that matches the actual rock core geometry.Substitute the numerical values in the three-dimensional tensor data volume of nuclear magnetic resonance MRI/ _T2 into the nodes in the three-dimensional regular cubic network model ( The connecting line value between two nodes is the hole/throat radius R.) For every hole/throat radius R in the network, take the _average of the MRI/T binaries of the two adjacent nodes in the network. By randomly moving the coordinates of each node in the space within the sphere, an irregular network-space structure is generated, and random variations in hole-throat length occur.From the network structure, By randomly removing some connections, we obtain void network models with different connectivity (number of ligands) features.

(4) Digitized rock core pore network model. What is reflected by the nuclear magnetic resonance MRI/ _T2 data volume of rock core samples is the relative size value of the inter-rock pore space and does not reflect the size of the rock hole/throat radius. Therefore, a trial-and-error method is employed here to estimate the conversion factor (R _i =αA _i ) between the rock hole-throat radius R _i and the corresponding MRI/ _T2 data amplitude A _i , thereby obtaining the digital A pore network model of the refined rock core can be constructed. Obtain the pore network model of the digitized rock core by calculating the hole-throat radius R _i based on the method of step (3) and through the conversion factor α, assuming the initial values of the conversion factors. Then, a single-phase stable seepage pore network simulation algorithm is adopted to calculate the permeability in the digitized pore network model of the constructed rock core. Verify that it matches the measured value of If not, the pore volume is adjusted by adjusting the conversion factor α until the permeability in the digitized rock core pore network model essentially matches the actual measured rock core permeability. Re-calculate the network hole/throat radius R _i , re-build the pore network model of the digitized rock core, and calculate its permeability value. This yields a digitized rock core pore network model that corresponds to the actual rock core sample.

For a cylindrical rock core sample that is 2.5 cm in diameter and 5 cm in length, there are over 1 million nodes in the corresponding digitized pore network model of the rock core. In this case, it is difficult to solve such a problem with the general CPU sparse matrix equation solving algorithm, so we have to perform the calculation and solving with the GPU algorithm.

The digitized rock core network model constructed by the above method can apply the microscopic infiltration mechanism to the rock core scale porosity network to conduct numerical simulation studies, and It is possible to carry out comparative verification with the macroscopic experimental results obtained in the experiment, and to perform a larger-scale scale-up analysis on the basis of the laboratory rock core analysis, so that the laboratory rock core analysis and the macroscopic analysis can be performed. It is an important complement to scale oil and gas mineral numerical simulations, and builds a bridge between laboratory rock core analysis and macroscopic oil and gas mineral numerical simulations.

2. To implement the unsteady fluid-structure interaction (FSI) infiltration mathematical model and numerical simulation method applied to the pore network model of the digitized rock core.

On the basis of the constructed digitized pore network model of the rock core, the Reynolds stress-transport equation of fluid mechanics, the finite volume method, to correctly explain the multiphase fluid transfer process and the interaction relationship between the fluid and the rock framework. and the fluid-structure interaction (FSI) infiltration mechanism, derive the unsteady fluid-structure interaction (FSI) multiphase infiltration mathematical model applied to the pore network model of the digitized rock core, and A numerical simulation method is presented to form a suite of nuclear magnetic resonance MRI/ _T2 -based digitized rock core modeling and unsteady fluid-structure interaction (FSI) infiltration simulations.

であることを得る。
NguyenとChoi（Nguyen, Q.H., & Choi, S.B., A new approach for an analytical solution of unsteady laminar flow in dispensing processes. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2010, 224(6), 1231-1243.）の研究によって、円形管束内における流体の軸方向速度v_rの分布が、下記の通り、即ち、

であることを得る。
式の中で、Rが円形管束の半径であり、Δpが管束両端の圧力差であり、lが管束長さであり、μが流体粘度であり、tが流動時間であり、ρが流体密度である。式2を採用して、管束半径R=1ミクロンと10ミクロンである時に、それぞれニュートン流体（水、密度ρ=1000kg/m³、粘度μ=1mPa・s）の流動時の速度を計算すると、次の事項、即ち、時間tがそれぞれ7e-7秒と5e-5秒に達する時に、流体の速度断面が、安定流動時の放物線分布に達することを発見できる。式2によって、微視的な穴・スロートチャンネルにおいて、流体の流動がより速く安定状態に達することを発見できるので、ポアズイユの式を、多孔質媒体中における流体の非定常流動に相変わらず運用できることを表す。（NguyenとChoiは、式2でマイクロナノメートル級の穴・スロートチャンネルにおける流体の流動を分析しなかった） (1) Unsteady flow of fluid in hole-throat channel The unsteady flow model of the pore scale is based on the hypothesis that the pressure drop is due to inertia and fluid viscosity. Assuming that the flow is laminar, by simplifying the Navier-Stokes equations, the axial velocity v _r of the fluid in the circular tube bundle is given by:

get to be
Nguyen and Choi (Nguyen, QH, & Choi, SB, A new approach for an analytical solution of unsteady laminar flow in dispensing processes. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 2010, 224(6 ), 1231-1243.), the distribution of the axial velocity v _r of the fluid in a circular tube bundle is as follows:

get to be
In the formula, R is the radius of the circular tube bundle, Δp is the pressure difference across the tube bundle, l is the tube bundle length, μ is the fluid viscosity, t is the flow time, and ρ is the fluid density is. Adopting Equation 2 to calculate the velocity during flow of a Newtonian fluid (water, density ρ=1000 kg/m ³ , viscosity μ=1 mPa s) when the tube bundle radii are R=1 and 10 microns respectively, It can be found that the velocity cross-section of the fluid reaches a parabolic distribution in stable flow when the time t reaches 7e-7 s and 5e-5 s respectively. 2 allows us to discover that the fluid flow reaches a steady state faster in the microscopic hole-throat channel, so we show that Poiseuille's equation can still be applied to the unsteady flow of fluids in porous media. show. (Nguyen and Choi did not analyze fluid flow in micronanometer-scale hole-throat channels in Equation 2.)

(2) Unsteady single-phase liquid fluid-structure interaction (FSI) seepage mathematical model A finite volume method is adopted to draw conclusions. First, the control volume in the pore network refers to all the pore space contained under the condition that any node in the pore network is centered and bounded by the half point of all the tube bundles connected to it, i.e., the pore space If the total number of nodes in the network is N, there are N control volumes. A single-phase liquid infiltration process satisfies the following hypotheses. 1) the fluid cannot pass through the tube bundle boundary of the void network and reach outside the tube bundle, 2) the fluid is a continuous-phase and Newtonian fluid, and 3) the fluid mass-energy Neglecting transformations and 4) liquid and rock frameworks have constant compressibility. In fact, rocks and fluids have micro-compressibility in the seepage process, which makes the pressure propagation unstable and causes pressure wave propagation during the flow process. The following is the derivation process of the unsteady single-phase liquid fluid-structure interaction (FSI) seepage mathematical model.

式の中で、ρが流体密度であり（単位：kg/m³）、pが制御ボリューム中心ノード位置の空隙率流体圧力であり、単位がPaであり、Rがパイプ半径であり、V_pがパイプ体積であり（単位：m³）、C_pが空隙圧縮係数であり（単位：Pa^-1）、C_ρが液体圧縮係数である（単位：Pa^-1）。前記式に対して、変量を分離した後、下記を得る。

同じ道理で、穴・スロート半径Rが下記を満たす。

式の中で、Vp0、R₀、ρ_0は、それぞれ初期時刻（p₀圧力の下で）の空隙体積法、穴・スロート半径と液体密度である。マクローリン級数によって式4を展開した後、前の二項目を取ると、下記を得る。

Fluid density and rock pore space volume changes with pressure change satisfy the following equations (ignoring hole/throat length changes due to pressure change):

In the equation, ρ is the fluid density (unit: kg/m ³ ), p is the porosity fluid pressure at the control volume center node location, the unit is Pa, R is the pipe radius, and V _p is the pipe volume (in m ³ ), C _p is the void compression coefficient (in Pa ⁻¹ ), and C _ρ is the liquid compression coefficient (in Pa ⁻¹ ). For the above equation, after separating the variables, we have:

By the same reason, the hole/throat radius R satisfies the following.

In the equations, Vp0, _R0 , _{ρ0 are the void volume, hole-throat radius and liquid density at the initial time (under p0 pressure), respectively} . After expanding Equation 4 by the Maclaurin series, taking the previous two terms, we get

式の中で、Vが、任意のノードを中心とする制御ボリューム体積を表し、ρvdAが、任意の時刻における、流体が制御ボリュームに流入したり、又はこれらか流出したりする質量流量を表す。レイノルズ応力輸送方程式は、ソース・シンク項がない状況の下で、制御面経由正味出力の当該物理量の流量及び制御ボリューム内の流体品質の地元時間に対する変化率を表す。任意の制御ボリューム内において、

となる。
式の中で、4aと4bによって、

となる。 In a void network, the law of conservation of mass for flow within an arbitrary control volume is expressed by the Reynolds stress-transport equation.

In the equations, V represents the control volume volume centered at any node, and ρvdA represents the mass flow rate at which fluid enters or leaves the control volume at any time. The Reynolds stress-transport equation expresses the rate of change of the physical quantity of interest in the net output through the control surface and the fluid quality in the control volume with respect to local time in the absence of source-sink terms. Within any control volume,

becomes.
In the equation, by 4a and 4b,

becomes.

式の中で、l_ijが隣のノードiとjとの間の穴・スロート長さである。 Since both C _p and C _ρ are sufficiently small numbers, ignoring the higher order terms in Equation 7b, we get

In the equation, l _ij is the hole-throat length between neighboring nodes i and j.

式の中で、A_ijが隣のノードiとjとの間の穴・スロートの断面積である。パイプの変形と圧力変動によって定常層流速度分布が壊される可能性があるが、式2の分析結果から見ると、微視的な穴・スロートチャンネル内の流体の流動速度断面が定常放物線速度分布まで速く回復することを発見できる。従って、任意の穴・スロートチャンネルの流量は、下記を満たす。

式9aと9bから、下記導きだすことができる。

式の中で、<p>_ijが隣のノードiとjとの間の平均圧力であり、Δp_ijが隣のノードiとjとの間の管束内流体圧力差である。式6、8及び10によって、空隙ネットワークモデルの非定常単相液体流動数学モデルを得る。 Express the net flow rate of fluid in and out of the hole-throat channel in any direction within the control volume as follows:

In the equation, A _ij is the cross-sectional area of the hole/throat between adjacent nodes i and j. Although the pipe deformation and pressure fluctuations may destroy the steady laminar flow velocity distribution, the analysis results of Eq. can be found to recover quickly up to Therefore, the flow rate for any hole/throat channel satisfies the following:

From equations 9a and 9b it can be derived that:

In the equation, <p> _ij is the average pressure between neighboring nodes i and j, and Δp _ij is the intra-bundle fluid pressure difference between neighboring nodes i and j. Equations 6, 8 and 10 give the unsteady single-phase liquid flow mathematical model of the void network model.

式の中で、C_t=C_ρ+C_pであり（C_tが総合圧縮係数である）、方程式両端の流体密度が消去されて、Δtが時間ステップである。（式の中で、q_ijが管束内流体の体積流量であり、g_ijが管束内液体の水力伝導率である）。式11が拡散方程式であり、当該方程式は、空隙ネットワークモデルでの流体プロセス中の圧力拡散（伝播）を述べる。C_t=C_ρ+C_p=0である場合、式11がキルヒホッフの法則とラプラス方程式（安定単相浸透数学モデル）に退化する。

In the equation, C _t =C _ρ +C _p (where C _t is the overall compression factor), fluid densities at both ends of the equation are eliminated, and Δt is the time step. (In the equation, q _ij is the volumetric flow rate of the fluid in the tube bundle and g _ij is the hydraulic conductivity of the liquid in the tube bundle). Equation 11 is the diffusion equation, which describes the pressure diffusion (propagation) during the fluid process in the pore network model. If C _t =C _ρ +C _p =0, Equation 11 degenerates to Kirchhoff's law and Laplace's equation (stable single-phase seepage mathematical model).

式の中で、Mがガス分子のモル質量であり、Zがガス偏差因子であり、pが制御ボリューム中心ノードの圧力であり、R_gがガス常数であり、Tがケルビン温度である。空隙ネットワークモデルにおけるガスの流動プロセスは、相変わらずレイノルズ応力輸送方程式（式6）を満たす。この時、式6の左側第一項を、下記の関係式で表すことができる。

式4、5と12から、下記導きだすことができる。

その中で、

となる。 (3) Unsteady single-phase gas-fluid-structure interaction (FSI) seepage mathematical model For a real gas, the calculation formulas for gas density ρ _g and gas compressibility factor C _g are as follows.

In the equation, M is the molar mass of the gas molecule, Z is the gas deviation factor, p is the pressure at the control volume center node, R _g is the gas constant, and T is the temperature in Kelvin. The gas flow process in the pore network model still satisfies the Reynolds stress transport equation (equation 6). At this time, the first term on the left side of Equation 6 can be expressed by the following relational expression.

From equations 4, 5 and 12, the following can be derived.

among them,

becomes.

ガス流動について、レイノルズ応力輸送方程式（式6）の左側第二項目が下記の関係式を満たす。

この時、ガスが流れて制御面を通り抜ける時の正味流量方程式は、下記の通りである。

式6、式15と式17の連立方程式によって、空隙ネットワークに適用される非定常単相ガス流体-構造連成（FSI）浸透数学モデルの構造式を得る。 As before, compared to C _p , C _p C _ρ is a sufficiently small number, and after ignoring that term, Equation 14 gives:

For gas flow, the second term on the left side of the Reynolds stress-transport equation (Equation 6) satisfies the following relationship:

The net flow equation as the gas flows through the control surface is then:

The system of equations 6, 15 and 17 gives the structural formula of the unsteady single-phase gas-fluid-structure interaction (FSI) seepage mathematical model applied to the pore network.

式18における非線形項p/Zμ_gに対する処理は、下記の通りである。低圧条件（ガス圧力が10MPaを下回る）の下で、ガスが理想気体であると仮定すると、この時にガス偏差因子Zが1である。ガス粘度μ_gが、圧力から僅か小さい影響を受けて、μ_gが常数であると考えても良い。数値解法を採用する場合、式18の中の圧力pについて、前の時間ステップの圧力値を、直接に代入して、次の時間ステップの圧力場を計算することができる。この時に、式18の両側の圧力を消去することができる。高圧条件（ガス圧力が10Mpaを上回る）の下で、p/Zμ_gが常数に近似的に等しい。以上の二つの条件を満たす時に、式18を下記のように簡略化することができる。

The processing for the nonlinear term p/Zμ _g in Equation 18 is as follows. Assuming that the gas is an ideal gas under low-pressure conditions (gas pressure below 10 MPa), the gas deviation factor Z is 1 at this time. Gas viscosity μ _g is slightly affected by pressure, and μ _g may be considered a constant. When adopting the numerical solution, for the pressure p in Equation 18, the pressure value of the previous time step can be directly substituted to calculate the pressure field of the next time step. At this time, the pressures on both sides of Equation 18 can be eliminated. Under high pressure conditions (gas pressure above 10 Mpa), p/ _Zμg is approximately equal to a constant. When the above two conditions are satisfied, Equation 18 can be simplified as follows.

式の中で、g_gijが管束内ガス流動の水力伝導率である。この時、ガス圧縮率因子C_gが、岩石の空隙圧縮率因子C_pより遥かに大きいので、上式の右側項に空隙圧縮率因子が無視されている。式11と式19を比べる時に、非定常単相液体流体コウ構造連成（FSI）浸透（式11）と非定常単相ガス流体コウ構造連成（FSI）浸透（式18と式19）は、異なる制約条件の下で、同じ数学表現式を有することを発見できる。

In the equation, g _gij is the hydraulic conductivity of gas flow in the tube bundle. At this time, since the gas compressibility factor C _g is much larger than the pore compressibility factor C _p of the rock, the pore compressibility factor is ignored in the right term of the above equation. When comparing Eqs. 11 and 19, the unsteady single-phase liquid-fluid co-structural interaction (FSI) infiltration (Eq. 11) and the unsteady single-phase gas-fluid co-structural interaction (FSI) infiltration (Eqs. 18 and 19) are , have the same mathematical expression under different constraints.

Constraint Unsteady Fluid-Structure Interaction (FSI) Multiphase Infiltration Mathematical Model Multiphase here refers to both oil-water, gas-water and oil-gas phases. An oil-water phase elimination process will now be described as an example. The water-based oil removal process satisfies the following assumptions. 1) There is at most one fluid interface between different fluids in any hole or throat, 2) In the flow process, mixed phases and mixed phases do not occur in the process, 3) In holes and throats, 4) ignoring the effects of gravity;

ノードiを中心とする任意の制御ボリューム内において、混合流体の密度計算式は下記の通りである。

式の中で、ρ_Iが注入流体の密度であり、ρ_Dが被排除流体の密度であり、V_ijIとV_ijDがそれぞれノードiとノードjとの間の穴・スロートチャンネル内における注入流体と被排除流体の体積であり、V_IとV_Dがノードiを中心とする制御ボリューム内における注入流体と被排除流体の体積である。計算プロセスを簡略化する為に、ここで、ノードiを中心とする制御ボリューム内における混合流体の密度が、ノードiとjとの間の穴・スロートチャンネル内における混合流体の密度に近いと仮定するので、制御ボリューム内における混合流体の有効密度ρ_effが次の通りであると定義する。 In the hole-throat channel between any neighboring nodes i and j during the exclusion process, the density calculation formula for the mixed fluid is:

In an arbitrary control volume centered at node i, the density calculation formula for the mixed fluid is as follows.

where ρ _I is the density of the injected fluid, ρ _D is the density of the rejected fluid, and V _ijI and V _ijD are the injected fluid in the hole/throat channel between node i and node j, respectively. and the volume of the excluded fluid, and V _I and V _D are the volumes of the injected and excluded fluids in the control volume centered at node i. To simplify the computational process, we now assume that the density of the mixed fluid in the control volume centered at node i is close to the density of the mixed fluid in the hole-throat channel between nodes i and j. So we define the effective density ρ _eff of the mixed fluid in the control volume as:

式の中で、S_IとS_Dがそれぞれiを中心とする制御ボリューム内の注入流体と被排除流体の飽和度である。制御ボリューム内に注入流体に満ちる時に、S_I=1となり、制御ボリューム内における流体密度がρ_eff=ρ_Iとなる。被排除流体に満ちる時に、制御ボリューム内における流体密度がρ_eff=ρ_Dとなる。空隙ネットワークモデルにおいて、任意の制御ボリューム内において、混合流体の圧縮性は、下記の関係式を満たす。

式の中で、C_effρが、混合流体の有効圧縮率因子であり、ρ_eff0が、初期条件（圧力p₀）の下での混合流体の密度であり、C_IとC_Dが、それぞれ注入流体と被排除流体の圧縮率因子である。非混合相排除プロセス中において、任意の時刻に、空隙ネットワークモデルにおける多相流体の流動プロセスは、皆レイノルズ応力輸送方程式を満たす。非定常単相流体（ガスと液体）流体－構造連成（FSI）浸透数学モデルの導き出しプロセスにおいて採用される方法と類似して、非混合相排除プロセスの下での流体－構造連成（FSI）多相浸透数学モデルを得ることができる。

In the equations, S _I and S _D are the saturation of injected and rejected fluids in the control volume centered at i, respectively. When the control volume is filled with injection fluid, S _I =1 and the fluid density in the control volume is ρ _eff =ρ _I . When filled with displaced fluid, the fluid density in the control volume is ρ _eff =ρ _D . In the void network model, the compressibility of mixed fluids satisfies the following relationship within an arbitrary control volume.

In the equation, C _effρ is the effective compressibility factor of the mixed fluid, ρ _eff0 is the density of the mixed fluid under the initial condition (pressure p ₀ ), and C _I and C _D are the injection is the compressibility factor of the fluid and the fluid to be displaced. During the immixed phase exclusion process, at any time, the multiphase fluid flow processes in the pore network model all satisfy the Reynolds stress-transport equation. Fluid-structure interaction (FSI) under a non-mixed phase exclusion process, analogous to the method employed in the derivation process of unsteady single-phase (gas and liquid) fluid-structure interaction (FSI) seepage mathematical models. ) a multiphase seepage mathematical model can be obtained.

式の中で、g_ijが管束内流体の水力伝導率であり、p_cijがノードiとjとの間の穴・スロートチャンネルの毛管力であり（当該穴・スロートチャンネルに、非混合相流体界面が存在する時に、p_cij=2γcosθ/R_ij,γが非混合相流体の間の界面張力であり、θが湿潤接触角であり、穴・スロートチャンネルが単相流体である時に、p_cij=0となる）、μ_effがノードiとjとの間の混合流体の有効粘度であり（μ_eff=μ_IX_ij+μ_D(1-X_ij)、X_ijが、穴・スロートチャンネル内における非混合相流体界面の相対位置であり（0≦X_ij≦1）、μ_Iが注入流体の粘度であり、μ_Dが被排除流体の粘度である。制御ボリューム内に注入流体がある時に、μ_eff=μ_Iとなる。制御ボリューム内に被排除流体しかない時に、μ_eff=μ_Dとなる。二相流体界面が存在するある穴・スロートチャンネルに、Δp_ij<p_cijとなる時に、相応の径・スロートの水力伝導率g_ij＝0となるので、この時、当該流体界面が「ロック」される。空隙ネットワークのすべての制御ボリューム内に、皆単相の流体（ガス、油又は水）しかない時、式23は、式11又は式19に退化する。

In the equation, g _ij is the hydraulic conductivity of the fluid in the tube bundle, and p _cij is the capillary force in the hole/throat channel between nodes i and j (a non-mixed-phase fluid When an interface exists, p _cij =2γ cos θ/R _ij , where γ is the interfacial tension between non-mixed-phase fluids, θ is the wetting contact angle, and when the hole-throat channel is a single-phase fluid, p _cij = 0), μ _eff is the effective viscosity of the mixed fluid between nodes i and j (μ _eff =μ _I X _ij +μ _D (1-X _ij ), X _ij is the hole-throat channel is the relative position of the non-mixed-phase fluid interface in (0 ≤ X _ij ≤ 1), μ _I is the viscosity of the injected fluid, and μ _D is the viscosity of the excluded fluid, with the injected fluid in the control volume When μ _eff =μ _I. When there is only the fluid to be displaced in the control volume, μ _eff =μ _D. For a hole-throat channel with a two-phase fluid interface, Δp _ij <p _cij At this time, the corresponding fluid interface is “locked” because the corresponding diameter-throat hydraulic conductivity g _ij = 0. All single-phase fluids (gas, (oil or water), Equation 23 degenerates to Equation 11 or Equation 19.

Equation 23 applies to oil-water, gas-water and oil-gas two-phase infiltration. During a multiphase infiltration process, when one phase fluid is gas, the gas compressibility factor C _g is calculated by Equation 13. Assuming the gas to be an ideal gas, the gas compressibility factor C _g =1/p.

(5) Unsteady fluid-structure interaction (FSI) multiphase infiltration numerical simulation method The fluid-structure interaction (FSI) multiphase infiltration mathematical model (Eq. ) can be viewed as a further push. Here, the numerical simulation method will be described using Equation 23 as an example. First, the nonlinear capillary pressure term in Equation 23 is linearized.

これで、式23aを下記の式に書き直す。

式25と式23が完全に同等であるが、数学的モデルの形式が非定常単相浸透数学モデル（式11と19）に一層近い。陰的数値シミュレーション手法を採用して、ソース・シンク項Q_iを導入して（あるノード内に外部流体注入がある可能性があるか、又は流体が空隙ネットワークモデルから流出する可能性があることを表す。さもなければ、Q_i＝0となる）、式25に対して離散処理を行う。

Now rewrite equation 23a as:

Although Equations 25 and 23 are completely equivalent, the form of the mathematical model is closer to the unsteady single-phase infiltration mathematical model (Equations 11 and 19). An implicit numerical simulation approach is adopted to introduce a source-sink term _Qi (that there may be an external fluid injection within a node or that fluid may flow out of the pore network model). (otherwise Q _i =0) and perform discrete processing on Equation 25.

式の中で、上付き文字tが現在時刻の下での流動状態を表し（流体水力伝導率と圧力）、上付き文字t+Δtが次の時刻の流動状態を表す。式26における異なる時刻の流動状態について分離と合併を行うことによって（現在の流動状態量を、方程式の右側に置き、次の時刻の流動状態量を、方程式の左側に置く）、下記を得る。

In the equation, the superscript t represents the flow state under the current time (fluid hydraulic conductivity and pressure), and the superscript t+Δt represents the flow state at the next time. By separating and merging the flow states at different times in Eq.

After scanning all network nodes in the pore network model, all control volumes in the pore network model are similar to Eq. A system of equations consisting of N equations can be obtained. The system of equations subscripts allow the system of equations to be reduced to matrix form [A] ^t+Δt [X] ^t+Δt =[B] ^t .

その中で、

式の中で、[A]^t+Δtが、N×Nサイズの流体水力伝導率関連スパースマトリックスであり（Nが、空隙ネットワークモデルのノード数である）、[X]^t+Δtと[B]^tが、長さがNである二つのベクトルであり、[X]^t+Δtが次の時刻の圧力場ベクトルであり、[B]^tが前の時刻の圧力場と境界条件関連ベクトルである。以上のマトリックス方程式の解を求めることによって、現在時刻の空隙ネットワークモデルにおける流体圧力場分布を得ることができる。

among them,

In the equation, [A] ^{t + Δt} is the fluid hydraulic conductivity-related sparse matrix of N × N size (N is the number of nodes in the pore network model), and [X] ^{t + Δt} and [B ] ^t are two vectors of length N, [X] ^t+Δt is the pressure field vector at the next time, and [B] ^t is the pressure field and boundary condition related vector at the previous time. be. The fluid pressure field distribution in the void network model at the current time can be obtained by obtaining the solution of the above matrix equation.

During the flow process, the interface between the two phases is driven with time, so by choosing a time step, all the appropriate displacements Δx=Δt(q _ij /A _ij ), where A _ij is the cross-sectional area of the hole/throat channel. "Appropriate amount" here means that the number of operations must be reduced as much as possible under the premise of guaranteeing accuracy. Here, a fixed time step can be employed. A fixed time step refers to dividing the seepage transient into equal length time intervals Δt _k (k=1, 2, . . . ). When Δt _k is small enough, it can be considered that the pressure field change in the next time step has stability and linearity. At this time, the pressure field can be solved in partly the same way as for single-phase flow.

In the simulation process, a variable time-step approach can be employed to speed up the numerical simulation calculations. The elimination process is carried out until the following conditions are met. Under each time step, we compute the new fluid interface position after the interface movement and, when the entire void network space is filled with the invading fluid or reaches some saturation value, all in the stop type After updating the hydraulic conductivity of the hole and throat channels and the saturation of each phase fluid in the model, the pressure field is solved. During the simulation process, it is necessary to update the hydraulic conductivity after fluid interface migration occurs, and in the worst case, the fluid interface is "locked", so during the solution process, one severely pathological matrix (Eq. 27) may occur, greatly increasing the computational difficulty. During the simulation process, the above matrix is solved using the GPU Algebraic Multigrid Generalized Minimum Residual Algorithm (AMG-GMRES) in the AMGCL algorithm library. The flow rate Q remains unchanged during the entire simulation process.

By combining the digitized rock core pore network model with the fluid-structure interaction (FSI) infiltration numerical simulation method, the present invention can simulate and reproduce the process of laboratory rock core flow experiments. Digitalized single-phase and multi-phase (oil-water, gas-water and oil-gas) fluid-structure interaction (FSI) infiltration numerical simulation analysis and testing in rock core experiments can be realized. Multiphase infiltration simulation and analysis by combining the fluid-structure interaction (FSI) infiltration numerical simulation method presented in this invention with related methods such as the digitized rock core network model based on CT scans be able to.

In the simulation technique proposed in the present invention, the two-phase fluid interfacial movement is calculated explicitly after calculating the pressure field distribution at the current time in the void network model by an implicit method. In the digitized rock core pore network model construction and multiphase infiltration simulation method presented in this invention, by improving the model scale, accuracy and computational efficiency through GPU computational acceleration technology, based on the same infiltration mathematical model From pore network, indoor rock core and physical model experiments, a multi-scale leap of the shaft-single-well oil mineral model scale can be realized.

The invention has had the following beneficial effects.

1. The pore network of digitized rock core samples in laboratory experiments by adopting the data of MRI and T2 spectra in high-precision nuclear magnetic resonance imaging and combining the interpolation algorithm and the random irregular pore network model. By constructing the model, the model overcomes the problem that the scale of early digital rock cores is small, which is disadvantageous to the implementation of multiphase infiltration analysis. To more fully analyze the impact of rock microscopic pore/throat inhomogeneity and macroscopic inhomogeneity on infiltration processes, with the advantage of being directly proportional to the size and characteristics of the sample. can help improve the research capacity of oil and natural gas extraction rate solutions,
2. By analyzing the unsteady flow characteristics of fluids in a single microscopic circular channel, we prove that the Poiseuille equation and the parabolic velocity profile still apply to the pore-scale unsteady seepage process. and
3. Obtained a single-phase liquid-fluid-structure interaction (FSI) infiltration mathematical model by analyzing the micro-compressibility features of the liquid and rock framework and their interaction relationships in the single-phase liquid infiltration model. Furthermore, in the single-phase gas infiltration model, the micro-compressibility of the rock framework and the strong compressibility of the gas are simultaneously considered, and based on the coupled action relationship between gas compression and solid rock framework shape change, single-phase gas infiltration model further refine the pore-scale fluid-structure interaction (FSI) infiltration theory and method by obtaining a phase gas-fluid-structure interaction (FSI) infiltration mathematical model;
4. In the pore network model, the multiphase (oil-water, gas-water and oil-gas) infiltration mathematics is Furthermore, by introducing fluid-structure interaction (FSI) into the pore network multiphase infiltration mathematical model (considering the micro-compressibility of liquid and rock framework and the compressibility of gas ), by constructing an unsteady fluid-structure interaction (FSI) multiphase seepage mathematical model and deriving a corresponding implicit numerical simulation method for the mathematical model, the pore-scale fluid-structure interaction (FSI) multi-phase Further improvement of phase infiltration theory and numerical simulation method,
5. By combining the digitized rock core pore network model with the unsteady fluid-structure interaction (FSI) multiphase infiltration numerical simulation method, actual experiments of rock core flow (exclusion, self-absorption, etc.) in laboratory experiments The process and results can be simulated and reproduced, and in addition, the veracity and reliability of the simulation results are ensured by directly comparing and analyzing the analysis results of magnetic resonance MRI images of numerical simulations and rock core flow experiments.
6. Through GPU algorithms, the model building speed and multiphase infiltration simulation speed can be greatly improved, so that infiltration simulation studies can be conducted on the same rock sample under different displacement velocities (or pressures), and indoor Compared with rock core experiments, it has advantages such as short cycle, fast speed, convenient operation and accurate results, and is an important supplement for indoor rock core flow experiments.
7. Adopting the modeling and simulation method proposed in this invention for some indoor rock core flow experiment analysis, which must be carried out under high temperature and high pressure, and the operation is relatively complicated. Then, based on temperature and pressure, only by adjusting the rock framework and fluid related parameters, it is possible to realize multi-phase infiltration simulation research and analysis under high temperature and high pressure, and to obtain rock core flow analysis results conveniently and quickly. what you can get and
8. Similarly, combine the unsteady fluid-structure interaction (FSI) infiltration numerical simulation method presented in this invention with related methods such as the digitized rock core/void network model based on CT scans. to be able to perform multiphase infiltration simulations and analyses.

4 is a technical flow chart of the present invention; 1 is a flow chart of a nuclear magnetic resonance MRI scan; FIG. 2 is an illustration of nuclear magnetic resonance MRI two-dimensional images at different positions of a rock core; FIG. 2 is a schematic diagram of a visualized image of three-dimensional tensor data based on nuclear magnetic resonance MRI/ _T2 . 1 is a schematic diagram of a method for constructing a space irregular network; FIG. Schematic representation of a pore network model of a digitized rock core. FIG. 1 is an illustration of a two-phase exclusion process in a digitized rock core model; FIG. 2 is an illustration 2 of a two-phase exclusion process in a digitized rock core model; FIG. 3 is an illustration 3 of a two-phase exclusion process in a digitized rock core model; FIG. 4 is an illustrative diagram 4 of the two-phase exclusion process in a digitized rock core model; FIG. 5 is an illustrative diagram 5 of a two-phase exclusion process in a digitized rock core model;

In order to understand the technical features, objects and beneficial effects of the present invention more clearly, the technical solutions of the present invention are now described as follows. The described embodiments are clearly only some embodiments of the present invention, not all of them, and should not be taken as limitations on the possible scope of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in the field under the premise of not performing creative work are all within the protection scope of the present invention.

In this embodiment, as shown in Fig. 1, it is a digitized multi-phase fluid-structure interaction (FSI) infiltration numerical simulation method in an indoor rock core experiment, and the technical flow includes: After the rock core sample is first selected and dried, the hole/throat length, porosity and permeability of the rock core are measured.In addition, the rock core sample is evacuated to measure the saturation simulating formation water or formation crude oil; performing nuclear magnetic resonance MRI/T _2- scan measurements on rock core samples to obtain two-dimensional images of the rock core at different cross-sections; to obtain the three-dimensional data volume of the rock core MRI/ _T2 , build a three-dimensional irregular pore network model, and convert the MRI/ _T2 three-dimensional data volume into the three-dimensional irregular pore network model. By substituting into the node, through the conversion factor α, calculate the hole-throat radius of each neighboring node in the three-dimensional irregular pore network model, and simulate the permeability of the digitized pore network model of the rock core. and additionally adjusting the magnitude of the conversion factor α to ensure that the permeability of the pore network is close to the measured permeability of the rock core; Building a pore network model and combining the digitized pore network model of the rock core with the unsteady fluid-structure interaction (FSI) multiphase seepage numerical simulation method for rock core fluid-structure interaction (FSI) seepage Including six steps of performing numerical simulations.

In this embodiment, the specific implementation process of the method is as follows.

1. Digitized pore network model of rock core based on nuclear magnetic resonance MRI/T ₂ (1) Two-dimensional images of rock core are acquired by nuclear magnetic resonance MRI/T ₂ data.
The nuclear magnetic resonance measurement in the present invention adopts a MesoMR23-60H-I type nuclear magnetic resonance apparatus, and adopts a reverse pulse sequence and a Carr-Purcell-Meiboom-Gill pulse sequence to measure nuclear magnetic resonance signals. . A SIRT (simultaneous iterative reconstruction technique) inversion algorithm is used to generate the _T2 distribution. A low magnetic field nuclear magnetic resonance rock core analysis system (MesoMR-060H-HTHP-I) with a magnetic field strength of 0.5t is adopted for measurement. Main test parameters include main frequency (21.326MHz), echo interval (TE=0.2ms), polarization time step (TW=3000ms) and number of echoes (NECH=8000). In the rock core scanning process, by placing a standard rock core with a length of 5 cm and a diameter of 2.5 cm in a steady magnetic field and applying gradient magnetic fields in all three directions of x, y, and z, , gradient magnetic field strength=magnetic field strength difference at both ends of gradient magnetic field/length of gradient magnetic field. When the sample is taken, we first match the resonances in the layer and apply a phase code gradient to the magnetic field, then remove the phase code gradient and then apply a frequency code gradient to generate one symbol per voxel. write. This process is called coding or space positioning. After applying a radio frequency pulse to a layer, the MR signal of that layer is received. From this, after obtaining the size of each voxel MR signal in the layer by decoding, according to the correspondence with the code of each voxel in the layer, the size of the voxel signal is displayed in the corresponding pixel of the screen. . A specific scanning flow chart for nuclear magnetic resonance imaging is shown in FIG. The magnitude of the signal is displayed in different grayscale grades, with the greater the signal the greater the pixel brightness, while the smaller the signal the less pixel brightness. Two-dimensional rock core MRI images at different positions can be acquired by performing section processing on different positions of the rock core. Depending on the scanning accuracy of the experimental instrument, an appropriate section position and number of sections (e.g., 6 sections, Fig. 3) can be selected within the processing range of the experimental instrument to obtain different sections of the rock core, as shown in Fig. 4. Acquire nuclear magnetic resonance MRI two-dimensional images at the location. Save slice position coordinates and pixel data in TXT text. The pixel point spacing within a single two-dimensional photograph and the spacing between two two-dimensional photographs are determined by the resolution of the scanner.

(2) Acquisition of nuclear magnetic resonance MRI/T2 three-dimensional tensor data volume of rock core Under normal conditions, the rock pit and throat length/variation range in actual deposits is 50-300 microns, and the average pore・Throat length/about 100-150 microns (Bernabe, Y., Li, M., Tang, YB, & Evans, B., Pore space connectivity and the transport properties of rocks. Oil & Gas Science and Technology, 2016, 71(4), 50). By interpolating the data volume in magnetic resonance imaging by the rough range of the rock core hole and throat lengths, the data volume size is matched to the requirements of the microscopic hole and throat features of the reaction rock. make it fit. Algorithms such as trilinear interpolation and Kriging interpolation can be used as the interpolation algorithm. The basic parameters of the interpolation are determined by the actual length of the rock core and the spatial position of the 2D image slices. The MRI/T ₂ 3D tensor data volume acquired after interpolation essentially satisfies the data scale for building a digitized rock pore network model (as shown in Figure 14).

(3) Construction of space irregular structure pore network model Based on the scale of the rock core, set the size of the model, the number of ligands, and the average hole/throat length. First, adopting the C++ language and the regular cubic network structure of the matrix computing library Eigen structure to generate an X×Y×Z three-dimensional regular cubic network. From this, we set the total number of nodes in the network model to be (X-1*(Y-1)*(Z-1), each node representing one gap, and the throat between nodes. (a uniform circular pipe) connected through (a model two-dimensional cross section is shown in Figure 6) and the other part is filled with solid particulate matter.In the network constructed here, around each node representing a void, , the six throats are connected to each other, the throat length is set to be a constant l=150 μm, the throat radius is set to be unit 1, and the side lengths in the x, y, and z directions of the model are respectively Let L _x = (X-1) x l, _Ly = (Y-1) x l and L _z = (Z-1) x l Record the coordinates of each node in the network model. Realize perfect connections through round tubes between all grid nodes in (at this time, the number of arbitrary node ligands z=6), and set the radius ratio of the gap to the throat to 1. _Ly and L _z and remove all points in the yoz plane of each layer whose distance from the center point is greater than 0.5L _y , so that the model is a plunger-like model that matches the actual rock core shape Substitute the numerical values in the three-dimensional tensor data volume of nuclear magnetic resonance MRI/ _T2 into the nodes in the three-dimensional regular cubic network model (the connecting line value between two nodes is the hole/throat radius R ).For every hole/throat radius value in the network, take the average value of the numerical values at the two neighboring nodes.From this, the coordinates of each node in the network are randomly moved in space within the sphere by , generate an irregular network space structure.As shown in Fig. 5, by randomly removing some connections from the network structure, we obtain a void network model with different connectivity (number of ligands) features .

(4) Digitized pore network model of rock core With the above method, a digitized pore network model of rock core (Fig. 6) based on nuclear magnetic resonance MRI/T ₂ can be obtained. What is reflected by the MRI data volume of rock core samples is the relative size value of the inter-rock pore space and does not directly reflect the size of the rock hole/throat radius. Therefore, a trial-and-error method is employed here to estimate the rock hole/throat radius value and to estimate the rock hole/throat radius value by obtaining the conversion factor between the rock's actual hole/throat radius and the MRI data volume, and the pore network model. Can be used for construction. Obtaining the pore network model of the digitized rock core by assuming initial values of the transformation coefficients and calculating the hole-throat radius R _i through the transformation coefficients, and then applying the single-phase stable infiltration pore network simulation algorithm is employed to calculate the permeability in the constructed digitized pore network model of the rock core, and then the simulated permeability is verified to match the measured permeability of the rock core. If not, the hole-throat radius R By recalculating _i , rebuilding the digitized rock core pore network model, and calculating its permeability value, the digitized rock core pore space corresponding to the actual rock core sample Get the network model.

を満たす（gが水力伝導率であり、pが圧力である）。ラプラス方程式を、空隙ネットワークモデルに応用することによって、空隙ネットワークモデルにおける、流体定常浸透の時に質量保存の法則：Σ_jq_ij=0を満たすことを得る。この関係によって、空隙ネットワークにおけるすべてのノードを走査することによって、線性連立方程式又はスパースマトリックス方程式[A][X]=[B]を構築することができる。本研究において、共役勾配を採用して、前記マトリックス方程式の解を求めることによって、ネットワークモデルにおいて流体が流動する流体圧力場を得ることができる。これから入口端と出口端との圧力差を通して、モデルにおける流体の流入と流出の流量及び、モデルの浸透率を計算する。 Calculation method of permeability in porous network model (Tang, YB, Li, M., Bernabe, Y., & Zhao, JZ, Viscous fingering and preferential flow paths in heterogeneous porous media. Journal of Geophysical Research: Solid Earth, 2020, 125(3), e2019JB019306): According to Kirchhoff's law, the total flow rate of the inflow and outflow of the fluid in the node is zero. As the fluid moves in the pore network, the Laplace equation:

(g is the hydraulic conductivity and p is the pressure). By applying the Laplace equation to the pore network model, we obtain that the pore network model satisfies the law of conservation of mass: Σ _j q _ij =0 at steady fluid seepage. This relationship allows us to construct a linear system of equations or a sparse matrix equation [A][X]=[B] by scanning all the nodes in the void network. In this work, conjugate gradients can be adopted to obtain the fluid pressure field in which the fluid flows in the network model by solving the matrix equation. From this, through the pressure difference between the inlet and outlet ends, the inflow and outflow flow rates of the fluid in the model and the permeability of the model are calculated.

2. Fluid-structure interaction (FSI) infiltration numerical simulation based on the digitized rock core network model. At Q, the oil is injected into the model from the center entry point on the left edge of the rock core and the oil in the model is eliminated. After setting the pressure at the outlet end of the right end of the model to be atmospheric pressure 0.1 MPa, and then substituting the above condition into Equation 27, the unsteady fluid-structure interaction (FSI) multiphase infiltration proposed in the present invention A mathematical model and numerical simulation method are adopted to carry out simulation analysis of water-induced oil rejection in rock cores. Figures 7-11 are illustrative diagrams of the two-phase exclusion process in the digitized rock core model constructed by the present invention. shows an image that occupies the Among them, the red color is the pipe occupied by the invasion of the excluded phase fluid, and the excluded phase fluid is not displayed.

The foregoing has illustrated and described the basic principles, main features and advantages of the present invention. Those skilled in the art will appreciate the following: the present invention is not limited by the foregoing embodiments, and what is described in the foregoing embodiments and specification is the principle of the invention, and the Various changes and improvements may be made to the present invention without departing from the spirit and scope of the invention, and all such changes and improvements shall be included within the scope of the present invention for which protection is sought. The scope of the invention for which protection is sought is defined by the appended claims and their equivalents.

Claims (10)

as follows, i.e.
Step 1: First select the rock core sample, after drying, measure the hole/throat length, porosity and permeability of the rock core, additionally evacuate the rock core sample, simulating formation water or formation crude oil for saturation;
Step two: performing nuclear magnetic resonance MRI/T ₂ scan measurements on the rock core sample to obtain two-dimensional images of the rock core at different cross-sections;
Step 3: Performing interpolation on the 2D image to obtain a 3D data volume of the rock core MRI/ _T2 ,
Step 4: constructing a three-dimensional random void network model and substituting the MRI/T ₂ three-dimensional data volume into the nodes of the three-dimensional random void network model;
Step 5: Calculate the hole and throat radius of each neighboring node in the three-dimensional irregular pore network model through the conversion coefficient α, simulate the permeability of the digitized rock core pore network model, and In addition, by adjusting the magnitude of the conversion factor α, we ensure that the permeability of the pore network is close to the measured permeability of the rock core, and in addition, the digitized rock core pore network model and constructing
Step 6: Combining the digitized rock core pore network model with the unsteady fluid-structure interaction (FSI) multiphase infiltration numerical simulation method to perform the rock core fluid-structure interaction (FSI) infiltration numerical simulation. A digitized multiphase fluid-structure interaction (FSI) infiltration numerical simulation method for laboratory rock core experiments characterized by six steps. The above step 2 specifically includes the following: 1. Use a nuclear magnetic resonance apparatus to perform a nuclear magnetic resonance scan on the rock core sample after vacuum evacuation and saturation numerical simulation; ₂ , section processing for each location of the rock core; After obtaining 2D rock core nuclear magnetic resonance MRI images at different section positions of the rock core by performing a nuclear magnetic resonance scan on the rock core by selecting an appropriate section position and number of sections within the processable range of , saving the section position coordinates and the pixel data of the two-dimensional image into the two-dimensional image TXT text. Phase fluid-structure interaction (FSI) infiltration numerical simulation method. The step 3 specifically uses the following content, namely, the interpolation algorithm, in addition to the two-dimensional image TXT text obtained in step 2 and the rock core hole and throat length in step 1, By performing interpolation processing on the data volume in magnetic resonance imaging, that is, the nuclear magnetic resonance MRI two-dimensional image of the rock core, the data volume scale meets the requirements of the microscopic hole/throat characteristic surface of the reaction rock. , and thereby acquiring a three-dimensional tensor data volume of rock core MRI/ _T2 . Fluid-structure interaction (FSI) seepage numerical simulation method. wherein said step specifically includes the following substeps:
S401: Set the size of the model, the number of ligands and the average hole/throat length according to the scale of the rock core,
S402: Adopting the three-dimensional regular cubic network structure of the computer programming language and matrix computing library Eigen structure to generate an X×Y×Z three-dimensional regular cubic network;
S403: The total number of nodes in the cubic network is (X - 1 x (Y - 1) x (Z - 1), each node represents one gap, and the nodes are connected through the throat, be installed so that the other part is a rock framework;
S404: 1, in the constructed cubic network model, six throats are connected to each other around each node representing the void, the throat length is the rock average hole-throat length <l>, x of the cubic network model, Let the side lengths in the y and z directions be L _x =(X-1)<l>, _Ly =(Y-1)<l> and L _z =(Z-1)<l> respectively, and the cube Between all grid nodes in the network model, realization of perfect connections through round tubes, and the ratio of the radii of the _voids to the _throats of 1; and removing all points in the yoz plane of each layer whose distance from the center point is greater than _Ly , so that the cube network model becomes a plunger-like model that actually matches the rock core shape things and
S405: 1, Substituting the numerical values in the three-dimensional tensor data volume of nuclear magnetic resonance MRI/ _T2 to each node in the three-dimensional regular cubic network model (the connecting line value between two nodes is the hole/throat ), 2. For every hole/throat radius R in the network, take the mean of the MRI/ _T binaries at the two neighboring nodes; 3. Coordinate each node in the model to the spherical surface by randomly moving within the inner space, creating an irregular network space structure and causing random variations in the hole/throat length;4; Digitized multiphase fluid in laboratory rock core experiments according to claim 1, characterized in that it comprises four steps of obtaining pore network models with different connectivity features by randomly removing - Structural interaction (FSI) infiltration numerical simulation method. Step 5 is specifically the following steps: 1. Assuming that the initial value of the conversion factor is α, calculate the hole-throat radius R _i according to the method of step 4 and through the conversion factor α. After obtaining the digitized pore network model of the rock core by and then verifying that the simulated permeability matches the measured permeability of the rock core2, and if not, the permeability in the digitized rock core pore network model is verified. , recalculate the hole-throat radius R _i of the pore network, adjusting the conversion factor α, until it essentially matches the actual measured permeability of the rock core. including two steps of obtaining a digitized rock core pore network model corresponding to an actual rock core sample by rebuilding the pore network model and calculating its permeability values. A digitized multi-phase fluid-structure interaction (FSI) infiltration numerical simulation method for laboratory rock core experiments according to claim 1. Step 6 specifically includes the following substeps:
S601: Analyzing the unsteady flow of the fluid in the hole-throat channel based on the rock core-sized pore network coupled unsteady flow model to analyze the fluid axial velocity v _r distribution in the circular tube bundle;
S602: Constructing a non-stationary single-phase liquid fluid-structure interaction (FSI) infiltration mathematical model employing the finite volume method and the Renault transport equation, based on the hypothetical conditions to which the single-phase liquid infiltration process should fit;
S603: Unsteady single- _phase gas-fluid- _structure interaction ( FSI) building a mathematical model of penetration,
S604: Utilizing unsteady single-phase liquid fluid-structure interaction (FSI) seepage mathematical model and unsteady single-phase gas fluid-structure interaction (FSI) seepage mathematical model, and mixing fluid characteristic parameters in pore network model constructing a fluid-structure interaction (FSI) multiphase seepage mathematical model under the process of unmixed exclusion, in conjunction with
S605: Combining a digitized rock core pore network model with an unsteady fluid-structure interaction (FSI) multiphase infiltration numerical simulation method to perform a rock core fluid-structure interaction (FSI) infiltration numerical simulation. A digitized multi-phase fluid-structure interaction (FSI) seepage numerical simulation method for laboratory rock core experiments according to claim 1, characterized in that it comprises three steps. The unsteady single-phase liquid-fluid-structure interaction (FSI) seepage mathematical model has the following formula:

Digitized multi-phase fluid-structure interaction (FSI) seepage numerical simulation method for laboratory rock core experiments according to claim 6, characterized in that C _t =C _ρ + C _p , C _t is the overall compression coefficient, C _p is the void compression coefficient (unit: Pa ⁻¹ ), C _ρ is the liquid compression coefficient (unit: Pa ⁻¹ ), and Δt is is the time step, q _ij is the volumetric flow rate of the fluid in the bundle, g _ij is the hydraulic conductivity of the liquid in the bundle, <p> _ij is the average pressure between neighboring nodes i and j, _Δpij is the intra-bundle fluid pressure difference between neighboring nodes i and j, _Vpi0 is the void volume of node i at the initial time, and _R0ij is the hole between nodes i and j at the initial time is the radius of the throat tube bundle, l _ij is the hole-throat length between adjacent nodes i and j, n is the number of nodes connected to and following the central node i of the control volume, and μ is is the fluid viscosity). The unsteady single-phase gas-fluid-structure interaction (FSI) seepage mathematical model has the following formula:

Digitized multi-phase fluid-structure interaction (FSI) seepage numerical simulation method for laboratory rock core experiments according to claim 6, characterized in that g _gij is the tube bundle is the hydraulic conductivity of the internal gas flow, μ _g is the gas viscosity, p/Z μ _g is the nonlinear term, p is the porosity fluid pressure at the control volume center node location, and R _g is the gas constant ). The fluid-structure interaction (FSI) infiltration mathematical model has the following formula:

Digitized multi-phase fluid-structure interaction (FSI) seepage numerical simulation method for indoor rock core experiments according to claim 6, characterized in that g _ij is the tube bundle is the hydraulic conductivity of the internal fluid, p _cij is the capillary force in the hole-throat channel between nodes i and j, μ _eff is the effective viscosity of the mixed fluid between nodes i and j, and C _I and _CD are the compressibility factors of the injected and excluded fluids, respectively, and _SI and _SD are the saturations of the injected and excluded fluids, respectively, in the control volume centered at i). The unsteady fluid-structure interaction (FSI) multi-phase infiltration numerical simulation method specifically includes the following steps:
S701: First, by linearizing the nonlinear capillary pressure term in the fluid-structure interaction (FSI) multiphase seepage mathematical model,

obtaining a linearized fluid-structure interaction (FSI) multiphase infiltration mathematical model;

S702: Adopting an implicit numerical simulation method to introduce the source-sink term Q _i to perform discretization on the linearized fluid-structure interaction (FSI) multiphase seepage mathematical model,

to obtain (in the formula, the superscript t represents the flow state under the current time, and the superscript t+Δt represents the flow state at the next time),
S703: By separating and merging the flow states at different times in the step 7023 formula,

to obtain
S704: After scanning all the network nodes in the pore network model, all the control volumes in the pore network model are arranged by substituting the current and next flow states into the formula of step S703. from the following matrix:
After forming [A] ^{t + Δt} [X] ^{t + Δt} = [B] ^t (where [A] ^{t + Δt} is a fluid hydraulic conductivity-related sparse matrix of size N × N, where N is the number of nodes in the void network model, [X] ^t+Δt and [B] ^t are two vectors of length N, and [X] ^t+Δt is the pressure field at the next time vector, where [B] ^t is the pressure field and boundary condition related vector at the previous time), by solving the above matrix using the GPU algebraic multigrid generalized minimum residual algorithm, the current Obtaining the fluid pressure field distribution in the pore network model of time,
S705: During the process of unsteady fluid-structure interaction (FSI) multi-phase infiltration numerical simulation, by adopting the method of fixed time step or variable time step to generate displacement at the fluid interface within the step width, Under each time step, we calculated the new fluid interface position after interface movement and updated the hydraulic conductivity of all holes and throat channels in the global void network model and the saturation of each phase fluid in the model. 6. After solving the fluid pressure field distribution until the entire void network space is filled with the infiltrating fluid or reaches a certain saturation value. A digitized multiphase fluid-structure interaction (FSI) infiltration numerical simulation method for the laboratory rock core experiment described in .

Images