WO2021243765A1 - 基于gpu矩阵的离散元流固耦合数值模拟方法及系统 - Google Patents

基于gpu矩阵的离散元流固耦合数值模拟方法及系统 Download PDF

Info

Publication number
WO2021243765A1
WO2021243765A1 PCT/CN2020/097584 CN2020097584W WO2021243765A1 WO 2021243765 A1 WO2021243765 A1 WO 2021243765A1 CN 2020097584 W CN2020097584 W CN 2020097584W WO 2021243765 A1 WO2021243765 A1 WO 2021243765A1
Authority
WO
WIPO (PCT)
Prior art keywords
connection
fluid
matrix
pore
particle
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/CN2020/097584
Other languages
English (en)
French (fr)
Inventor
刘春�
朱遥
刘辉
乐天呈
施斌
王宝军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University
Nanjing Tech University
Original Assignee
Nanjing University
Nanjing Tech University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University, Nanjing Tech University filed Critical Nanjing University
Publication of WO2021243765A1 publication Critical patent/WO2021243765A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Definitions

  • the invention relates to a particle discrete element simulation, in particular to a fluid-solid coupling numerical simulation of a two-dimensional discrete accumulation body.
  • Discrete element can solve this kind of problem very well, but several commonly used methods at present have certain shortcomings, such as DEM-LBM calculation is too large, it is difficult to carry out large-scale numerical simulation calculation; DEM-SPH in stability and accuracy There are still shortcomings in the above; DEM-CFD has a relatively large amount of numerical simulation calculations on the pore scale.
  • a fluid network is generally created on the basis of a collection of particles to simulate a solid to simulate the flow of a fluid.
  • the fluid grid is based on the collection of particles that are in contact with each other and can form a circular closed system, and the pores of the particles are filled with fluid. Therefore, when solid particles migrate due to external forces, the pores in the system are also changing.
  • the present invention provides a discrete element fluid-structure coupling numerical simulation method and system based on a GPU matrix to realize the dynamic simulation of a large-scale fluid-structure coupling system.
  • a discrete element fluid-structure coupling numerical simulation method based on GPU matrix including the following steps:
  • Step 1 Establish a pore network model, including the establishment of a grid number matrix: establish a preliminary array of neighboring particles, and each particle is the center particle relative to its neighboring particles.
  • the neighbor matrix N is first obtained, and the neighbor matrix N is sorted in the row direction according to the connection angle, and converted into the connection matrix C.
  • Step 2 Establish a fluid unit transport model:
  • Step 2.1 set the initial fluid element state and boundary conditions in the initial model.
  • Step 2.2 Calculate the density ⁇ i of each fluid unit through the relationship between water pressure, temperature and density, where i represents the number of the fluid unit.
  • Step 2.4 calculate the pore throat permeability coefficient in the model: the pore throat permeability coefficient is related to the length and width of the channel, the pore throat permeability coefficient:
  • k(d) represents the pore throat permeability coefficient
  • k is the set value
  • L represents the distance between the centers of the two particles
  • R 1 and R 2 represent the radii of the two particles
  • minR represents the minimum radius of all particles
  • l represents Length of pore throat channel.
  • Step 2.5 calculate the percolation flow of a single fluid unit: within a time step dt of the discrete element calculation, a certain fluid unit pore i and the surrounding fluid unit pore n undergo fluid migration.
  • the fluid migration obeys Darcy's law, which defines the fluid inflow If it is positive, a single pore throat causes the volume seepage flow due to the pressure gradient:
  • P n and P i are the pressure fluid and pore i means the pore n
  • R a and R b are respectively composed of two particles of the pore throat radius of a and b
  • q i denotes seepage discharge
  • A represents the channel pore throat Width
  • J represents the pressure gradient coefficient
  • dT represents the time step.
  • Step 2.6 calculate the total transport volume:
  • the total transport volume of pore i with the surrounding fluid unit in one time step is:
  • the density of each fluid unit V i is the pore volume
  • M i is the mass of each fluid cell in the pores.
  • Step 2.7 calculate the mass of the fluid unit after migration: the total mass of pore i after the first time step:
  • M'i represents the mass of the fluid unit after fluid migration
  • M i represents the mass before migration
  • Q n represents the total migration volume
  • Step 2.8 update the state of the fluid unit: the density of the fluid unit after a time step dt
  • Step 3 fluid-solid coupling: pore water pressure acts on the particles, and the movement of the particles changes the volume of the pores, thereby achieving fluid-solid coupling.
  • Step 3.1 the solid particles are rebalanced: due to the movement of the fluid, the water pressure in each fluid unit changes, so that the force between each fluid unit on the solid particles also changes.
  • the particles are affected by the fluid unit, and the The displacement distance is:
  • x i is the displacement distance in a time step
  • P i is the pressure of the fluid cell
  • l i is the center of the particle and the fluid unit contact length
  • F i is the contact force between the particles
  • Step 3.2 Update pore volume: Since a change in position of the solid particles, the volume of the fluid cell has changed, V 'i is the volume of the pores of the solid particles after change.
  • Step 3.3 fluid unit status update: calculate the density of the fluid unit after the pore volume changes, and its density after a time step V '1 pore volume of the solid particles after change.
  • Step 3.4 repeat step 1 to step 3.3 to realize the two-dimensional fluid-structure coupling numerical simulation.
  • step 1 the method of establishing the pore network model in step 1:
  • Step 1.1 construct the initial particle discrete element stacking model, the total number of particles is m, and the particles are numbered continuously from 1 to m.
  • Step 1.2 in the discrete element calculation, the neighbor matrix N is first obtained.
  • the neighbor matrix N records the number of the neighboring particles around each numbered particle. Because the number of neighboring particle units of each particle unit is different, the total length is used for the areas where the row length is inconsistent. Add 1 to the number of particles to fill.
  • Step 1.3 input the maximum aperture Dmax to calculate the actual connected particles. Based on the neighbor matrix N, if the distance between the neighboring particles and the central particle is less than Dmax, it is recorded as 1, and it is judged to be connected. Again, it is 0, and it is judged as not connected, and the connection is the pore throat channel between the two particles.
  • Step 1.4 Regarding particles as particle units, a connection matrix C is established in the connection matrix according to the angle Na of the connection line between each particle unit and its neighboring particle units with respect to the positive x-axis direction, and the connection matrix C is arranged in order according to the angle of each connection.
  • Step 1.5 find the position index of the reverse connection, and rearrange the connection matrix C in the row direction so that the number in the first column is smaller than the number in the second column.
  • Step 1.6 in order to arrange the connection matrix C in sequence, multiply the number in the first column of the connection matrix C by a very large positive number plus the number in the second column, and then add the sum of the obtained numbers to each row from small to large
  • the bottom arrangement makes the paired connections close to each other, that is, the reverse connection is next to the positive connection, and the index of the reverse connection corresponding to each connection in the connection matrix C is found through index changes, and the reverse connection index matrix CI is established.
  • Step 1.7 find the position index of the next connection: to obtain the next connection of a certain connection of the pore, only need to get the index of the connection whose anti-connection angle is reduced, that is, the index number of the anti-connection of a certain connection of the pore minus 1.
  • the connection corresponding to the index is the next connection, and the next connection matrix CN corresponding to the connection matrix C is established, so that the position index of the next connection is successfully found.
  • Step 1.8 search for pores: starting from the initial connection, use the next connection index CN to search the particle unit connection forward, and back to the original particle unit to obtain a pore.
  • step 1.9 starting from the particle unit numbered 1, each connection uses this connection as the initial connection, and searching the particle unit connection forward until it returns to the original particle unit, thus successfully describing all the pores in the accumulation body.
  • Step 1.10 delete duplicate records, get all pores.
  • Step 1.11 build a pore network model based on the searched connections and pores.
  • a GPU matrix-based discrete element fluid-structure coupling numerical simulation fast calculation system including an input module, a pore network model module, a fluid unit migration model module, and a fluid-structure interaction module.
  • the input module is used to input particle-related data, and set the initial fluid unit state and boundary conditions in the initial model.
  • the pore network model module is used to establish a grid number matrix according to particle-related data; at the same time, it is used to establish a preliminary neighboring particle array, each particle is the central particle relative to its neighboring particle; in the discrete element calculation, the neighbor matrix N is first obtained , Sort the neighbor matrix N in the row direction according to the connection angle, and transform it into the connection matrix C; use the connection matrix C to control the search algorithm: find the inverse connection index matrix CI of the connection matrix C and the next connection index matrix CN; then Starting from the initial connection, the index matrix CN of the next connection is used to search the node forward until it returns to the original point, and all the pores of the accumulation body are obtained.
  • the fluid unit migration model module is used to calculate the density of each fluid unit through the relationship between water pressure, temperature and density; calculate the mass of each fluid unit; calculate the pore throat permeability coefficient in the model; calculate the seepage flow of a single fluid unit; Calculate the total migration; calculate the mass of the fluid unit after migration and update the status of the fluid unit.
  • the fluid-solid interaction module is used to calculate the displacement distance of the particles in a time step under the action of the fluid unit; update the pore volume; calculate the density of the fluid unit after the pore volume changes, and its density after a time step.
  • the present invention has the following beneficial effects:
  • the invention realizes that all discrete element fluid-structure coupling calculations are completed in the form of pure matrix calculation based on GPU, and cleverly utilizes the reverse connection and the next connection, and successfully proposes a discrete element pore search method, which successfully solves the problem of fluid-structure coupling.
  • the difficulty of searching for pores caused by the continuous change of pores in the process is conducive to the further development of discrete element simulation research on fluid-solid coupling.
  • a pore fluid unit is proposed to describe the pore fluid medium, which provides a fast and concise method for simulating the fluid-solid coupling of the pore medium.
  • the invention can realize the dynamic simulation of a large-scale fluid-structure coupling system, and significantly improves the calculation speed and the number of calculations of the discrete element fluid-structure coupling numerical simulation.
  • Figure 1 is a flow chart of establishing a pore network model.
  • Figure 2 shows the establishment of a fluid unit migration model.
  • Figure 3 is a schematic diagram of the fluid-solid coupling process.
  • Figure 4 shows the particle accumulation model
  • Figure 5 shows the pore network model
  • Figure 6 is a schematic diagram of a single pore connection.
  • Figure 7 is a schematic diagram of the pore throat model.
  • Figure 8 is a schematic diagram of fluid migration.
  • Figure 9 is a schematic diagram of the force on the central particle.
  • a discrete element fluid-structure coupling numerical simulation method based on GPU matrix including the following steps:
  • Step 1 Establish a pore network model, including the establishment of a grid number matrix, as shown in Figure 1: A preliminary array of neighboring particles is established, and each particle is the center particle relative to its neighboring particles.
  • the neighbor matrix N is first obtained, and the neighbor matrix N is sorted in the row direction according to the connection angle, and converted into the connection matrix C.
  • Step 10 Construct the initial particle discrete element stacking model, the total number of particles is m, and the particles are numbered continuously from 1 to m.
  • Step 11 the neighbor matrix N is first obtained.
  • the neighbor matrix N records the number of the neighboring particles around each numbered particle. Since the number of neighboring particle units of each particle unit is different, the total length is used for the areas where the row length is inconsistent. Add 1 to the number of particles to fill.
  • Step 12 Input the maximum aperture Dmax to calculate the actual connected particles. Based on the neighbor matrix N, if the distance between the neighboring particles and the central particle is less than Dmax, it is recorded as 1, and it is judged to be connected. Again, it is 0, and it is judged as not connected, and the connection is the pore throat channel between the two particles.
  • Step 13 Regarding the particles as particle units, in the connection matrix, a connection matrix C is established according to the angle Na of the connection line between each particle unit and its neighboring particle units relative to the positive x-axis direction, and the connection matrix C is arranged in order according to the angle of each connection. As shown in Table 1:
  • Na records the angle (counterclockwise) of the connection line between each unit and its neighbors relative to the positive direction of the x-axis. So far, the connections have been sorted, and the next step is to find the pores.
  • Step 14 Find the position index of the reverse connection, and rearrange the connection matrix C in the row direction so that the number in the first column is smaller than the number in the second column, as shown in Table 2:
  • connection matrix is rearranged in the row direction
  • Step 15 In order to arrange the connection matrix C in sequence, multiply the number in the first column of the connection matrix C by a very large positive number (the number in this method is 1e7) plus the number in the second column, and the obtained number The sum of each row is arranged downwards from small to large, so that the paired connections are next to each other, that is, the reverse connection is next to the positive connection, and the index of the reverse connection corresponding to each connection in the connection matrix C is found through the index change. Establish the anti-connection index matrix CI.
  • Step 16 Find the position index of the next connection: In order to describe a certain pore, all connections that form the pore must be obtained. Therefore, it is necessary to know that any connection corresponds to the next connection of a certain pore until the connection is returned, so that a pore is successfully described.
  • the next connection of the 1-2 connection is the 2-3 connection
  • the 2-3 connection is the reverse connection of the 1-2 connection.
  • the next connection matrix CN corresponding to the connection matrix C is established, so that the position index of the next connection is successfully found.
  • Step 17 search for pores: starting from the initial connection, use the next connection index CN to search the particle unit connection forward, and return to the original particle unit to obtain a pore.
  • each connection uses this connection as the initial connection, and the particle unit connection is searched forward until it returns to the original particle unit, thus successfully describing all the pores in the accumulation body.
  • Step 19 Delete duplicate records and get all pores.
  • Step 110 Establish a pore network model based on the searched connections and pores.
  • Step 2 Establish a fluid unit migration model.
  • Step 20 the initial model needs to set the initial fluid unit state (temperature, pressure) and boundary conditions.
  • Step 23 Calculate the pore throat permeability coefficient in the model: the pore throat permeability coefficient is related to the channel length and width, as shown in Figure 7, the pore throat permeability coefficient:
  • k(d) represents the pore throat permeability coefficient
  • k is the set value
  • L represents the distance between the centers of the two particles
  • R 1 and R 2 represent the radii of the two particles
  • minR represents the minimum radius of all particles
  • l represents Length of pore throat channel.
  • Step 24 Calculate the percolation flow of a single fluid unit: within a time step dt of the discrete element calculation, a certain fluid unit pore i and the surrounding fluid unit pore n (the number of P n is the number of fluid units adjacent to P i) occurs Fluid migration, as shown in Figure 8, fluid migration obeys Darcy’s law, which defines the fluid inflow as positive, and the volume seepage flow of a single pore throat due to the pressure gradient:
  • P n and P i are the pressure fluid and pore i means the pore n
  • R a and R b are respectively composed of two particles of the pore throat radius of a and b
  • q i denotes seepage discharge
  • A represents the channel pore throat Width
  • J represents the pressure gradient coefficient
  • dT represents the time step.
  • Step 25 Calculate the total transport volume:
  • the total transport volume of pore i with the surrounding fluid unit in one time step is:
  • the density of each fluid unit V i is the pore volume
  • M i is the mass of each fluid cell in the pores.
  • Step 26 calculate the mass of the fluid unit after the migration: the total mass of pore i after the first time step:
  • the density V i is the volume of the pores
  • M i is the mass of each fluid cell in the pores
  • M 'i represents the mass of the fluid means a fluid migration occurs
  • M i represents the mass before the migration
  • Q n is the total amount of migration.
  • Step 27 update the state of the fluid unit: the density of the fluid unit after a time step dt
  • Step 3 fluid-solid coupling, as shown in Figure 3: pore water pressure acts on the particles, and the movement of the particles changes the volume of the pores, so as to achieve fluid-solid coupling.
  • the interaction between the fluid unit and the solid unit is realized.
  • Step 30 the solid particles are rebalanced: due to the movement of the fluid, the water pressure in each fluid unit changes, so that the force between each fluid unit on the solid particles also changes, and the particles are acted on by the fluid unit, as shown in Figure 9,
  • the displacement distance in a time step is:
  • x i is the displacement distance in a time step
  • P i is the pressure of the fluid cell
  • l i is the center of the particle and the fluid unit contact length
  • F i is the contact force between the particles
  • Step 31 the update pore volume: Since a change in position of the solid particles, the volume of the fluid cell has changed, V 'i is the volume of the pores of the solid particles after change.
  • Step 32 Update the fluid unit status: Calculate the density of the fluid unit after the pore volume changes, and its density after a time step V '1 pore volume of the solid particles after change.
  • Step 33 Repeat steps 0 to 32 to realize the two-dimensional fluid-structure coupling numerical simulation.
  • the system constructed according to the foregoing GPU matrix-based discrete element fluid-structure coupling numerical simulation method includes an input module, a pore network model module, a fluid unit migration model module, and a fluid-structure interaction module, where:
  • the input module is used to input particle-related data, and set the initial fluid unit state and boundary conditions in the initial model.
  • the pore network model module is used to establish a grid number matrix according to particle-related data; at the same time, it is used to establish a preliminary neighboring particle array, each particle is the central particle relative to its neighboring particle; in the discrete element calculation, the neighbor matrix N is first obtained , Sort the neighbor matrix N in the row direction according to the connection angle, and transform it into the connection matrix C; use the connection matrix C to control the search algorithm: find the inverse connection index matrix CI of the connection matrix C and the next connection index matrix CN; then Starting from the initial connection, the index matrix CN of the next connection is used to search the node forward until it returns to the original point, and all the pores of the accumulation body are obtained.
  • the fluid unit migration model module is used to calculate the density of each fluid unit through the relationship between water pressure, temperature and density; calculate the mass of each fluid unit; calculate the pore throat permeability coefficient in the model; calculate the seepage flow of a single fluid unit; Calculate the total migration; calculate the mass of the fluid unit after migration and update the status of the fluid unit.
  • the fluid-solid interaction module is used to calculate the displacement distance of the particles in a time step under the action of the fluid unit; update the pore volume; calculate the density of the fluid unit after the pore volume changes, and its density after a time step.
  • the pore system is simplified to a network composed of individual particles (as shown in Figures 4 and 5) by simplifying each particle to a point position to establish a pore network model.
  • particle unit 1 1-2 connections, The 2-16 connection and the 16-1 connection form a triangle.
  • N the neighbor matrix N and the unit coordinates, obtain the matrix Na of the connection angle between each unit and the neighbor unit; start the connected unit from 1, the angle is the counterclockwise angle between the connection and the horizontal line, and start sorting from 0 degrees to generate a sorted Connect the array C.
  • the connection has a direction.
  • C first determine the anti-connection position index of each connection, such as getting the position of the anti-connection 2-1 connection of the 1-2 connection, and then determine the next connection in the counterclockwise direction of the anti-connection 2-1 connection. That is 2-16 connection. From the 2-16 connection and then looking forward, we get a 16-1 connection, and finally we get a triangle 1-2-16 in a clockwise direction.
  • use CI2 to record the reverse connection index (such as 2-1) of the C connection (such as 1-2), and use the CN to record the next connection in C (2-16). Starting from the set initial connection, use CN to calculate continuously to obtain triangles and quadrilaterals.
  • the repeated records are removed to obtain all the pores, and the pore network model is established; the fluid in the pores migrate through the pore throats, and the flow rate is determined by the nature of the rock and soil, the diameter of the pore throat, etc.; according to the volumetric strain and volumetric model of the pores
  • the pore water pressure is determined by the quantity; the pore water pressure acts on the particles, and the movement of the particles changes the volume of the pores, thereby achieving fluid-solid coupling.
  • the invention significantly improves the calculation speed of the two-dimensional fluid-solid coupling numerical simulation, provides a fast and concise method for simulating the fluid-solid coupling of a porous medium, and can realize the dynamic simulation of a large-scale fluid-solid coupling system.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于GPU矩阵的离散元流固耦合数值模拟方法及系统，包括建立孔隙网络模型。建立流体单元运移模型：孔隙内的流体单元通过颗粒单元之间的孔喉相连接，其渗透服从达西定律，其流量由岩土体性质、孔喉直径等确定，水压与温度以及密度有着一个函数关系式。流体固体耦合作用：孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用。本发明可实现大规模流固耦合系统的动态模拟，显著提高了离散元流固耦合数值模拟的计算速度以及计算数量。

Description

基于GPU矩阵的离散元流固耦合数值模拟方法及系统 技术领域
本发明涉及一种颗粒离散元模拟，尤其是对二维离散堆积体的流固耦合数值模拟。
背景技术
岩土工程、地质工程和能源工程中的大量重要的工程问题均涉及到复杂流固耦合过程，如隧道突水突泥、注浆过程、砂基液化、地面沉降、水力压裂等。由于有限元法或者无单元法等常规的数值模拟方法无法充分考虑流体—固体的相互作用以及固体颗粒迁移流失后造成的孔隙率以及渗流的渐进变化，所以很难真实模拟流固耦合过程中渗流以及渗流破坏规律。
离散元能很好的解决该类问题，但目前常用的几种方法都存在一定的缺陷，如DEM-LBM计算量过大，难以进行大规模的数值模拟计算；DEM-SPH在稳定性与精度上还有所欠缺；DEM-CFD在孔隙尺度上的数值模拟计算量较大。在进行流固耦合的数值模拟时，一般都是在颗粒集合模拟固体的基础上创建流体网络，来模拟流体的流动。流体网格是以相互接触且能够形成环向闭合系统的颗粒集合为依托，颗粒的孔隙间充满流体。所以当固体颗粒由于外力作用发生迁移时，系统内的孔隙也在发生变化，此时如何在流固耦合这个动态的过程中获取离散元堆积体中所有的孔隙成为了一大难题。同时，如何实现计算规模的提高以满足大尺度的流固耦合要求，也是目前亟待解决的问题。
发明内容
发明目的：为了克服现有技术中存在的不足，本发明提供一种基于GPU矩阵的离散元流固耦合数值模拟方法及系统，实现大规模流固耦合系统的动态模拟。
技术方案：为实现上述目的，本发明采用的技术方案为：
一种基于GPU矩阵的离散元流固耦合数值模拟方法，包括以下步骤：
步骤1，建立孔隙网络模型，包括建立网格编号矩阵：建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒。在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C。利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN。然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙。
步骤2，建立流体单元运移模型：
步骤2.1，在初始模型中设定初始流体单元状态以及边界条件。
步骤2.2，通过水压与温度以及密度关系式计算各流体单元的密度ρ i，i表示流体单元的编号。
步骤2.3，计算各流体单元的质量M i＝ρ iV i，V i为孔隙的面积。
步骤2.4，计算模型中的孔喉渗透系数：孔喉渗透系数与通道长度以及宽度有关，孔喉渗透系 数：
Figure PCTCN2020097584-appb-000001
其中，k(d)表示孔喉渗透系数，k为设定的数值，L表示两颗粒圆心之间的距离，R 1、R 2表示两颗粒的半径，minR表示所有颗粒的最小半径，l表示孔喉通道长度。
步骤2.5，计算单个流体单元的渗流量：在离散元计算的一个时间步dt内，某个流体单元pore i与周围流体单元pore n发生流体运移，流体运移服从达西定律，定义流体流入为正，单个孔喉由于压力梯度导致体积渗流量：
q i＝K(d)*A*J*dT
J＝(P n-P i)/min(R a,R b)
其中，P n和P i分别是流体单元pore n和pore i的压力，R a和R b分别是组成该孔喉的两颗粒a和b的半径，q i表示渗流量，A表示孔喉通道宽度，J表示压力梯度系数，dT表示时间步。
步骤2.6，计算总运移量：pore i在一个时间步内与周围流体单元进行的运移总量为：
Figure PCTCN2020097584-appb-000002
其中，各流体单元的密度
Figure PCTCN2020097584-appb-000003
V i为孔隙的体积,M i为孔隙内各流体单元的质量。
步骤2.7，计算运移后流体单元质量：pore i在第一个时间步后的总质量：
M′ i＝M i+Q n
其中，M′ i表示流体单元发生流体运移后的质量，M i表示运移前的质量，Q n表示总运移量。
步骤2.8，更新流体单元状态：一个时间步dt后流体单元的密度
Figure PCTCN2020097584-appb-000004
流体单元的压力通过函数关系式：P′ i＝f(ρ′ i，T i)进行更新，其中，P′ i表示，M′ i表示。
步骤3，流体固体耦合作用：孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用。
步骤3.1，固体颗粒重新平衡：由于流体的运移，各流体单元内的水压发生改变，从而各流体单元对固体颗粒之间的力也发生改变，颗粒受流体单元作用，在一个时间步内的位移距离为：
Figure PCTCN2020097584-appb-000005
其中，x i为在一个时间步内的位移距离，P i为流体单元的水压，l i为流体单元与中心颗粒的接触长度，F i为颗粒间的接触作用力，m 0为固体颗粒的质量。
步骤3.2，孔隙体积更新：由于固体颗粒的位置改变，所以流体单元的体积发生了改变，V′ i为固体颗粒改变之后孔隙的体积。
步骤3.3，流体单元状态更新：计算孔隙体积变化后的流体单元密度，在一个时间步后其密度
Figure PCTCN2020097584-appb-000006
V′ 1为固体颗粒改变之后孔隙的体积。压力为P″ i＝f(ρ″ i，T 1)，其中，ρ″ i表示流体单元的密度，P″ i表示流体单元的压力。
步骤3.4，重复步骤1到步骤3.3，即可实现二维流固耦合数值模拟。
优选的：步骤1中建立孔隙网络模型的方法：
步骤1.1，构建初始的颗粒离散元堆积模型，颗粒总数为m，颗粒由1开始连续编号到m。
步骤1.2，在离散元计算中首先获得邻居矩阵N，邻居矩阵N记录了各编号颗粒周围与其相邻的颗粒编号,由于每个颗粒单元的邻居颗粒单元数不同，行长度不一致的地方则用总颗粒数加1来填充。
步骤1.3，输入最大孔径Dmax计算获得实际连接颗粒，基于邻居矩阵N，若相邻颗粒与中心颗粒的距离小于Dmax则记为1，判断为连接。反正则为0，判断为不连接，该连接则为两颗粒之间的孔喉通道。
步骤1.4，将颗粒作为颗粒单元，在连接矩阵中根据每个颗粒单元与其邻近颗粒单元连接线相对于x轴正方向的角度Na建立连接矩阵C，连接矩阵C按照各连接的角度大小依次排列。
步骤1.5，找到反连接的位置索引，将连接矩阵C按行方向重新排列，使得第一列的数字小于第二列的数字。
步骤1.6，为了连接矩阵C依次按顺序排列，将连接矩阵C第一列的数乘以一个极大的正数加上其第二列的数，把得到的数字之和从小到大将每一行向下排列，使得成对连接相互紧邻，即反连接在正连接的下一个，并通过索引变化找出连接矩阵C中每一个连接对应的反连接的索引，建立反连接索引矩阵CI。
步骤1.7，找到下一个连接的位置索引：获得孔隙的某一连接的下一个连接，只需得到其反连接角度减小的连接的索引，即孔隙的某一连接的反连接的索引数减去1，该索引对应的连接就是下一个连接，建立对应于连接矩阵C的下一个连接矩阵CN，这样就成功的找出下一个连接的位置索引。
步骤1.8，搜索孔隙：从初始连接出发，利用下一个连接索引CN向前查找颗粒单元连接，回到直到回到初始颗粒单元，获得了一个孔隙。
步骤1.9，从编号为1的颗粒单元开始，各个连接以本连接为初始连接，向前查找颗粒单元连接，直至回到初始颗粒单元，这样就成功描述了堆积体内的所有孔隙。
步骤1.10，删除重复的记录，得到全部孔隙。
步骤1.11，根据搜索到的连接与孔隙，建立孔隙网络模型。
优选的：步骤2.2中水压与温度以及密度关系式：P i＝f(ρ i，T i)，其中，P i 表示流体单元的压力，ρ i表示各流体单元的密度，T i表示流体单元的温度。
一种基于GPU矩阵的离散元流固耦合数值模拟快速计算系统，包括输入模块、孔隙网络模型模块、流体单元运移模型模块、流固耦合作用模块，其中，
所述输入模块用于输入颗粒相关数据，在初始模型中设定初始流体单元状态以及边界条件。
所述孔隙网络模型模块用于根据颗粒相关数据建立网格编号矩阵；同时，用于建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒；在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C；利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN；然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙。
所述流体单元运移模型模块用于通过水压与温度以及密度关系式计算各流体单元的密度；计算各流体单元的质量；计算模型中的孔喉渗透系数；计算单个流体单元的渗流量；计算总运移量；计算运移后流体单元质量以及更新流体单元状态。
所述流固耦合作用模块用于计算颗粒受流体单元作用，在一个时间步内的位移距离；孔隙体积更新；计算孔隙体积变化后的流体单元密度，在一个时间步后其密度。
本发明相比现有技术，具有以下有益效果：
本发明实现所有的离散元流固耦合计算均基于GPU以纯矩阵计算形式完成，并巧妙的利用了反连接与下一个连接，成功提出了离散元孔隙搜索的方法，成功解决了在流固耦合的过程中由于孔隙在不断变化而造成的搜索孔隙难的问题，有利于进一步开展对于流固耦合的离散元模拟研究。并且提出了以孔隙流体单元来描述孔隙流体介质，为模拟孔隙介质的流固耦合提供了一种快速简洁的方法。该发明可实现大规模流固耦合系统的动态模拟，显著提高了离散元流固耦合数值模拟的计算速度以及计算数量。
附图说明
图1为建立孔隙网络模型流程图。
图2为建立流体单元运移模型。
图3为流固耦合过程示意图。
图4为颗粒堆积模型。
图5为孔隙网络模型。
图6为单个孔隙连接示意图。
图7为孔喉模型示意图。
图8为流体运移示意图。
图9为中心颗粒受力示意图。
具体实施方式
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
一种基于GPU矩阵的离散元流固耦合数值模拟方法，包括以下步骤：
步骤1，建立孔隙网络模型，包括建立网格编号矩阵，如图1所示：建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒。在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C。利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN。然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙。
步骤10，构建初始的颗粒离散元堆积模型，颗粒总数为m，颗粒由1开始连续编号到m。
步骤11，在离散元计算中首先获得邻居矩阵N，邻居矩阵N记录了各编号颗粒周围与其相邻的颗粒编号,由于每个颗粒单元的邻居颗粒单元数不同，行长度不一致的地方则用总颗粒数加1来填充。
步骤12，输入最大孔径Dmax计算获得实际连接颗粒，基于邻居矩阵N，若相邻颗粒与中心颗粒的距离小于Dmax则记为1，判断为连接。反正则为0，判断为不连接，该连接则为两颗粒之间的孔喉通道。
步骤13，将颗粒作为颗粒单元，在连接矩阵中根据每个颗粒单元与其邻近颗粒单元连接线相对于x轴正方向的角度Na建立连接矩阵C，连接矩阵C按照各连接的角度大小依次排列，如表1所示：
cList  
1 2
1 3
1 6
1 5
1 15
1 16
2 11
2 4
2 3
2 1
2 16
2 17
2 10
3 4
表1 连接矩阵
Na记录每个单元与其邻近单元连接线相对于x轴正方向的角度(逆时针)，至此已经对连接进行好了排序，接下来就是查找孔隙的步骤。
步骤14，找到反连接的位置索引，将连接矩阵C按行方向重新排列，使得第一列的数字小于第二列的数字，如表2所示：
cList2  
1 2
1 3
1 6
1 5
1 15
1 16
2 11
2 4
2 3
1 2
2 16
2 17
2 10
3 4
表2 连接矩阵按行方向重新排列
步骤15，为了连接矩阵C依次按顺序排列，将连接矩阵C第一列的数乘以一个极大的正数(本 方法的数字采用1e7)加上其第二列的数，把得到的数字之和从小到大将每一行向下排列，就可以使得成对连接相互紧邻，即反连接在正连接的下一个，并通过索引变化找出连接矩阵C中每一个连接对应的反连接的索引，建立反连接索引矩阵CI。
步骤16，找到下一个连接的位置索引：为了描述某一孔隙，必须获得形成该孔隙的所有连接。所以需要得知任意一个连接对应于某一孔隙的下一个连接，直至回到该连接，这样就成功的描述了一个孔隙。
如图6所示，1-2连接的下一个连接为2-3连接，2-3连接为1-2连接的反连接2-1在C矩阵里的角度减小的上一个连接。
所以为了获得孔隙的某一连接的下一个连接，只需得到其反连接角度减小的连接的索引，即孔隙的某一连接的反连接的索引数减去1，该索引对应的连接就是下一个连接，建立对应于连接矩阵C的下一个连接矩阵CN，这样就成功的找出下一个连接的位置索引。
步骤17，搜索孔隙：从初始连接出发，利用下一个连接索引CN向前查找颗粒单元连接，回到直到回到初始颗粒单元，获得了一个孔隙。
步骤18，从编号为1的颗粒单元开始，各个连接以本连接为初始连接，向前查找颗粒单元连接，直至回到初始颗粒单元，这样就成功描述了堆积体内的所有孔隙。
步骤19，删除重复的记录，得到全部孔隙。
步骤110，根据搜索到的连接与孔隙，建立孔隙网络模型。
步骤2，建立流体单元运移模型，孔隙内的流体单元通过颗粒单元之间的孔喉相连接，其渗透服从达西定律，其流量由岩土体性质、孔喉直径等确定，水压与温度以及密度有着一个函数关系式：P i＝f(ρ i，T i)，通过该函数的引入，实现了密度差异—压力差异驱动孔隙流体运动，如图2所示：
步骤20，初始模型需设定初始流体单元状态(温度，压力)以及边界条件。
步骤21，经过试验，水压与温度以及密度有着一个函数关系式：：P i＝f(ρ i，T i)，其中，P i表示流体单元的压力，ρ i表示各流体单元的密度，T i表示流体单元的温度。通过水压与温度以及密度关系式计算各流体单元的密度ρ i。
步骤22，计算各流体单元的质量M i＝ρ iV i，V i为孔隙的面积。
步骤23，计算模型中的孔喉渗透系数：孔喉渗透系数与通道长度以及宽度有关，如图7所示，孔喉渗透系数：
Figure PCTCN2020097584-appb-000007
其中，k(d)表示孔喉渗透系数，k为设定的数值，L表示两颗粒圆心之间的距离，R 1、R 2表示两颗粒的半径，minR表示所有颗粒的最小半径，l表示孔喉通道长度。
步骤24，计算单个流体单元的渗流量：在离散元计算的一个时间步dt内，某个流体单元pore i与周围流体单元pore n(P n的数量为P i相邻的流体单元数目)发生流体运移，如图8所示，流体运移服从达西定律，定义流体流入为正，单个孔喉由于压力梯度导致体积渗流量：
q i＝K(d)*A*J*dT
J＝(P n-P i)/min(R a,R b)
其中，P n和P i分别是流体单元pore n和pore i的压力，R a和R b分别是组成该孔喉的两颗粒a和b的半径，q i表示渗流量，A表示孔喉通道宽度，J表示压力梯度系数，dT表示时间步。
步骤25，计算总运移量：pore i在一个时间步内与周围流体单元进行的运移总量为：
Figure PCTCN2020097584-appb-000008
其中，各流体单元的密度
Figure PCTCN2020097584-appb-000009
V i为孔隙的体积,M i为孔隙内各流体单元的质量。
步骤26，计算运移后流体单元质量：pore i在第一个时间步后的总质量：
M′ i＝M i+Q n
其中，密度
Figure PCTCN2020097584-appb-000010
V i为孔隙的体积,M i为孔隙内各流体单元的质量，M′ i表示流体单元发生流体运移后的质量，M i表示运移前的质量，Q n表示总运移量。
步骤27，更新流体单元状态：一个时间步dt后流体单元的密度
Figure PCTCN2020097584-appb-000011
流体单元的压力通过函数关系式：P′ i＝f(ρ′ i，T i)进行更新，其中，P′ i表示流体单元的压力，M′ i表示流体单元的质量。
步骤3，流体固体耦合作用，如图3所示：孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用，通过这一步骤实现了流体单元与固体单元之间的交互从而实现流体与固体的耦合。
步骤30，固体颗粒重新平衡：由于流体的运移，各流体单元内的水压发生改变，从而各流体单元对固体颗粒之间的力也发生改变，颗粒受流体单元作用，如图9所示，在一个时间步内的位移距离为：
Figure PCTCN2020097584-appb-000012
其中，x i为在一个时间步内的位移距离，P i为流体单元的水压，l i为流体单元与中心颗粒的接触长度，F i为颗粒间的接触作用力，m 0为固体颗粒的质量。
步骤31，孔隙体积更新：由于固体颗粒的位置改变，所以流体单元的体积发生了改变，V′ i为固体颗粒改变之后孔隙的体积。
步骤32，流体单元状态更新：计算孔隙体积变化后的流体单元密度，在一个时间步后其密度
Figure PCTCN2020097584-appb-000013
V′ 1为固体颗粒改变之后孔隙的体积。压力为P″ i＝f(ρ″ i，T 1)，其中，ρ″ i表示流体单元更新后的密度，P″ i表示流体单元更新后的压力。
步骤33，重复步骤0到步骤32，即可实现二维流固耦合数值模拟。
根据上述一种基于GPU矩阵的离散元流固耦合数值模拟方法构建的系统，包括输入模块、孔隙网络模型模块、流体单元运移模型模块、流固耦合作用模块，其中，
所述输入模块用于输入颗粒相关数据，在初始模型中设定初始流体单元状态以及边界条件。
所述孔隙网络模型模块用于根据颗粒相关数据建立网格编号矩阵；同时，用于建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒；在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C；利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN；然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙。
所述流体单元运移模型模块用于通过水压与温度以及密度关系式计算各流体单元的密度；计算各流体单元的质量；计算模型中的孔喉渗透系数；计算单个流体单元的渗流量；计算总运移量；计算运移后流体单元质量以及更新流体单元状态。
所述流固耦合作用模块用于计算颗粒受流体单元作用，在一个时间步内的位移距离；孔隙体积更新；计算孔隙体积变化后的流体单元密度，在一个时间步后其密度。
基于matlab矩阵计算，通过把各个颗粒简化为点位，将孔隙系统简化为各个颗粒构成的网络(如图4、5所示)从而建立孔隙网络模型，从颗粒单元1出发，1-2连接、2-16连接和16-1连接形成一个三角形。根据邻居矩阵N和单元坐标，获取每个单元与邻居单元间连接角度的矩阵Na；将连接的单元从1开始，角度为连接与水平线的逆时针夹角，从0度开始排序，生成排序的连接数组C。连接具有方向，在C中先确定每个连接的反连接位置索引，如得到1-2连接的反连接2-1连接的位置，然后再确定反连接2-1连接逆时针方向下一个连接，即2-16连接。由2-16连接再往前查找，又得到16-1连接，最后沿顺时针方向得到了三角形1-2-16。在计算时，采用CI2记录C连接(如1-2)的反向连接索引(如2-1)，用CN记录C中连接的下一个连接(2-16)。从设定的初始连接开始，利用CN不断向前推算，获得三角形，四边形。最后将重复的记录去除，得到所有的孔隙，建立孔隙网络模型；孔隙中的流体通过孔喉来运移，其流量由岩土体性质、孔喉直径等确定；根据孔隙的体积应变和体积模量来确定孔隙水压力；孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用。
本发明显著提高了二维流固耦合数值模拟的运算速度，为模拟孔隙介质的流固耦合提供了一种快速简洁的方法，可实现大规模流固耦合系统的动态模拟。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

  1. 一种基于GPU矩阵的离散元流固耦合数值模拟方法，其特征在于，包括以下步骤：
    步骤1，建立孔隙网络模型，包括建立网格编号矩阵：建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒；在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C；利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN；然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙；
    步骤2，建立流体单元运移模型：
    步骤2.1，在初始模型中设定初始流体单元状态以及边界条件；
    步骤2.2，通过水压与温度以及密度关系式计算各流体单元的密度ρ i，i表示流体单元的编号；
    步骤2.3，计算各流体单元的质量M i＝ρ iV i，V i为孔隙的面积；
    步骤2.4，计算模型中的孔喉渗透系数：孔喉渗透系数与通道长度以及宽度有关，孔喉渗透系数：
    Figure PCTCN2020097584-appb-100001
    其中，k(d)表示孔喉渗透系数，k为设定的数值，L表示两颗粒圆心之间的距离，R 1、R 2表示两颗粒的半径，minR表示所有颗粒的最小半径，l表示孔喉通道长度；
    步骤2.5，计算单个流体单元的渗流量：在离散元计算的一个时间步dt内，某个流体单元pore i与周围流体单元pore n发生流体运移，流体运移服从达西定律，定义流体流入为正，单个孔喉由于压力梯度导致体积渗流量：
    q i＝K(d)*A*J*dT
    J＝(P n-P i)/min(R a,R b)
    其中，P n和P i分别是流体单元pore n和pore i的压力，R a和R b分别是组成该孔喉的两颗粒a和b的半径，q i表示渗流量，A表示孔喉通道宽度，J表示压力梯度系数，dT表示时间步；
    步骤2.6，计算总运移量：pore i在一个时间步内与周围流体单元进行的运移总量为：
    Figure PCTCN2020097584-appb-100002
    其中，各流体单元的密度
    Figure PCTCN2020097584-appb-100003
    V i为孔隙的体积,M i为孔隙内各流体单元的质量；
    步骤2.7，计算运移后流体单元质量：pore i在第一个时间步后的总质量：
    M′ i＝M i+Q n
    其中，M′ i表示流体单元发生流体运移后的质量，M i表示运移前的质量，Q n表示总运移量；
    步骤2.8，更新流体单元状态：一个时间步dt后流体单元的密度
    Figure PCTCN2020097584-appb-100004
    流体单元的压力通 过函数关系式：P′ i＝f(ρ′ i，T i)进行更新，其中，ρ′ i表示流体单元的密度，M′ i表示流体单元的质量；
    步骤3，流体固体耦合作用：孔隙水压力作用于颗粒上，而颗粒运动又改变孔隙的体积，从而实现流固耦合作用；
    步骤3.1，固体颗粒重新平衡：由于流体的运移，各流体单元内的水压发生改变，从而各流体单元对固体颗粒之间的力也发生改变，颗粒受流体单元作用，在一个时间步内的位移距离为：
    Figure PCTCN2020097584-appb-100005
    其中，x i为在一个时间步内的位移距离，P i为流体单元的压力，l i为流体单元与中心颗粒的接触长度，F i为颗粒间的接触作用力，m 0为固体颗粒的质量；
    步骤3.2，孔隙体积更新：由于固体颗粒的位置改变，所以流体单元的体积发生了改变，V′ i为固体颗粒改变之后孔隙的体积；
    步骤3.3，流体单元状态更新：计算孔隙体积变化后的流体单元密度，在一个时间步后其密度
    Figure PCTCN2020097584-appb-100006
    V′ 1为固体颗粒改变之后孔隙的体积；压力为P″ i＝f(ρ″ i，T 1)，其中，ρ″ i表示流体单元的密度，P″ i表示流体单元的压力；
    步骤3.4，重复步骤1到步骤3.3，即可实现二维流固耦合数值模拟。
  2. 根据权利要求1所述基于GPU矩阵的离散元流固耦合数值模拟方法，其特征在于：步骤1中建立孔隙网络模型的方法：
    步骤1.1，构建初始的颗粒离散元堆积模型，颗粒总数为m，颗粒由1开始连续编号到m；
    步骤1.2，在离散元计算中首先获得邻居矩阵N，邻居矩阵N记录了各编号颗粒周围与其相邻的颗粒编号,由于每个颗粒单元的邻居颗粒单元数不同，行长度不一致的地方则用总颗粒数加1来填充；
    步骤1.3，输入最大孔径Dmax计算获得实际连接颗粒，基于邻居矩阵N，若相邻颗粒与中心颗粒的距离小于Dmax则记为1，判断为连接；反正则为0，判断为不连接，该连接则为两颗粒之间的孔喉通道；
    步骤1.4，将颗粒作为颗粒单元，在连接矩阵中根据每个颗粒单元与其邻近颗粒单元连接线相对于x轴正方向的角度Na建立连接矩阵C，连接矩阵C按照各连接的角度大小依次排列；
    步骤1.5，找到反连接的位置索引，将连接矩阵C按行方向重新排列，使得第一列的数字小于第二列的数字；
    步骤1.6，为了连接矩阵C依次按顺序排列，将连接矩阵C第一列的数乘以一个极大的正数加上其第二列的数，把得到的数字之和从小到大将每一行向下排列，使得成对连接相互紧邻，即反连接在正连接的下一个，并通过索引变化找出连接矩阵C中每一个连接对应的反连接的索引，建 立反连接索引矩阵CI；
    步骤1.7，找到下一个连接的位置索引：获得孔隙的某一连接的下一个连接，只需得到其反连接角度减小的连接的索引，即孔隙的某一连接的反连接的索引数减去1，该索引对应的连接就是下一个连接，建立对应于连接矩阵C的下一个连接矩阵CN，这样就成功的找出下一个连接的位置索引；
    步骤1.8，搜索孔隙：从初始连接出发，利用下一个连接索引CN向前查找颗粒单元连接，回到直到回到初始颗粒单元，获得了一个孔隙；
    步骤1.9，从编号1的颗粒单元开始，各个连接以本连接为初始连接，向前查找颗粒单元连接，直至回到初始颗粒单元，这样就成功描述了堆积体内的所有孔隙；
    步骤1.10，删除重复的记录，得到全部孔隙；
    步骤1.11，根据搜索到的连接与孔隙，建立孔隙网络模型。
  3. 根据权利要求2所述基于GPU矩阵的离散元流固耦合数值模拟方法，其特征在于：步骤2.2中水压与温度以及密度关系式：P i＝f(ρ i，T i)，其中，P i表示流体单元的压力，ρ i表示各流体单元的密度，T i表示流体单元的温度。
  4. 一种基于权利要求1所述基于GPU矩阵的离散元流固耦合数值模拟方法的系统，其特征在于：包括输入模块、孔隙网络模型模块、流体单元运移模型模块、流固耦合作用模块，其中，
    所述输入模块用于输入颗粒相关数据，在初始模型中设定初始流体单元状态以及边界条件；
    所述孔隙网络模型模块用于根据颗粒相关数据建立网格编号矩阵；同时，用于建立初步邻近颗粒数组，每个颗粒相对于其邻近颗粒为中心颗粒；在离散元计算中首先获得邻居矩阵N，将邻居矩阵N在行方向上按连接角度进行排序，并转化为连接矩阵C；利用连接矩阵C控制搜索算法：找出连接矩阵C的反连接索引矩阵CI以及下一个连接的索引矩阵CN；然后由初始连接出发，利用下一个连接的索引矩阵CN向前查找节点，直到回到原始点，获得堆积体的全部孔隙；
    所述流体单元运移模型模块用于通过水压与温度以及密度关系式计算各流体单元的密度；计算各流体单元的质量；计算模型中的孔喉渗透系数；计算单个流体单元的渗流量；计算总运移量；计算运移后流体单元质量以及更新流体单元状态；
    所述流固耦合作用模块用于计算颗粒受流体单元作用，在一个时间步内的位移距离；孔隙体积更新；计算孔隙体积变化后的流体单元密度，在一个时间步后其密度。
PCT/CN2020/097584 2020-06-05 2020-06-23 基于gpu矩阵的离散元流固耦合数值模拟方法及系统 Ceased WO2021243765A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202010503940.8A CN111507024B (zh) 2020-06-05 2020-06-05 基于gpu矩阵的离散元流固耦合数值模拟方法及系统
CN202010503940.8 2020-06-05

Publications (1)

Publication Number Publication Date
WO2021243765A1 true WO2021243765A1 (zh) 2021-12-09

Family

ID=71878795

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/097584 Ceased WO2021243765A1 (zh) 2020-06-05 2020-06-23 基于gpu矩阵的离散元流固耦合数值模拟方法及系统

Country Status (2)

Country Link
CN (1) CN111507024B (zh)
WO (1) WO2021243765A1 (zh)

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114239367A (zh) * 2021-12-31 2022-03-25 西南石油大学 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
CN114386302A (zh) * 2021-12-31 2022-04-22 西南石油大学 一种非定常流固耦合多相渗流模型构建方法
CN114638178A (zh) * 2022-03-25 2022-06-17 中国海洋石油集团有限公司 一种基于离散元的疏松砂岩微粒运移模拟方法
CN114912333A (zh) * 2022-07-19 2022-08-16 上海索辰信息科技股份有限公司 一种轻流体下和结构的声振耦合的模拟方法
CN115422818A (zh) * 2022-11-03 2022-12-02 北京云庐科技有限公司 基于云服务的离散元并行实时仿真边坡预警系统及方法
CN115726382A (zh) * 2022-11-28 2023-03-03 中国地质大学（武汉） 一种确定基坑落底式止水帷幕接缝渗漏的方法及装置
CN116611359A (zh) * 2023-05-10 2023-08-18 中国地质大学（武汉） Co2注入过程中矿物颗粒动态运移和堵塞数值模拟方法
CN120874206A (zh) * 2025-09-28 2025-10-31 中铁二十局集团第二工程有限公司 一种反复渗流作用下浸水路基渗透破坏判别方法

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112417715B (zh) * 2020-10-26 2022-11-04 山东大学 一种真三轴伺服加载状态下的岩体破裂模拟方法及系统
CN112818611B (zh) * 2021-01-28 2023-08-22 南京大学 一种单裂隙岩石水力压裂过程流固耦合的数值模拟方法
CN113221431B (zh) * 2021-05-14 2022-05-06 湖北理工学院 基于颗粒离散元与格子Boltzmann的压缩渗透试验数值模拟方法
CN113405966B (zh) * 2021-06-08 2022-08-23 浙江广天构件集团股份有限公司 一种水泥基材料颗粒堆积体系孔径分布计算方法
CN113656656B (zh) * 2021-07-21 2024-01-26 南京南力科技有限公司 用于宽级配离散元颗粒系统的高效邻居检索方法及系统
CN119397870B (zh) * 2024-10-14 2025-09-23 清华大学 基于cfd-dem耦合的电池喷发颗粒物沉积传热计算方法
CN119312462B (zh) * 2024-12-17 2025-04-04 山东科技大学 一种基于固液衔接的岩体注浆模拟方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180371874A1 (en) * 2016-03-11 2018-12-27 Halliburton Energy Services, Inc. Fracture network fluid flow simulation with enhanced fluid-solid interaction force determination
CN109977505A (zh) * 2019-03-13 2019-07-05 南京大学 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
CN110263362A (zh) * 2019-04-25 2019-09-20 南京大学（苏州）高新技术研究院 基于孔隙密度流的岩土体离散元流固耦合数值模拟方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006072566A (ja) * 2004-08-31 2006-03-16 Mizuho Information & Research Institute Inc 流体構造連成解析方法及び流体構造連成解析プログラム
CN105507893A (zh) * 2015-12-07 2016-04-20 西南石油大学 预测页岩储层体积改造后生产效果的方法
CN107462508B (zh) * 2017-08-16 2018-10-02 西南石油大学 一种多场耦合渗流多功能实验装置及测试方法
CN108729908B (zh) * 2018-05-21 2021-11-09 中国石油大学(华东) 一种基于孔隙网络模型的致密油流动模拟及渗透率预测方法
CN108830020B (zh) * 2018-07-12 2019-05-31 西南石油大学 一种模拟海上油田微压裂增注裂缝扩展的方法
CN109684685B (zh) * 2018-12-10 2019-12-17 青岛海洋地质研究所 一种多分支井水合物降压开采条件下的产能及储层稳定性分析方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20180371874A1 (en) * 2016-03-11 2018-12-27 Halliburton Energy Services, Inc. Fracture network fluid flow simulation with enhanced fluid-solid interaction force determination
CN109977505A (zh) * 2019-03-13 2019-07-05 南京大学 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
CN110263362A (zh) * 2019-04-25 2019-09-20 南京大学（苏州）高新技术研究院 基于孔隙密度流的岩土体离散元流固耦合数值模拟方法

Cited By (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2023099406A (ja) * 2021-12-31 2023-07-13 西南石油大学 室内岩石コア実験における、デジタル化された多相流体－構造連成（ｆｓｉ）浸透数値シミュレーション手法
CN114386302A (zh) * 2021-12-31 2022-04-22 西南石油大学 一种非定常流固耦合多相渗流模型构建方法
CN114239367B (zh) * 2021-12-31 2022-12-30 西南石油大学 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
CN114386302B (zh) * 2021-12-31 2023-02-10 西南石油大学 一种非定常流固耦合多相渗流模型构建方法
CN114239367A (zh) * 2021-12-31 2022-03-25 西南石油大学 一种室内岩心的数字化多相流固耦合渗流数值模拟方法
CN114638178A (zh) * 2022-03-25 2022-06-17 中国海洋石油集团有限公司 一种基于离散元的疏松砂岩微粒运移模拟方法
CN114912333A (zh) * 2022-07-19 2022-08-16 上海索辰信息科技股份有限公司 一种轻流体下和结构的声振耦合的模拟方法
CN115422818A (zh) * 2022-11-03 2022-12-02 北京云庐科技有限公司 基于云服务的离散元并行实时仿真边坡预警系统及方法
CN115422818B (zh) * 2022-11-03 2023-02-10 北京云庐科技有限公司 基于云服务的离散元并行实时仿真边坡预警系统及方法
CN115726382A (zh) * 2022-11-28 2023-03-03 中国地质大学（武汉） 一种确定基坑落底式止水帷幕接缝渗漏的方法及装置
CN116611359A (zh) * 2023-05-10 2023-08-18 中国地质大学（武汉） Co2注入过程中矿物颗粒动态运移和堵塞数值模拟方法
CN116611359B (zh) * 2023-05-10 2025-08-05 中国地质大学（武汉） Co2注入过程中矿物颗粒动态运移和堵塞数值模拟方法
CN120874206A (zh) * 2025-09-28 2025-10-31 中铁二十局集团第二工程有限公司 一种反复渗流作用下浸水路基渗透破坏判别方法

Also Published As

Publication number Publication date
CN111507024B (zh) 2023-10-13
CN111507024A (zh) 2020-08-07

Similar Documents

Publication Publication Date Title
WO2021243765A1 (zh) 基于gpu矩阵的离散元流固耦合数值模拟方法及系统
CN114218840B (zh) 河口航道水沙运动及其地形演变整体建模及可视化系统
CN110263362B (zh) 基于孔隙密度流的岩土体离散元流固耦合数值模拟方法
CN109977505B (zh) 基于gpu矩阵计算的离散元孔隙系统快速搜索方法
CN105224704B (zh) 一种海底不稳定性的评价方法
CN115688590B (zh) 一种基于深度学习的多矿物相岩心渗透率预测方法
CN106991244B (zh) 一种基于图论的裂隙网络连通性及渗流计算的方法
CN117077577B (zh) 一种适用于低渗透裂缝性油藏的快速模拟及优化方法
CN108446502B (zh) 一种利用完整二维浅水方程组获得流域单位线的方法
CN111046567B (zh) 一种基于Godunov格式的城市排水管网水流数值模拟方法
CN104991999A (zh) 一种基于二维sph的溃坝洪水演进模拟方法
CN108197345B (zh) 一种基于耦合理论的长线输水结构动力时程分析建模方法
CN112818611B (zh) 一种单裂隙岩石水力压裂过程流固耦合的数值模拟方法
CN117763895A (zh) 一种混凝土重力坝动力响应敏感性分析方法及系统
CN110147561A (zh) 一种含天然裂缝致密油气储层体积压裂缝网预测方法
CN107145671B (zh) 一种油藏数值模拟方法及系统
CN109918821A (zh) 一种迎风守恒型河道漫溢出流数值模拟方法
CN114154430A (zh) 一种压裂油藏co2驱油流动模拟方法
CN105404758A (zh) 一种基于有限单元法的固体连续介质变形的数值模拟方法
CN103149600A (zh) 一种基于优化控制点确定水文地质参数的自动配线方法
CN111625954A (zh) 一种降雨型滑坡模型trigrs的并行优化方法及系统
Lippmann et al. Analysis of the response of geothermal reservoirs under injection and production procedures
CN109783991B (zh) 一种无已知底滑面的滑坡滑动过程模拟方法
CN118350319B (zh) 河口最大浑浊带变动趋势预测方法、设备、介质及产品
CN116796579A (zh) 二氧化碳驱压裂数模一体化方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20939002

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20939002

Country of ref document: EP

Kind code of ref document: A1