CN106021828A - Fluid simulation method based on grid-boltzmann model - Google Patents

Abstract

The invention discloses a fluid simulation method based on a grid-boltzmann model. The simulation method comprises the following steps: gridding an image of a porous medium, expressing a grid point (x, y) with fi (x, y, t) and taking a kinematic velocity as the particle distribution corresponding to ci particle; judging if the moving direction ci of the particle is oriented to the boundary of fixed wall, ordering the particle distribution fi (x, y, t) to execute reverse function or collision function of the grid-boltzmann model, and then acquiring the fluid density p' (x, y) and fluid speed u' (x, y) according to the particle distribution fi (x, y, t) after evolution; finishing till meeting the evolution ending condition. According to the fluid simulation method provided by the invention, the image of the porous medium is gridded and the fluid in the grid is dispersed as the moving particles, so that the fluid density and speed can be acquired according to the particle distribution fi (x, y, t) and the simulation running efficiency can be increased.

Description

A kind of fluid simulation method based on grid-Boltzmann model

Technical field

The invention belongs to hydrodynamics and computer high-performance computing sector, more particularly, to one based on grid-glass 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 fluid flowing is widely present in the fields such as oil exploitation, environmental conservation, biochip and Chemical Engineering, example As, during oil exploitation, pressure influence permeability and the recovery ratio of oil reservoir, the shadow of research pressure-driven fluid flow Recovery method can be played directive function by sound.Therefore, numerical method and computer technology is utilized to be simulated solving above-mentioned neck The practical problem in territory has highly important meaning.

In actual applications, owing to the structure of porous media is complicated, BORDER PROCESSING difficulty, simultaneously user to extensive, Efficient and Fast simulation requirement improves constantly, and this makes to simulate fluid flowing in porous media and faces significant challenge.Patent literary composition Offer CN 103776739A and disclose the prediction side of a kind of Robertson-Si Difu fluid free-boundary problem in porous media Method.But the method needs to relate to the minimum pore radius of porous media, maximum pore radius and capillary tube in the calculation Multiple parameter such as straight length so that calculation procedure is relatively complicated；Secondly, the method is only capable of mould starting pressure, it is impossible to prediction The speed of fluid and the follow-up impact of flowing of pressure versus flow body, it was predicted that limited in one's ability.

Summary of the invention

For disadvantages described above or the Improvement requirement of prior art, the invention provides a kind of based on grid-Boltzmann mould The fluid simulation method of type, its object is to by by porous media is image gridding, and by the fluid discretization in grid For motion particle, thus according to particle be distributed f_i(x, y t) obtain fluid density and speed, improve the operational efficiency of simulation.

It is an aspect of this invention to provide that provide the analogy method of a kind of fluid, comprise the following steps:

S1. (x y), and obtains described simulation to obtain simulation lattice I after the image to porous media carries out gridding process Grid I (x, y) the corresponding fluid viscosity ν of nondimensionalization, fluid density ρ (x, y, t), fluid velocity u (x, y, t) and Gu Bi Border；Wherein, x, y represent respectively simulation lattice I (x, y) abscissa of upper mesh point and vertical coordinate, x is 1～Nx, y be 1～ Ny；

S2. make moment t=0, utilize grid-Boltzmann model to set up the divergent density function of fluid, discrete velocity letter Number and equilibrium distribution function；Described divergent density function isDescribed Discrete velocity function isDescribed equilibrium state divides Cloth function isf_i(x, y, t) table Show that (x, y) place, movement velocity is c to mesh point_iThe particle corresponding to particle distribution, i represent particle rapidity direction numbering；Represent particle distribution f_i(x, y, equilibrium state t), ω_iFor weight coefficient；

S3. particle rapidity c is judged_iThe direction of motion whether towards Gu Bi border, be then by the reverse letter of grid-Boltzmann model Number, makes speed c_iCorresponding particle distribution f_i(x, y, t) at moment t+ δ_tValue be the distribution of reverse with i particle, otherwise perform grid-glass The collision function of the graceful model of Wurz Wherein, time stepδ_xFor the length of side of mesh point, c_sFor the velocity of sound, θ_iRevolve counterclockwise along the positive direction of x-axis for direction i The angle turned, slack time

S4. by the divergent density function after developing and discrete velocity function, it is thus achieved that moment t+ δ_tTime stream Body density ρ ' (x, y) and fluid velocity u ' (x, y)；Divergent density function after described evolution isDiscrete velocity function after evolution is

S5. judging whether to meet evolution termination condition, be, simulation terminates, and otherwise makes moment t=t+ δ_t, return step S3。

Preferably, in described step S1, the fluid viscosity of nondimensionalization Kinematic viscosity system for fluid Number, L₀For the characteristic length of fluid, u₀For the characteristic velocity of fluid, and u₀=L₀/ 1 second；

After described step S5, also include, it is thus achieved that the actual speed of fluid

As it is further preferred that the characteristic length of described fluidOr For simulation lattice I (x, y) Overall width,For simulation lattice I (x, total height y).

Preferably, in described step S1, the image of porous media is carried out the method for gridding process particularly as follows: will 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, it is thus achieved that described simulation Grid I (x, y).

As it is further preferred that the method for described binaryzation is Two-peak method, Otsu method, Bernsen method or cycle threshold Method.

As it is further preferred that the method for described removal noise is connection domain method.

Preferably, in described step S2, during i=0, c_i=(0,0)；During i=1～8, θ_iBe respectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, corresponding particle distribution f_i(x, y, t) in particle in time step δ_tDistance c of interior motion_i δ_tFor from mesh point I, (x, y) to the distance of another mesh point in its eight neighborhood；

In described step S3, described reverse function is: f_i(x, y, t)=f_i+2(x, y, t), i=1,2,5 or 6；f_i(x, Y, t)=f_I-2(x, y, t), i=3,4,7 or 8.

As it is further preferred that in described step S2, weight coefficient ω_iMeet:

$\left\{\begin{array}{cc}{\mathrm{\ω}}_i=4/9& i=0,\\ {\mathrm{\ω}}_i=1/9& i=1~4,\\ {\mathrm{\ω}}_i=1/36& i=5~8.\end{array}\right.$

Preferably, the described evolution termination condition in described step S5 is: moment t whether >=incremental time T, or fluid Speed u ' (whether x y) meetsErr is speed convergence condition.

As it is further preferred that the described incremental time T in described step S5 is 10000 δ_t～1000000 δ_t。

As it is further preferred that described speed convergence condition err in described step S5 is 0～0.1.

Preferably, after described step S5, also include fluid velocity u ' (x, y) visualization.

The method have the advantages that

1, image gridding by by porous media of the present invention, using the region at porous media place as Gu Bi border, Achieve the setting for the complicated calculations operating mode in porous media；

2, the present invention is by discrete for the fluid in the grid particle turning to motion so that interparticle collision and particle are with solid The effect on wall border can parallel processing, thus improve simulation computational efficiency and calculate speed；

3, can run due to the calculating of the particle of mesh points different in the present invention, the present invention may utilize OpenACC etc. simultaneously Method carries out parallel optimization, is suitable to perform in multiple accelerator platform parallelization.

Accompanying drawing explanation

Fig. 1 is the direction of motion schematic diagram that the present invention forms the particle of fluid；

Fig. 2 is that the present invention forms the Particles Moving of fluid towards Gu Bi border reverse schematic diagram；

Fig. 3 is the schematic diagram that the present invention forms that the particle of fluid moves to neighbor mesh points；

Fig. 4 is the embodiment of the present invention 1 flow chart；

Fig. 5 is the porous media image schematic diagram in the embodiment of the present invention 1 after binaryzation；

Fig. 6 is pressure boundary and the Gu Bi border schematic diagram of the embodiment of the present invention 1；

Fig. 7 is that the embodiment of the present invention 1 fluid velocity visualizes schematic diagram；

Fig. 8 is to utilize OpenACC technology that the embodiment of the present invention 1 is carried out the flow chart of hardware optimization；

Fig. 9 is that in different platform, the operation embodiment of the present invention 1 iterative steps is the time-consuming cartogram of calculating when 10000 times.

Detailed description of the invention

In order to make the purpose of the present invention, technical scheme and advantage clearer, below in conjunction with drawings and Examples, right The present invention is further elaborated.Should be appreciated that specific embodiment described herein only in order to explain the present invention, and It is not used in the restriction present invention.If additionally, technical characteristic involved in each embodiment of invention described below The conflict of not constituting each other just can be mutually combined.

The invention provides a kind of fluid simulation method, comprise the following steps；

S11. set up plane right-angle coordinate, and the image lattice of porous media a size of is turned to the gray scale of Nx × Ny Image, gray level image described in binaryzation, and remove the noise on gray level image by methods such as connected domains, it is thus achieved that described simulation lattice (x, y), wherein, x, y represent that (x, y) goes up abscissa and the vertical coordinate of mesh point to I, and x is 1～Nx, and y is 1～Ny to I respectively；Described The method of binaryzation is Two-peak method, Otsu method, Bernsen method or cycle threshold method；

S12. with corresponding described simulation lattice I, (x, the matrix of binaryzation y) is distinguished the region at fluid place and is situated between with porous The region at matter place, generally can represent the region at porous media place with numerical value 1, and numerical value 0 represents the region at fluid place；

S13. according to simulation lattice I (x, y), the size of the image of porous mediaThe actual density of fluid The kinematic viscosity coefficient of fluidAnd pressure boundary, it is thus achieved that the fluid viscosity ν of nondimensionalization, fluid velocity u (x, y, T), fluid density ρ (x, y, t) and Gu Bi border；

Wherein, the kinematic viscosity coefficient of fluid is relevant with the character of fluid and temperature, and the kinesiology of usual crude oil glues Property coefficient ν is at 2m²/ s～300m²Between/s, such as, No. 50 oil kinematic viscosity coefficient when 20 DEG C is about 50m²/s；

The fluid viscosity of nondimensionalizationL₀For the characteristic length of fluid, u₀For the characteristic velocity of fluid, generally makeOru₀=L₀/ 1 second, fluid velocity was usually 0, even u (x, y, t)=(0,0)；

Fluid density ρ (x, y, t) and then heel pressure border, Gu Bi border and simulation lattice I (x, y) in porous be situated between Two pairs of borders in simulation lattice in actual applications, are generally respectively set to pressure boundary and Gu Bi border by qualitative correlation；Example As, during oil exploitation, in petroleum pipeline fluid apply upward pressure P ' time, then porous media place region, Left margin and right margin are Gu Bi border；

Fluid in petroleum pipeline under the situation of not stressing, a usually atmospheric pressure P of its actual pressure₀Add upper fluid certainly The pressure gh, g of body represents that acceleration of gravity, h represent the degree of depth of fluid, and the pressure under stress is P₀+gh+P’；? During oil exploitation, it is assumed that the pressure suffered by a pressure boundary is p_out, the pressure suffered by another pressure boundary is p_in, The fluid density that assume that a wherein pressure boundary is 1；Again due to nondimensionalization fluid pressure p (x, y, t) and fluid (x, y t) meet p (x, y, t)=ρ (x, y, t) × c to density p_s ², c_sRepresent the velocity of sound, then the density of another pressure boundary is p_in/ p_out；Such as, during oil exploitation, in petroleum pipeline fluid apply upward pressure P ' time, ρ (x, 1, t) ≡ 1, ρ (x, Ny,t)≡p_in/p_out；

S2. make moment t=0, utilize grid-Boltzmann model to set up the divergent density function of fluid, discrete velocity letter Number and equilibrium distribution function；Described divergent density function isDescribed Discrete velocity function isDescribed equilibrium state divides Cloth function isf_i(x, y, t) table Show that (x, y) place, movement velocity is c to mesh point I_iThe particle corresponding to particle distribution, i represent particle rapidity direction numbering；Represent particle distribution f_i(x, y, equilibrium state t), ω_iFor weight coefficient；

For example, it is possible to represent the direction of 8 kinds of Particles Movings with 0～8, as shown in Figure 1；As i=0, c_i=(0,0)；Work as i When=1～8, particle distribution f_i(x, y, t) the direction c of corresponding Particles Moving_iIt is respectively the positive direction along x-axis to rotate counterclockwise Angle, θ_i, it is respectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, is present in particle distribution f_i(x, y, t) in Particle is in time step δ_tDistance c of interior motion_iδ_tFor from mesh point I (x, y) to another mesh point in its eight neighborhood away from From；Time stepδ_xFor the length of side of mesh point, c_sFor the velocity of sound, ω_iFor weight coefficient；Therefore discrete velocity meetsWeight coefficient ω_iMeet:

$\left\{\begin{array}{cc}{\mathrm{\ω}}_i=4/9& i=0,\\ {\mathrm{\ω}}_i=1/9& i=1~4,\\ {\mathrm{\ω}}_i=1/36& i=5~8;\end{array}\right.$

S3. particle rapidity c is judged_iThe direction of motion whether towards Gu Bi border, be, make correspondence particle distribution f_i(x, Y, t) at moment t+ δ_tValue is reverse particle distribution, otherwise performs the collision function of grid-Boltzmann model

Slack time

Understand for the technical scheme that i in above-mentioned steps S2 is 0～8, if now particle rapidity c_iDirection of motion i towards During Gu Bi border, as shown in Figure 2；

f_i(x, y, t)=f_i+2(x, y, t) i=1,2,5,6；f_i(x, y, t)=f_I-2(x, y, t) i=3,4,7,8；

Otherwise, particle distribution f it is present in_i(x, y, t) in particle can be in time step δ_tInterior one mesh point of motion away from From, as it is shown on figure 3, the mesh point f of correspondence direction 1,2,3,5,6_i(x, y, t) from center O respectively along the side shown in A, B, C, D, E To moving to adjacent mesh point；

S4. according to the divergent density function after developing and discrete velocity function

$\left\{\begin{array}{c}{\mathrm{\ρ}}^{\′}(x,y)=\underset{i}{\Σ}{f}_i(x+c_i{\mathrm{\δ}}_t{\mathrm{Cos\θ}}_i,y+c_i{\mathrm{\δ}}_t{\mathrm{Sin\θ}}_i,t+{\mathrm{\δ}}_t)\\ u^{\′}(x,y)=\frac{1}{{\mathrm{\ρ}}^{\′}(x,y)}\underset{i}{\Σ}{c}_i{f}_i(x+c_i{\mathrm{\δ}}_t{\mathrm{Cos\θ}}_i,y+c_i{\mathrm{\δ}}_t{\mathrm{Sin\θ}}_i,t+{\mathrm{\δ}}_t)\end{array}\right.$

Obtain moment t+ δ_tTime fluid density ρ ' (x, y) and fluid velocity u ' (x, y)；

S5. whether judge moment t >=incremental time T, or fluid velocity u ' (whether x y) meetsErr is speed convergence condition, is, simulation terminates, and otherwise makes moment t=t+ δ_t, return step S3；

Generally incremental time T is 10000 δ_t～1000000 δ_t, described speed convergence condition err is 0～0.1, and it specifically takes (x, size and computational accuracy y) are relevant, and the size of analog image is the biggest, and computational accuracy is the highest, then institute with analog image I for value The time that need to calculate is the longest, and therefore incremental time T is the biggest, and speed convergence condition err is the least.

The most finally can be by fluid velocity u ' (x, y) visualization, or the actual speed by acquisition fluid To make to observe more intuitively or further process.

Embodiment 1

The analogy method of the present embodiment 1 fluid based on grid-Boltzmann model specifically includes following steps, such as Fig. 4 Shown in:

S1. gather the image of porous media, the image of the porous media of the size of (=10m × 5m) is saved as Nx × The picture file of Ny pixel；Wherein, Nx and Ny is respectively the number of pixels (Nx=in the present embodiment of horizontal direction and vertical direction 400, Ny=200)；Picture file utilizes Two-peak method determine, and the threshold value of binaryzation is 127, and carries out noise reduction process, it is thus achieved that (x, y) as it is shown in figure 5, solid portion during wherein black region represents porous media, white portion represents fluid to analog image I Region；With the matrix of binaryzation to should analog image；Entry of a matrix element only includes 0 and 1, for representing the calculating of porous media Operating mode, wherein, 1 represents the region at porous media place, and numerical value 0 represents the region at fluid place；

S2. according to practical problem, zoning, initial condition, boundary condition, physical parameter etc. are determined, and according to matrix, Zoning is carried out stress and strain model；

The size of zoning is generally identical with the size Nx × Ny of the picture file in step S1, and concrete needs are according to many The size of the image of hole medium, the parallel processing speeds of computer and the computational accuracy of requirement are arranged；

Boundary condition in the present embodiment includes pressure boundary and Gu Bi border, except porous media place region with Outward, the most also selecting up-and-down boundary is Gu Bi border, and left margin is entrance boundary, corresponding pressure p_in=36.67, Right margin is outlet border, corresponding pressure p_out=33.33；

Physical parameter in the present embodiment includes the actual density of fluidThe kinesiology of fluid Viscosity coefficientThe actual speed of fluidCharacteristic length Characteristic density ρ₀ =1 × 10³kg/m³, characteristic velocity u₀=L₀/ 1 second=10m/s；

According to formulaObtain this enforcement The initial condition of example, including: the fluid viscosity ν of nondimensionalization, fluid velocity u (x, y, t)=0, fluid density ρ (x, y, t)=1 (x=2～(Nx-1)) and boundary condition ρ (1, y, t) ≡ 1, ρ (1, Nx, t) ≡ 1.1；c_sFor the velocity of sound in fluid, need full The following condition of foot, u_max/c_s＜ ＜ 1, wherein u_maxIt is the maximum of the fluid velocity estimated, takes c in the present embodiment_s=5.77；

S2. moment t=0 is made,X=2 ～(Nx-1)；Wherein, f_i(x, y t) represent that (x, y) place, movement velocity is c to mesh point I_iThe particle corresponding to particle distribution, I represents the direction numbering of particle rapidity；As i=0, c_i=(0,0)；As i=1～8, particle rapidity c_iDirection be respectively Positive direction along x-axis rotates 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4 counterclockwise, and particle is in time step δ_tIn Distance c of motion_iδ_tFor from mesh point I, (x, y) to the distance of another mesh point in its eight neighborhood；Time step δ_xFor the length of side of mesh point, (δ in the present embodiment_x=2.5 × 10^-3, δ_t=2.5 × 10^-4)；In Practical Calculation, also can be according to Stress and strain model, storage allocation, for storing maroscopic quantity and the distribution function of each mesh point, and according to physical parameter, to calculating Maroscopic quantity and equilibrium distribution function in region initialize；

S3. equilibrium distribution function is made

Wherein, ω_iFor Weight coefficient；And according to entry of a matrix element value corresponding to mesh point, it is judged that particle rapidity c_iThe direction of motion whether towards Gu Bi limit Boundary；It is to make the particle distribution f of correspondence_i(x, y, t) value is reverse particle distribution, otherwise performs collision function

Slack timeWeight coefficient ω_iMeet: ω₀=4/9, ω₅=ω₆=ω₇=ω₈=1/9, ω₁=ω₁= ω₁=ω_i=1/36；During Practical Calculation, the collision in a time step can be performed on each mesh point Function, and the Gu Bi border at boundary net point and pressure boundary, exercise boundary processes；

S4. according to equation group

$\left\{\begin{array}{c}{\mathrm{\ρ}}^{\′}(x,y)=\underset{i}{\Σ}{f}_i(x+c_i{\mathrm{\δ}}_t{\mathrm{Cos\θ}}_i,y+c_i{\mathrm{\δ}}_t{\mathrm{Sin\θ}}_i,t+{\mathrm{\δ}}_t)\\ u^{\′}(x,y)=\frac{1}{{\mathrm{\ρ}}^{\′}(x,y)}\underset{i}{\Σ}{c}_i{f}_i(x+c_i{\mathrm{\δ}}_t{\mathrm{Cos\θ}}_i,y+c_i{\mathrm{\δ}}_t{\mathrm{Sin\θ}}_i,t+{\mathrm{\δ}}_t)\end{array}\right.$

Obtain moment t+ δ_tTime, fluid density ρ on each mesh point ' (x, y) and fluid velocity u ' (x, y)；

S5. whether judge moment t >=incremental time T, or fluid velocity u ' (whether x y) meetsErr, for calculating 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 taken_t, M_max is the iterations upper limit, M_max=1000000 in the present embodiment, is that simulation terminates, and otherwise returns step S3；In the present embodiment, have passed through Iteration evolution (the i.e. t=630000 δ of 630000 steps_t=157, corresponding actual time is 157s) after reached intended calculating Error；

S6. use the poster processing soft to output fluid velocity u ' (x, y) utilizes line of vector integration method to visualize, As it is shown in fig. 7, wherein the black lines with arrow represents the flow field that fluid flows, 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 utilizes OpenACC technology at CPU+ In GPU platform, serial LB code is carried out parallel optimization, analyze can be parallel in serial LB code code segment, and according to analysis As a result, for can be parallel code segment add OpenACC line identifier, analyze the data transmission bottle neck on CPU+GPU platform, and root According to analysis result, add OpenACC data transmission mark, as shown in Figure 8.

Such as, matrix can be read at CPU end, carry out stress and strain model 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 utilizes OpenACC Parallelizing Techniques to transform so that it is Can perform in multiple accelerator platform parallelization, different accelerator platforms runs when, only need to use same set of generation Code is it is achieved that thus considerably reduce multi-platform development and maintenance cost.

The bump flow process of distribution function in step s3, fluid density ρ in step S4 ' (x, y) with fluid speed (difference calculation process runs in step S5 may utilize OpenACC technology on accelerator end parallelization ground to degree u ' for x, statistic processes y) Performing, each calculation procedure can be distributed a GPU thread and carry out operation independent, thus be improve the arithmetic speed of simulation.

Meanwhile, for the visualization problem in step S6, in order to solve the data transmission between CPU internal memory and GPU video memory Problem, the data transmission to CPU+GPU platform is optimized: in the beginning of code, utilizes OpenACC to state and applies for aobvious Deposit space, only when output data file perform once from GPU video memory to the data transmission procedure of CPU internal memory.So can show Write the data transmission period reduced between CPU internal memory and GPU video memory, improve computational efficiency.Meanwhile, this optimization does not interferes with generation Code parallelization on multi thread CPU platform runs.

The arithmetic speed test result of the present embodiment is as follows:

Test platform and the test parameter of the present embodiment are as shown in table 1.Test platform includes multi thread CPU and accelerates platform Accelerating platform with CPU+GPU, operation program includes fluid in the Lattice Boltzmann Method porous media based on C++ of serial Fluid-flow analogy in the Lattice Boltzmann Method porous media based on OpenACC that flow simulating program is multi-platform with support Program.The parameter that in test process, two programs are used when calculating is identical.

Table 1 test platform and test parameter

In fig .9, give the calculating that the program of the present embodiment runs on CPU time-consuming, give also utilization simultaneously The Lattice Boltzmann program of OpenACC parallel optimization accelerates operation on platform and CPU+GPU acceleration platform in multi thread CPU Calculate time-consuming.The calculating be given in figure is time-consumingly the average result of 10 tests, to ensure the stability of test result.Permissible Finding out, the calculating using the OpenACC concurrent program of four thread CPU acceleration platforms is time-consumingly 241.54 seconds, and calculating is time-consumingly The 48.78% of serial program, and using the calculating of the OpenACC concurrent program of CPU+GPU acceleration platform is time-consumingly 22.95 seconds, Calculate the 4.63% of time-consuming only serial program.

Enforcement step and test result according to the present invention are it can be seen that present invention utilizes grid-Boltzmann model Convection cell is simulated, thus improves operational efficiency.The method may utilize OpenACC technology to bump flow mistake therein Journey, BORDER PROCESSING process, statistics maroscopic quantity process and calculating error control conditioning process carry out parallelization transformation, thus are applicable to Multi thread CPU accelerates platform and CPU+GPU accelerates platform, further increases the computational efficiency of fluid flow problem, the most also Make a program can be simultaneously run in different acceleration on platform, thus 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 and the foregoing is only presently preferred embodiments of the present invention, not in order to Limit the present invention, all any amendment, equivalent and improvement etc. made within the spirit and principles in the present invention, all should comprise Within protection scope of the present invention.

Claims (9)

1. a fluid simulation method, it is characterised in that comprise the following steps:

S1. (x y), and obtains described simulation lattice I to obtain simulation lattice I after the image to porous media carries out gridding process (x, y) the corresponding fluid viscosity ν of nondimensionalization, fluid density ρ (x, y, t), fluid velocity u (x, y, t) and Gu Bi border； Wherein, x, y represent that (x, y) goes up abscissa and the vertical coordinate of mesh point to simulation lattice I, and x is 1～Nx, and y is 1～Ny respectively；

S2. make moment t=0, utilize grid-Boltzmann model set up the divergent density function of fluid, discrete velocity function with And equilibrium distribution function；Described divergent density function isDescribed discrete Velocity function isDescribed equilibrium distribution letter Number is(x, y t) represent fi (x, y) place, movement velocity is c to mesh point_iThe particle corresponding to particle distribution, i represent particle rapidity direction numbering；Represent particle distribution f_i(x, y, equilibrium state t), ω_iFor weight coefficient；

S3. particle rapidity c is judged_iThe direction of motion whether towards Gu Bi border, be then by the reverse function of grid-Boltzmann model, Make speed c_iCorresponding particle distribution f_i(x, y, t) at moment t+ δ_tValue be the distribution of reverse with i particle, otherwise perform grid-Bohr The collision function of the most graceful model Wherein, time stepδ_xFor the length of side of mesh point, c_sFor the velocity of sound, θ_iRotate counterclockwise along the positive direction of x-axis for direction i Angle, slack time

S4. by the divergent density function after developing and discrete velocity function, it is thus achieved that moment t+ δ_tTime fluid density ρ ' (x, y) and stream Body speed u ' (x, y)；Divergent density function after described evolution is Discrete velocity function after evolution is

S5. judging whether to meet evolution termination condition, be, simulation terminates, and otherwise makes moment t=t+ δ_t, return step S3.

2. mould fluid simulation method as claimed in claim 1, it is characterised in that in described step S1, the stream of nondimensionalization Body viscosity For the kinematic viscosity coefficient of fluid, L₀For the characteristic length of fluid, u₀For the characteristic velocity of fluid, and u₀=L₀/ 1 second；

After described step S5, also include, it is thus achieved that the actual speed of fluid

3. fluid simulation method as claimed in claim 2, it is characterised in that the characteristic length of described fluidOr For simulation lattice I (x, overall width y),For simulation lattice I (x, total height y).

4. fluid simulation method as claimed in claim 1, it is characterised in that the figure in described step S1, to porous media As carrying out the method for gridding process particularly as follows: the image of porous media to be converted into the gray level image of gridding, binaryzation institute State gray level image, and remove noise, it is thus achieved that and described simulation lattice I (x, y).

5. fluid simulation method as claimed in claim 1, it is characterised in that in described step S2, during i=0, c_i=(0, 0)；During i=1～8, θ_iIt is respectively 0, pi/2, π, 3 pi/2s, π/4,3 π/4,5 π/4 and 7 π/4, corresponding particle distribution f_i(x,y,t) In particle in time step δ_tDistance c of interior motion_iδ_tFor from mesh point I, (x, y) to another mesh point in its eight neighborhood Distance；

In described step S3, described reverse function is: f_i(x, y, t)=f_i+2(x, y, t), i=1,2,5 or 6；f_i(x,y,t) =f_I-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 described step S2, weight coefficient ω_iMeet: ω₀=4/9, ω₁=ω₂=ω₃=ω₄=1/9, ω₅=ω₆=ω₇=ω₈=1/36.

7. fluid simulation method as claimed in claim 1, it is characterised in that the described evolution in described step S5 Termination condition is: moment t whether >=incremental time T, or fluid velocity u ' (whether x y) meetsErr is speed convergence condition.

8. fluid simulation method as claimed in claim 7, it is characterised in that the described incremental time T in described step S5 is 10000δ_t～1000000 δ_t。

9. fluid simulation method as claimed in claim 7, it is characterised in that the described speed convergence condition in described step S5 Err is 0～0.1.