CN106021828B - A kind of reservoir fluid analogy method based on Lattice Boltzmann model - Google Patents

A kind of reservoir fluid analogy method based on Lattice Boltzmann model Download PDF

Info

Publication number
CN106021828B
CN106021828B CN201610554629.XA CN201610554629A CN106021828B CN 106021828 B CN106021828 B CN 106021828B CN 201610554629 A CN201610554629 A CN 201610554629A CN 106021828 B CN106021828 B CN 106021828B
Authority
CN
China
Prior art keywords
fluid
particle
velocity
function
simulation
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.)
Expired - Fee Related
Application number
CN201610554629.XA
Other languages
Chinese (zh)
Other versions
CN106021828A (en
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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to CN201610554629.XA priority Critical patent/CN106021828B/en
Publication of CN106021828A publication Critical patent/CN106021828A/en
Application granted granted Critical
Publication of CN106021828B publication Critical patent/CN106021828B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/30—Circuit design
    • G06F30/36—Circuit design at the analogue level
    • G06F30/367—Design verification, e.g. using simulation, simulation program with integrated circuit emphasis [SPICE], direct methods or relaxation methods

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Microelectronics & Electronic Packaging (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention discloses a kind of fluid simulation method based on Lattice Boltzmann model.The analogy method includes, by the image gridding of porous media, uses fi(x, y, t) represents mesh point I (x, y) place, and movement velocity is ciParticle corresponding to particle distribution；Judge Particles Moving direction ciIt is to make particle be distributed f whether towards Gu Bi bordersi(x, y, t) performs reverse function, otherwise to fi(x, y, t) performs the collision function of Lattice Boltzmann model, is then distributed f according to the particle after evolutioni(x, y, t) obtains fluid density ρ ' (x, y) and fluid velocity u ' (x, y)；Until meeting evolution termination condition.The present invention is by by the image gridding of porous media, and by the discrete particle for turning to motion of the fluid in grid, so as to be distributed f according to particlei(x, y, t) obtains fluid density and speed, improves the operational efficiency of simulation.

Description

A kind of reservoir fluid analogy method based on grid-Boltzmann model
Technical field
The invention belongs to hydrodynamics and computer high-performance computing sector, grid-glass is based on more particularly, to one kind The fluid simulation method of the graceful model of Wurz.
Background technology
Rock and soil are all natural porous media, and its pore diameter is about between 1 μm~500 μm.Porous media The problem of modelling of middle flow of fluid is widely present in the fields such as oil exploitation, environmental protection, biochip and Chemical Engineering, example Such as, during oil exploitation, pressure influence the permeability and recovery ratio of oil reservoir, studies the shadow of pressure-driven fluid flow Sound can play directive function to recovery method.Therefore, simulated using numerical method and computer technology to solving above-mentioned neck The practical problem in domain has highly important meaning.
In actual applications, because the complicated of porous media, BORDER PROCESSING are difficult, at the same user to it is extensive, The requirement of efficient and Fast simulation is improved constantly, and this make it that flow of fluid faces significant challenge in simulation porous media.Patent text Offer the prediction that the A of CN 103776739 disclose a kind of free-boundary problem of Robertson-Si Difu fluids in porous media Method.But this method needs to be related to the minimum pore radius of porous media, maximum pore radius and capillary in the calculation The multiple parameters such as straight length so that calculation procedure is relatively complicated；Secondly, this method is only capable of mould starting pressure, it is impossible to pre- The speed and pressure versus flow body of fluid measured flow follow-up influence, and predictive ability is limited.
The content of the invention
For the disadvantages described above or Improvement requirement of prior art, grid-Boltzmann mould is based on the invention provides one kind The fluid simulation method of type, its object is to by by the image gridding of porous media, and by the fluid discretization in grid For the particle of motion, so as to be distributed f according to particlei(x, y, t) obtains fluid density and speed, improves the operational efficiency of simulation.
It is an aspect of this invention to provide that there is provided a kind of analogy method of fluid, comprising the following steps：
S1. the image of porous media is carried out obtaining simulation lattice I (x, y) after gridding processing, and obtains the simulation Fluid viscosity ν, fluid density ρ (x, y, t), fluid velocity u (x, y, t) and the Gu Bi of the corresponding nondimensionalizations of grid I (x, y) Border；Wherein, x, y represent the abscissa and ordinate of mesh point on simulation lattice I (x, y) respectively, and x is 1~Nx, y is 1~ Ny；
S2. moment t=0 is made, divergent density function, the discrete velocity letter of fluid are set up using grid-Boltzmann model Number and equilibrium distribution function；The divergent density function isIt is described Discrete velocity function isThe equilibrium state point Cloth function isfi(x,y,t) Mesh point (x, y) place is represented, movement velocity is ciParticle corresponding to particle distribution, i represent particle rapidity direction numbering；Represent particle distribution fiThe equilibrium state of (x, y, t), ωiFor weight coefficient；
S3. particle rapidity c is judgediThe direction of motion whether towards Gu Bi borders, be then by grid-Boltzmann model Reverse function, make speed ciCorresponding particle is distributed fi(x, y, t) is in moment t+ δtValue be to be distributed with the reverse particles of i, Otherwise the collision function of grid-Boltzmann model is performedWherein, when Between step-lengthδxFor the length of side of mesh point, csFor the velocity of sound, θiThe angle for the positive direction rotate counterclockwise for being direction i along x-axis, Slack time
S4. by the divergent density function and discrete velocity function after evolution, moment t+ δ is obtainedtWhen fluid density ρ ' (x, y) and fluid velocity u ' (x, y)；Divergent density function after the evolution isDiscrete velocity function after evolution is
S5. judge whether to meet evolution termination condition, be, simulation terminates, and otherwise makes moment t=t+ δt, return to step S3。
Preferably, in the step S1, the fluid viscosity of nondimensionalization For the kinematic viscosity system of fluid Number, L0For the characteristic length of fluid, u0For the characteristic velocity of fluid, and u0=L0/ 1 second；
After the step S5, in addition to, obtain the actual speed of fluid
As it is further preferred that the characteristic length of the fluidOr For simulation lattice I (x, y) Overall width,For simulation lattice I (x, y) total height.
Preferably, in the step S1, the method to the image progress gridding processing of porous media is specially：To be many The image of hole medium is converted into the gray level image of gridding, gray level image described in binaryzation, and removes noise, obtains the simulation Grid I (x, y).
As it is further preferred that the method for the binaryzation is Two-peak method, Otsu methods, Bernsen methods or cycle threshold Method.
As it is further preferred that the method for removing noise is connection domain method.
Preferably, in the step S2, during i=0, ci=(0,0)；During i=1~8, θiRespectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, correspondence particle distribution fiParticle in (x, y, t) is in time step δtInterior motion apart from ci δtFor the distance of another mesh point in from mesh point I (x, y) to its eight neighborhood；
In the step S3, the reverse function is：fi(x, y, t)=fi+2(x, y, t), i=1,2,5 or 6；fi(x, Y, t)=fI-2(x, y, t), i=3,4,7 or 8.
As it is further preferred that in the step S2, weight coefficient ωiMeet：
Preferably, the evolution termination condition in the step S5 is：Moment t whether >=incremental time T, or fluid Whether speed u ' (x, y) meetsErr is speed convergence condition.
As it is further preferred that the incremental time T in the step S5 is 10000 δt~1000000 δt。
As it is further preferred that the speed convergence condition err in the step S5 is 0~0.1.
Preferably, after the step S5, in addition to by fluid velocity u ' (x, y) visualize.
The invention has the advantages that：
1st, the present invention is by by the image gridding of porous media, using the region where porous media as Gu Bi borders, Realize the setting for the complicated calculations operating mode in porous media；
2nd, it is of the invention by the discrete particle for turning to motion of the fluid in grid so that interparticle collision and particle are with consolidating The effect on wall border can parallel processing, so as to improve the computational efficiency and calculating speed of simulation；
3rd, because the calculating of the particle of different mesh points in the present invention can be run simultaneously, the present invention can utilize OpenACC etc. Method carries out parallel optimization, suitable for being performed in a variety of accelerator platform parallelizations.
Brief description of the drawings
Fig. 1 constitutes the direction of motion schematic diagram of the particle of fluid for the present invention；
Fig. 2 for present invention composition fluid Particles Moving towards Gu Bi borders and reverse schematic diagram；
The schematic diagram that Fig. 3 moves for the particle of present invention composition fluid to neighbor mesh points；
Fig. 4 is the flow chart of the embodiment of the present invention 1；
Fig. 5 is the porous media image schematic diagram after binaryzation in the embodiment of the present invention 1；
Pressure boundaries and Gu Bi border schematic diagram of the Fig. 6 for the embodiment of the present invention 1；
Fig. 7 is that the fluid velocity of the embodiment of the present invention 1 visualizes schematic diagram；
Fig. 8 is the flow chart for carrying out hardware optimization to the embodiment of the present invention 1 using OpenACC technologies；
Fig. 9 is that the calculating run in different platform when the iterative steps of the embodiment of the present invention 1 are 10000 times takes statistical chart.
Embodiment
In order to make the purpose , technical scheme and advantage of the present invention be clearer, it is right below in conjunction with drawings and Examples The present invention is further elaborated.It should be appreciated that the specific embodiments described herein are merely illustrative of the present invention, and It is not used in the restriction present invention.As long as in addition, technical characteristic involved in each embodiment of invention described below Not constituting conflict each other can just be mutually combined.
The invention provides a kind of fluid simulation method, comprise the following steps；
S11. plane right-angle coordinate is set up, and the image lattice for the porous media for being by size turns to Nx × Ny gray scale Image, gray level image described in binaryzation, and the noise on gray level image is removed with methods such as connected domains, obtain the simulation lattice I (x, y), wherein, x, y represent the abscissa and ordinate of mesh point on I (x, y) respectively, and x is 1~Nx, and y is 1~Ny；It is described The method of binaryzation is Two-peak method, Otsu methods, Bernsen methods or cycle threshold method；
S12. the region where fluid is distinguished with the matrix of the binaryzation of the correspondence simulation lattice I (x, y) is situated between with porous Region where matter, can generally represent the region where porous media with numerical value 1, and numerical value 0 represents the region where fluid；
S13. according to simulation lattice I (x, y), the image of porous media sizeThe actual density of fluid The kinematic viscosity coefficient of fluidAnd pressure boundary, obtain the fluid viscosity ν of nondimensionalization, fluid velocity u (x, y, T), fluid density ρ (x, y, t) and Gu Bi borders；
Wherein, the kinematic viscosity coefficient of fluid and the property and temperature of fluid are relevant, and the kinematics of usual crude oil is glued Property coefficient ν is in 2m2/ s~300m2Between/s, for example, kinematic viscosity coefficient of No. 50 oil at 20 DEG C is about 50m2/s；
The fluid viscosity of nondimensionalizationL0For the characteristic length of fluid, u0For the characteristic velocity of fluid, generally orderOru0=L0/ 1 second, fluid velocity was usually 0, even u (x, y, t)=(0,0)；
The porous Jie of fluid density ρ (x, y, t) and Gu Bi borders then in heel pressure border and simulation lattice I (x, y) Two pairs of borders in simulation lattice in actual applications, are generally respectively set to pressure boundary and Gu Bi borders by qualitative correlation；Example Such as, during oil exploitation, the fluid in petroleum pipeline is applied when upwarding pressure P ', then region where porous media, Left margin and right margin are Gu Bi borders；
Fluid in petroleum pipeline is under the situation of not stressing, and its actual pressure is usually an atmospheric pressure P0Plus fluid certainly The pressure gh, g of body represent acceleration of gravity, and h represents the depth of fluid, and the pressure under stress is P0+gh+P’； During oil exploitation, it is assumed that the pressure suffered by a pressure boundary is pout, the pressure suffered by another pressure boundary is pin, The fluid density that assume that wherein one pressure boundary is 1；Again due to the fluid pressure p (x, y, t) and fluid of nondimensionalization Density p (x, y, t) meets p (x, y, t)=ρ (x, y, t) × cs 2, csThe velocity of sound is represented, then the density of another pressure boundary is pin/ pout；For example, during oil exploitation, the fluid in petroleum pipeline is applied when upwarding pressure P ', ρ (x, 1, t) ≡ 1, ρ (x, Ny,t)≡pin/pout；
S2. moment t=0 is made, divergent density function, the discrete velocity letter of fluid are set up using grid-Boltzmann model Number and equilibrium distribution function；The divergent density function isIt is described Discrete velocity function isThe equilibrium state point Cloth function isfi(x,y,t) Mesh point I (x, y) place is represented, movement velocity is ciParticle corresponding to particle distribution, i represent particle rapidity direction compile Number；Represent particle distribution fiThe equilibrium state of (x, y, t), ωiFor weight coefficient；
For example, can be with 0~8 direction for representing 8 kinds of Particles Movings, as shown in Figure 1；As i=0, ci=(0,0)；Work as i When=1~8, particle distribution fiThe direction c of (x, y, t) corresponding Particles MovingiRespectively along the positive direction rotate counterclockwise of x-axis Angle, θi, it is respectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, is present in particle distribution fiIn (x, y, t) Particle is in time step δtInterior motion apart from ciδtFor another mesh point in from mesh point I (x, y) to its eight neighborhood away from From；Time stepδxFor the length of side of mesh point, csFor the velocity of sound, ωiFor weight coefficient；Therefore discrete velocity is metWeight coefficient ωiMeet：
S3. particle rapidity c is judgediThe direction of motion whether towards Gu Bi borders, be to make corresponding particle distribution fi(x, Y, t) in moment t+ δtValue is distributed for reverse particle, otherwise performs the collision function of grid-Boltzmann modelSlack time
Understood for i in above-mentioned steps S2 for 0~8 technical scheme, if now particle rapidity ciDirection of motion i directions During Gu Bi borders, as shown in Figure 2；
fi(x, y, t)=fi+2(x, y, t) i=1,2,5,6；fi(x, y, t)=fI-2(x, y, t) i=3,4,7,8；
Otherwise, it is present in particle distribution fiParticle in (x, y, t) can be in time step δtOne mesh point of interior motion away from From as shown in figure 3, the mesh point f of correspondence direction 1,2,3,5,6iThe side of (x, y, t) from center O respectively shown in A, B, C, D, E To moving to adjacent mesh point；
S4. according to the divergent density function and discrete velocity function after evolution
Obtain moment t+ δtWhen fluid density ρ ' (x, y) and fluid velocity u ' (x, y)；
S5. whether moment t is judged >=incremental time T, or whether fluid velocity u ' (x, y) meetsErr is speed convergence condition, is, simulation terminates, and otherwise makes moment t=t+ δt, return to step S3；
Usual incremental time T is 10000 δt~1000000 δt, the speed convergence condition err is 0~0.1, and it specifically takes It is worth size and the computational accuracy correlation with analog image I (x, y), the size of analog image is bigger, and computational accuracy is higher, then institute The time that need to be calculated is longer, therefore incremental time T is bigger, and speed convergence condition err is smaller.
S6. fluid velocity u ' (x, y) can finally be visualized, or by obtaining the actual speed of fluid To make more intuitively observation or further processing.
Embodiment 1
The analogy method of fluid of the present embodiment 1 based on grid-Boltzmann model specifically includes following steps, such as Fig. 4 It is shown：
S1. gather porous media image, by the image of the porous media of the size of (=10m × 5m) save as Nx × The picture file of Ny pixels；Wherein, Nx and Ny are respectively the number of pixels (Nx=in the present embodiment of horizontal direction and vertical direction 400, Ny=200)；Picture file is determined that the threshold value of binaryzation is 127 using Two-peak method, and carries out noise reduction process, acquisition Analog image I (x, y) is as shown in figure 5, wherein black region represents the solid portion in porous media, and white portion represents fluid Region；With the matrix of binaryzation to should analog image；The element of matrix only includes 0 and 1, the calculating for representing porous media Operating mode, wherein, 1 represents the region where porous media, and numerical value 0 represents the region where fluid；
S2. according to practical problem, zoning, primary condition, boundary condition, physical parameter etc. are determined, and according to matrix, Mesh generation is carried out to zoning；
The size of zoning is generally identical with size Nx × Ny of the picture file in step S1, specific to need according to many The size of the image of hole medium, the parallel processing speeds of computer and the computational accuracy of requirement are set；
Boundary condition in the present embodiment includes pressure boundary and Gu Bi borders, except the region where porous media with Outside, it is Gu Bi borders up-and-down boundary also to be selected in the present embodiment, and left margin is entrance boundary, corresponding pressure pin=36.67, Right margin is outlet border, corresponding pressure pout=33.33；
Physical parameter in the present embodiment includes the actual density of fluidThe kinematics of fluid Viscosity learns coefficientThe actual speed of fluidCharacteristic length Characteristic density ρ0 =1 × 103kg/m3, characteristic velocity u0=L0/ 1 second=10m/s；
According to formulaObtain this implementation The primary condition of example, including：Fluid viscosity ν, fluid velocity u (x, y, t)=0, the fluid density ρ (x, y, t)=1 of nondimensionalization (x=2~(Nx-1)) and boundary condition ρ (1, y, t) ≡ 1, ρ (1, Nx, t) ≡ 1.1；csFor the velocity of sound in fluid, it is necessary to full The following condition of foot, umax/cs＜ ＜ 1, wherein umaxIt is the maximum for the fluid velocity estimated, c is taken in the present embodiments=5.77；
S2. moment t=0 is made,X=2 ~(Nx-1)；Wherein, fi(x, y, t) represents mesh point I (x, y) place, and movement velocity is ciParticle corresponding to particle distribution, I represents the direction numbering of particle rapidity；As i=0, ci=(0,0)；As i=1~8, particle rapidity ciDirection be respectively Along the positive direction rotate counterclockwise 0 of x-axis, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, particle is in time step δtIt is interior Motion apart from ciδtFor the distance of another mesh point in from mesh point I (x, y) to its eight neighborhood；Time stepδx For the length of side of mesh point, (δ in the present embodimentx=2.5 × 10- 3, δt=2.5 × 10- 4)；, can also be according to net in actually calculating Lattice are divided, storage allocation, maroscopic quantity and distribution function for storing each mesh point, and according to physical parameter, to calculating area Maroscopic quantity and equilibrium distribution function in domain are initialized；
S3. equilibrium distribution function is made
Wherein, ωiFor Weight coefficient；And according to the element value of the corresponding matrix of mesh point, judge particle rapidity ciThe direction of motion whether towards Gu Bi sides Boundary；It is to make corresponding particle distribution fi(x, y, t) value is distributed for reverse particle, otherwise performs collision functionSlack timeWeight coefficient ωiMeet：ω0=4/9, ω5=ω6=ω7=ω8=1/9, ω1=ω1=ω1=ωi= 1/36；, can be on each mesh point in actual calculating process, the collision function in one time step of execution, and Gu Bi borders and pressure boundary at boundary net point, exercise boundary processing；
S4. according to equation group
Obtain moment t+ δtWhen, fluid density ρ ' (x, y) and fluid velocity u ' (x, y) on each mesh point；
S5. whether moment t is judged >=incremental time T, or whether fluid velocity u ' (x, y) meetsErr is calculation error, is, simulation terminates, and otherwise makes t=t+ δt, return Step S3；In the present embodiment, err=0.000001, incremental time T=M_max × δ are takent, M_max is the iterations upper limit, M_max=1000000 in the present embodiment, is that simulation terminates, otherwise return to step S3；In the present embodiment, it have passed through Iteration evolution (the i.e. t=630000 δ of 630000 stepst=157, the corresponding real time be 157s) after reached expected calculating Error；
S6. the fluid velocity u ' (x, y) of output is visualized using line of vector integration method using the poster processing soft, As shown in fig. 7, wherein the black lines with arrow represent the flow field of flow of fluid, according to formula Also the actual speed of fluid can be obtained
The present embodiment can be compiled as LB (serial grid-Boltzmann algorithm) code, and using OpenACC technologies in CPU+ In GPU platform, parallel optimization is carried out to serial LB codes, analyze can be parallel in serial LB codes code segment, and according to analysis As a result, it is that code segment that can be parallel adds OpenACC and line identifier, analyzes the data transmission bottle neck on CPU+GPU platforms, and root According to analysis result, addition OpenACC data transfers are identified, as shown in Figure 8.
For example, matrix can be read at CPU ends, carry out mesh generation and perform the initialization of flow field parameter and model parameter Journey, by the bump flow process of particle distribution function, flow field maroscopic quantity statistic processes, BORDER PROCESSING in Lattice Boltzmann Method Process, error control procedure etc. are computationally intensive and the intensive part called is transformed using OpenACC Parallelizing Techniques, make it It can be performed in a variety of accelerator platform parallelizations, when different accelerator platform operations, only need to use same set of generation Code is it is achieved that so as to considerably reduce multi-platform development and maintenance cost.
Fluid density ρ ' (x, y) and fluid speed in the bump flow process of distribution function in step s3, step S4 U ' (x, y) statistic processes is spent, the difference calculation process runs in step S5 can be using OpenACC technologies in the parallelization of accelerator end Perform, each calculation procedure can distribute a GPU thread and carry out operation independent, so as to improve the arithmetic speed of simulation.
Meanwhile, for the visualization problem in step S6, in order to solve the data transfer between CPU internal memories and GPU video memorys Problem, the data transfer to CPU+GPU platforms is optimized：In the beginning of code, stated using OpenACC and apply showing Space is deposited, is performed only in output data file once from GPU video memorys to the data transmission procedure of CPU internal memories.It can so show The data transmission period between reduction CPU internal memories and GPU video memorys is write, computational efficiency is improved.Meanwhile, this optimization does not interfere with generation Parallelization operation of the code on multi thread CPU platform.
The arithmetic speed test result of the present embodiment is as follows：
The test platform and test parameter of the present embodiment are as shown in table 1.Test platform includes multi thread CPU and accelerates platform Accelerate platform with CPU+GPU, operation program includes fluid in the serial Lattice Boltzmann Method porous media based on C++ Fluid-flow analogy in flow simulating program and the multi-platform Lattice Boltzmann Method porous media based on OpenACC of support Program.The parameter that two programs are used when calculating in test process is identical.
The test platform of table 1 and test parameter
In fig .9, give the calculating that the program of the present embodiment runs on CPU to take, while also giving utilization The Lattice Boltzmann program of OpenACC parallel optimizations accelerates platform and CPU+GPU to accelerate what is run on platform in multi thread CPU Calculate time-consuming.The calculating provided in figure is time-consuming be 10 tests average result, to ensure the stability of test result.Can be with Find out, accelerate the calculating of the OpenACC concurrent programs of platform to take as 241.54 seconds using four thread CPU, calculating, which takes, is The 48.78% of serial program, and accelerate the calculating of the OpenACC concurrent programs of platform to take as 22.95 seconds using CPU+GPU, It is only the 4.63% of serial program to calculate time-consuming.
It can be seen that present invention utilizes grid-Boltzmann model according to the implementation steps and test result of the present invention Fluid is simulated, so as to improve operational efficiency.This method can be using OpenACC technologies to bump flow mistake therein Journey, BORDER PROCESSING process, statistics maroscopic quantity process and calculation error control condition process carry out parallelization transformation, so as to be applied to Multi thread CPU accelerates platform and CPU+GPU to accelerate platform, further increases the computational efficiency of fluid flow problem, while A program allow while operating on different acceleration platforms, so as to significantly reduce the exploitation of multi-platform concurrent program And maintenance cost.
As it will be easily appreciated by one skilled in the art that the foregoing is merely illustrative of the preferred embodiments of the present invention, it is not used to The limitation present invention, any modifications, equivalent substitutions and improvements made within the spirit and principles of the invention etc., it all should include Within protection scope of the present invention.

Claims (8)

1. a kind of reservoir fluid analogy method, it is characterised in that comprise the following steps：
S1. the image of porous media is carried out obtaining simulation lattice I (x, y) after gridding processing, and obtains the simulation lattice I The fluid viscosity ν, fluid density ρ (x, y, t), fluid velocity u (x, y, t) and Gu Bi borders of (x, y) corresponding nondimensionalization； Wherein, x, y represent the abscissa and ordinate of mesh point on simulation lattice I (x, y) respectively, and x is 1~Nx, and y is 1~Ny；
S2. make moment t=0, using grid-Boltzmann model set up the divergent density function of fluid, discrete velocity function with And equilibrium distribution function；The divergent density function isIt is described discrete Velocity function isThe equilibrium distribution letter Number isfi(x, y, t) is represented Mesh point (x, y) place, movement velocity is ciParticle corresponding to particle distribution, i represent particle rapidity direction numbering；Represent particle distribution fiThe equilibrium state of (x, y, t), ωiFor weight coefficient；
S3. particle rapidity c is judgediThe direction of motion whether towards Gu Bi borders, be then by the anti-of grid-Boltzmann model To function, speed c is madeiCorresponding particle is distributed fi(x, y, t) is in moment t+ δtValue be to be distributed with the reverse particles of i, otherwise Perform the collision function of grid-Boltzmann modelWherein, when Between step-lengthδxFor the length of side of mesh point, csFor the velocity of sound, θiThe folder for the positive direction rotate counterclockwise for being direction i along x-axis Angle, slack time
S4. by the divergent density function and discrete velocity function after evolution, moment t+ δ is obtainedtWhen fluid density ρ ' (x, Y) with fluid velocity u ' (x, y)；Divergent density function after the evolution isDiscrete velocity function after evolution is
S5. judge whether to meet evolution termination condition, be, simulation terminates, and otherwise makes moment t=t+ δt, return to step S3；It is described Evolution termination condition is：Moment t whether >=incremental time T, or whether fluid velocity u ' (x, y) meetErr is speed convergence condition.
2. fluid simulation method as claimed in claim 1, it is characterised in that in the step S1, the fluid of nondimensionalization Viscosity For the kinematic viscosity coefficient of fluid, L0For the characteristic length of fluid, u0For the characteristic velocity of fluid, and u0 =L0/ 1 second；
After the step S5, in addition to, obtain the actual speed of fluid
3. fluid simulation method as claimed in claim 2, it is characterised in that the characteristic length of the fluidOr For simulation lattice I (x, y) overall width,For simulation lattice I (x, y) total height.
4. fluid simulation method as claimed in claim 1, it is characterised in that in the step S1, to the figure of porous media As the method for carrying out gridding processing is specially：The image of porous media is converted into the gray level image of gridding, binaryzation institute Gray level image is stated, and removes noise, the simulation lattice I (x, y) is obtained.
5. fluid simulation method as claimed in claim 1, it is characterised in that in the step S2, during i=0, ci=(0, 0)；During i=1~8, θiRespectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, correspondence particle distribution fi(x,y,t) In particle in time step δtInterior motion apart from ciδtFor another mesh point in from mesh point I (x, y) to its eight neighborhood Distance；
In the step S3, the reverse function is：fi(x, y, t)=fi+2(x, y, t), i=1,2,5 or 6；fi(x,y,t) =fI-2(x, y, t), i=3,4,7 or 8.
6. fluid simulation method as claimed in claim 5, it is characterised in that in the step S2, weight coefficient ωiMeet： ω0=4/9, ω1=ω2=ω3=ω4=1/9, ω5=ω6=ω7=ω8=1/36.
7. fluid simulation method as claimed in claim 1, it is characterised in that the incremental time T in the step S5 is 10000δt~1000000 δt。
8. fluid simulation method as claimed in claim 1, it is characterised in that the speed convergence condition in the step S5 Err is 0~0.1.
CN201610554629.XA 2016-07-15 2016-07-15 A kind of reservoir fluid analogy method based on Lattice Boltzmann model Expired - Fee Related CN106021828B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610554629.XA CN106021828B (en) 2016-07-15 2016-07-15 A kind of reservoir fluid analogy method based on Lattice Boltzmann model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610554629.XA CN106021828B (en) 2016-07-15 2016-07-15 A kind of reservoir fluid analogy method based on Lattice Boltzmann model

Publications (2)

Publication Number Publication Date
CN106021828A CN106021828A (en) 2016-10-12
CN106021828B true CN106021828B (en) 2017-09-08

Family

ID=57118931

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610554629.XA Expired - Fee Related CN106021828B (en) 2016-07-15 2016-07-15 A kind of reservoir fluid analogy method based on Lattice Boltzmann model

Country Status (1)

Country Link
CN (1) CN106021828B (en)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106709191A (en) * 2016-12-29 2017-05-24 中国石油大学(北京) Numerical simulation method and apparatus for seismic wave field
CN108267390B (en) * 2016-12-30 2020-07-10 中国石油天然气股份有限公司 Method for determining gas permeability of reservoir containing nanopores
CN108455733A (en) * 2018-01-22 2018-08-28 太原理工大学 A kind of biological film model construction method of film biological sewage processing
CN108345964B (en) * 2018-02-10 2021-09-21 北京师范大学 Water quality prediction method and system based on water quality model
CN111368487B (en) * 2020-03-17 2023-07-18 广西师范大学 Flow field processing method for simulating periodic movement of particles based on lattice Boltzmann model
CN111858066B (en) * 2020-07-30 2022-07-15 中国空气动力研究与发展中心超高速空气动力研究所 CPU + GPU heterogeneous parallel optimization method in pneumatic theory unified algorithm
CN112666059A (en) * 2020-12-14 2021-04-16 中国石油大学(华东) Method for determining gas-water relative permeability of porous medium in gas hydrate decomposition process
CN112613243B (en) * 2020-12-16 2023-10-20 中国科学院深圳先进技术研究院 Method, device and computer readable storage medium for hydrodynamic simulation
CN113268901B (en) * 2021-04-12 2023-12-22 东南大学 Lattice Boltzmann dynamic pressure gas bearing gap micro-flow simulation method
CN113125325B (en) * 2021-04-26 2021-11-16 东北石油大学 Coal rock fracture characteristic characterization and permeability simulation method
CN115795989B (en) * 2023-01-28 2023-06-16 安世亚太科技股份有限公司 Fluid motion simulation method, simulation terminal, electronic equipment and medium
CN115952753B (en) * 2023-03-14 2023-05-12 中国测绘科学研究院 Chip flow simulation method, system and equipment combining CA and LBM

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102945295B (en) * 2012-10-15 2015-09-02 浪潮(北京)电子信息产业有限公司 A kind of parallel acceleration method of Lattice Boltzmann Method and system
CN105278346B (en) * 2015-11-06 2018-06-12 北京航空航天大学 A kind of hot fluid emulation mode based on the double distributed models of discrete LATTICE BOLTZMANN

Also Published As

Publication number Publication date
CN106021828A (en) 2016-10-12

Similar Documents

Publication Publication Date Title
CN106021828B (en) A kind of reservoir fluid analogy method based on Lattice Boltzmann model
JP5629704B2 (en) Physics simulation on graphics processor
CN104268943A (en) Fluid simulation method based on Eulerian-Lagrangian coupling method
CN101727653B (en) Graphics processing unit based discrete simulation computation method of multicomponent system
CN107633123A (en) A kind of method accelerated for smoothed particle method simulation bleeding and processing
Fietz et al. Optimized hybrid parallel lattice Boltzmann fluid flow simulations on complex geometries
CN111695309B (en) High-performance large-scale fluid-solid coupling fluid simulation method based on statistical dynamics
CN1415090A (en) Method and program for simulating physical system using object-oriented programming
CN105354879A (en) Particle spring structure based universal garment three-dimensional model simulation method and system
CN102831275B (en) A kind of emulation mode of 3D fluid and system
WO2017150626A1 (en) Particle simulation device, particle simulation method, and particle simulation program
CN106687827A (en) Method for producing precipitation-strengthening-type martensitic stainless steel
CN110717269A (en) Fluid surface detail protection method based on grid and particle coupling
CN105760588A (en) SPH fluid surface reconstruction method based on second-layer regular grid
CN103426196B (en) A kind of animation modeling method under fluid environment
Ma et al. A parallel meshless dynamic cloud method on graphic processing units for unsteady compressible flows past moving boundaries
CN1932884A (en) Process type ground fast drawing method based on fractal hierarchical tree
CN104318601A (en) Human motion simulation method under fluid environment
CN108427605B (en) Acceleration method for realizing streamline simulation based on particle tracking algorithm
Washington et al. Micro‐mechanical simulation of geotechnical problems using massively parallel computers
Hung et al. Automatic clustering method for real-time construction simulation
CN112100939B (en) Real-time fluid simulation method and system based on computer loader
Gaudin et al. Optimising Hydrodynamics applications for the Cray XC30 with the application tool suite
Norman Targeting atmospheric simulation algorithms for large, distributed-memory, GPU-accelerated computers
Dai et al. Interactive smoke simulation and rendering on the GPU

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20170908

Termination date: 20190715