CN104866695A - GPU-accelerated fluid-structure coupling simulation method through immersion boundary and lattice Boltzmann methods - Google Patents

Abstract

The invention provides a GPU-accelerated fluid-structure coupling simulation method through immersion boundary and lattice Boltzmann methods. The GPU-accelerated fluid-structure coupling simulation method comprises the following steps: 1, establishing a structure body to be detected and an internal fluid immersion boundary body shape in professional pre-processing software, and dividing a boundary lattice for program read; reading boundary body shape data through a program, and setting a boundary condition of a point required to be measured and a program model parameter; calculating the boundary acting force applied to a fluid by using a finite element method according to a structure of the boundary at the current time and a reference structure; diffusing the boundary acting force to surrounding fluids through the Dirac function; solving a Navier-Stokes equation carrying an external force term by using the lattice Boltzmann method; interpolating the fluid speed into the boundary through the Dirac function to obtain a boundary movement speed to further update the boundary position; repeatedly performing the steps 3-6 till the updated boundary position reaches an calculated ending point. According to the GPU-accelerated fluid-structure coupling simulation method, structure bodies in any shapes can be processed, an underlying fluid lattice does not require to be reconstructed in the calculation process, and the calculation efficiency is high.

Description

A kind of immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method accelerated through GPU

Technical field

The invention belongs to the crossing domain of computational fluid dynamics and Solid State Mechanics of Computation, particularly high-performance calculation is in the application in above-mentioned field.

Background technology

Interaction between fluid and deformable structure, at nature ubiquity, is seen everywhere especially especially in biological tissue and organ, and typical example comprises insect wing, fin, human heart valve and vocal cords.Although said structure belongs to and separates plane and learn or the research object in physiology field, it is most important that they complete the function of specifying in three-dimensional kinematics character and the spatial form that changes at any time to biosome.Because the inherent complicacy in the drastic deformation that structure experiences and surrounding flow field, so correctly simulate one of fluid structurecoupling problem a major challenge remaining academia under three-dimensional condition up to now.

Traditional analog stream that is used for consolidates the method for coupled problem usually based on fit grid, and typical example comprises Arbitrary Lagrangian-Euler method and finite element method.The inferior position of said method is that bottoms stream grid needs to adjust according to locations of structures, and in order to prevent the drastic deformation of grid from needing to carry out special processing.Generally speaking mess generation is process of a high price, particularly needing calculating consuming time in a large number when planform is complicated to upgrade mesh topology, making to push classic method to practical application and having difficulties.

Therefore, be necessary very much for a kind of new fluid structurecoupling analogy method of efficiency exploitation, can not with under reducing prerequisite that computational accuracy is cost, shortening computing time, for Practical Project provides reliable guiding opinion.

Immersed boundary method has obtained extensively and successfully application as the numerical algorithm of the solid coupled problem of a kind of special analog stream since proposition.In this method, fluid and structure adopt Euler method and Lagrangian method to describe respectively, and be used for the orthogonal cartesian grid of discrete fluid and structure and unstructured grid respectively and there is no need to coincide together, and this just immersed boundary method compared to a large advantage of classic method, solve the difficult point that each time step all needs to re-construct according to structure position and form grid.The interaction of fluid and structure is realized by Dirac function, utilize this function the acting force be defined on structure grid node can be diffused to peripheral fluid net point, and also fluid velocity can be interpolated on structure net point, thus upgrade the position of structure.Except above-mentioned steps, the simulation of fluid structurecoupling can be completed very compactly without the need to extra additional treatments.

In immersed boundary method, the most basic is a bit how to solve the incompressible NS equation of viscosity.The method for solving of this equation is very abundant, such as Fast Fourier Transform (FFT), projecting method or finite volume method.In addition the very welcome Lattice Boltzmann Method that newly-developed gets up also is good selection.The NS method solved with classic method is different, adopting is situated between see that the Lattice Boltzmann Method at visual angle solves be distribution function the Boltzmann equation that meets.In the past during the decade, Lattice Boltzmann Method has been confirmed to be a kind of very powerful computational fluid dynamics instrument, and all obtains successful Application in its branch field.Compare with classic method, Lattice Boltzmann Method be easy to parallel and extra physical influence the concern that advantage obtains more and more scholar such as to facilitate the introduction of.

Since 2003, the microprocessor Design of semi-conductor industry starts towards two Main way development.First multinuclear path relies on the execution efficiency of strengthening serial code and multi-core technology to improve performance.Multinuclear starts from two nuclear models, and along with the doubles of the core that updates each time, the Core i7 processor of such as Intel has four cores.In contrast, above the handling capacity that many core paths having focused on then by concern improves parallel codes execution, compare with central processing unit, but many-core processor will comprise more much smaller core, the Geforce GTX280 graphic process unit of such as NVIDIA company comprises 240 cores.As comparing, the ratio of the peak value operational performance of graphic process unit in 2009 and central processing unit is 10:1, the peak throughput that the Geforce 780Ti of NVIDIA company in 2014 provides is 5500GFLOPS, and the data of the Ivy Bridge framework of Intel Company are 700GFLOPS.So large operational performance gap has allowed increasing scientific research personnel that the computation-intensive part in their programs is transferred to Graphics Processing Unit to perform.

Because such as need directly for chip programming as OpenGL and Direct3D technology, and until programmer in 2006 must visit its core by application programming interface, so the use of graphic process unit is very difficult.These application programming interface limit the type of the application program can write for chip, utilize graphic process unit to run the necessary technology of specific program, and then limit the wide-scale distribution of this technology so only have fraction programmer to grasp.All were released CUDA technology along with NVIDIA company in 2007 and change.Because evaded traditional application programming interface, pass through CUDA, present programmers can with they the class C/C++ language be familiar with come directly for chip programming, even if thus make multiple programming technology also be very easy to left-hand seat for beginner, advance the application of this technology in scientific domain widely.

Because Lattice Boltzmann Method is the explicit scheme be defined on orthogonal grid, and only need the information exchanging adjacent node, so it is very suitable for parallelization in graphic process unit.Be developed so far, the article about Lattice Boltzmann Method parallelization delivered is very abundant, and compared to central processing unit, travelling speed generally can improve two orders of magnitude.Although the parallelization of Lattice Boltzmann Method becomes better and approaching perfection day by day, but the research it being combined the solid coupled problem of analog stream with immersed boundary method does not still start to walk, and does not more have the article of the potential data parallelism of analysis immersed boundary method of thorough to deliver.

Summary of the invention

For the deficiency that prior art exists, the present invention proposes a kind of solid coupled simulation method of efficient stream based on GPU.

In order to solve the problems of the technologies described above, the present invention adopts following technical scheme:

Concrete operation step of the present invention is as follows:

1) in professional pre-processing software, set up the immersed boundary build of structure to be measured and internal flow thereof, and divide border grid and read for program;

2) program reading boundary build data, arrange boundary condition and the procedural model parameter of zoning;

3) Finite element arithmetic border is utilized to be applied to acting force on fluid according to border current time configuration and reference configuration;

4) acting force on border is diffused to peripheral fluid via Dirac function;

5) Navier Stokes equation of external force term is utilized Lattice Boltzmann Method to solve to carry;

6) fluid velocity is interpolated into via Dirac function translational speed border obtaining border equally, and then upgrades boundary position;

7) step 3-6 is repeated, until reach calculating terminal.

Described professional pre-processing software is Gambit or ICEM or ProE.

In the present invention, detailed process is as follows:

One, the interaction between immersed boundary method process fluid and structure is adopted.

Comprise the feature of following several respects specifically:

1), under two-dimensional condition, one dimension wire grid discrete topology body is adopted.Under three-dimensional condition, adopt destructuring triangular mesh discrete topology body.

2) under two-dimensional condition, consider that structure opposing stretches and bending ability simultaneously, be defined in tension on node and resistance to bending density is calculated as follows respectively:

${\left(F_s\right)}_l=\frac{K_s}{{\left({\mathrm{\Δ\α}}_l\right)}^2}\underset{m=1}{\overset{n_f-1}{\Σ}}\left(\right|X_{m+1}-X_m|-{\mathrm{\Δ\α}}_l)\frac{X_{m+1}-X_m}{|X_{m+1}-X_m|}({\mathrm{\δ}}_{ml}-{\mathrm{\δ}}_{m+1,l})\·$

${\left(F_b\right)}_l=\frac{K_b}{{\left({\mathrm{\Δ\α}}_l\right)}^4}\underset{m=2}{\overset{n_f-1}{\Σ}}(X_{m+1}+X_{m-1}-2X_m)(2{\mathrm{\δ}}_{ml}-{\mathrm{\δ}}_{m+1,l}-{\mathrm{\δ}}_{m-1,l})\·$

(F in above formula _s) _l(F _b) _lrepresent respectively to be defined in and be numbered tension on the node of l and unbending lines of force density; K _sand K _bbe respectively tensile index and bending stiffness; Δ α _lfor the former length of the node of numbering l; X _m+1, X _m-1and X _mbe respectively the node coordinate of numbering m+1, m-1 and m; δ _mlfor Kronecker symbol, as m=l, its value is 1, and in all the other situations, value is 0; n _ffor being used for total nodes of discrete topology body.

Under three-dimensional condition, the stress of Elasticity method computation structure body can be used.In the present invention, we only consider the drawing-resistant function of structure.First only there is the contribution of tension in the energy of structure, and the gross energy W of structure can be expressed as W=∫ w ^sdA, wherein w ^sfor drawing of structure can density, and dA is bin area.For the material following broad sense Hooke theorem, it draws and can density can be expressed as:

$w^s=\frac{E}{6}(I_1+\frac{1}{I_2+1}-1),$

Wherein E is elastic modulus, I ₁and I ₂be respectively the first and second strain invariants.Under three-dimensional condition, can prove that strain invariant can be calculated as follows:

${\mathrm{\λ}}_1^2+{\mathrm{\λ}}_2^2=a^2+b^2+c^2,{\mathrm{\λ}}_1^2{\mathrm{\λ}}_2^2=a^2{c}^2$

Wherein λ ₁and λ ₂be along the draw ratio in a pair mutually perpendicular direction in normal plane, and the calculating formula of a, b, c is:

$a=\frac{l}{l_0},$

In above formula, accompanying drawing 1 is shown in the definition of parameter, in above formula, geometric parameter is defined as follows: after structure reference configuration and distortion, configuration chooses two pieces of corresponding Triangular patch, in above-mentioned triangle, choose arbitrarily two limits respectively, on reference configuration, its length of side is respectively l ₀and l ₀', and after deformation on configuration its length of side be respectively l and l ', with be respectively the angle of above-mentioned two limits after reference configuration and distortion on configuration.。

The strain energy of structure can be calculated by above-mentioned steps, is finally defined in the acting force density F be numbered on m node _mthe principle of virtual work can be utilized to calculate:

$F_m=-\frac{\∂W\left(X_m\right)}{\∂X_m}\·$

3) after obtaining structure acting force density, we need to be diffused to structure peripheral fluid node, and by introducing Dirac function, said process can be described by following equation:

$f(x,t)=\underset{m}{\Σ}{F}_m\mathrm{\δ}(x-X_m){\mathrm{\Δs}}_k$

The acting force density that f (x, t) in above formula is subject to for t point x place fluid, respectively just to two and three dimensions condition, Δ s _kbe defined in be numbered arc length on the structure node of k and area respectively, and δ (x-X _m) be Dirac function, generally take following concrete form:

$\mathrm{\&delta;}(x-X_m)=\left\{\begin{array}{cc}\frac{1}{8}(3-2|x-X_m|+\sqrt{1+4|x-X_m|-4{(x-X_m)}^2}),& 0\&le;|x-X_m|<1\\ \frac{1}{8}(5-2|x-X_m|-\sqrt{-7+12|x-X_m|-4{(x-X_m)}^2}),& 1\&le;|x-X_m|<2,\\ 0,& |x-X_m|\&GreaterEqual;2.\end{array}\right.$

In order to obtain the change procedure of structure shape, we need the configuration upgrading structure according to the speed of structure net point, and under known flow field velocity distribution occasion, structure speed will utilize Dirac function by obtaining from flow field interpolation equally:

$\frac{{\mathrm{dX}}_m}{dt}=U_m=\underset{x}{\&Sigma;}u(x,t)\mathrm{\&delta;}(x-X_m)\mathrm{\&Delta;}h$

U in above formula _mfor the speed of structure net point, u (x, t) is the speed of t point x place fluid, and Δ h is mesh spacing.

Two, the incompressible modeling flow field in Lattice Boltzmann Method is adopted.

Lattice Boltzmann Method is a kind of novel method in Fluid Mechanics Computation field.It is often used to numerical solution unsteady flo w can not press NS equation.Its advantage is the ability presenting complicated physical phenomena, and these phenomenons comprise from polyphasic flow up to the chemical reaction between fluid and ambient substance.The method stems from the description of fluid molecule level and can directly the interactional information between molecule and its simple form be embedded into inside model.

Different from the classic method of the variablees such as the velocity pressure in direct solution NS equation, the starting point of Lattice Boltzmann Method be placed on function g (x, ξ, t) on the Boltzmann equation that meets.Function g (x, ξ, t) physical significance is that t finds speed to be positioned at the probability of the molecule of unit interval near ξ near x in unit volume, according to Molecule Motion Theory, under the thin condition approximate with collision term employing BGK of gas, its equation met is Boltzmann equation:

$\frac{\&part;g(x,\mathrm{\&xi;},t)}{\&part;t}+\mathrm{\&xi;}\&CenterDot;\&dtri;g(x,\mathrm{\&xi;},t)+f(x,t)\&CenterDot;\frac{\&part;g(x,\mathrm{\&xi;},t)}{\&part;\mathrm{\&xi;}}=-\frac{1}{\mathrm{\&lambda;}}\&lsqb;g(x,\mathrm{\&xi;},t)-g^{eq}(x,\mathrm{\&xi;},t)\&rsqb;,$

Wherein g ^eq(x, ξ, t) is equilibrium distribution function, for vector differential operator, λ is slack time.By above-mentioned equation discrete in the velocity space of time and space and finite dimensional, we can obtain Lattice Boltzmann equation:

$\begin{array}{c}{g}_{\mathrm{\&alpha;}}\left(x+{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t\right)-g_{\mathrm{\&alpha;}}\left(x,t\right)=-\frac{1}{\mathrm{\&tau;}}\left(g_{\mathrm{\&alpha;}}\left(x,t\right)-g_{\mathrm{\&alpha;}}^{eq}\left(x,t\right)\right)+\\ \left(1-\frac{1}{2\mathrm{\&tau;}}\right)\left\{3\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}-u\left(x,t\right)\&rsqb;+9\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\&CenterDot;u\left(x,t\right)\&CenterDot;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\&rsqb;\right\}\&CenterDot;f\left(x,t\right),\end{array}$

In formula, g _α(x, t) represents that t x place is along discrete velocity reversal ξ _αdistribution function, and g _α(x+ ξ _αδ _t, t+ δ _t) represent t+ δ _tmoment x+ ξ _αδ _tplace is along discrete velocity reversal ξ _αdistribution function. that t x place is along discrete velocity reversal ξ _αequilibrium distribution function.τ is nondimensional slack time, relevant with intermolecular viscosity coefficient.δ _tfor time step.The macroscopic physical quantity that we are concerned about such as density p (x, t), speed etc. can be obtained by each rank velocity moment of distribution function:

$\mathrm{\&rho;}(x,t)=\underset{\mathrm{\&alpha;}}{\&Sigma;}{g}_{\mathrm{\&alpha;}}(x,t),\mathrm{\&rho;}(x,t)u(x,t)=\underset{\mathrm{\&alpha;}}{\&Sigma;}{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}{g}_{\mathrm{\&alpha;}}(x,t)+f(x,t)/2.$

Pressure is then determined by state equation.

Distribution function Evolution represented by Lattice Boltzmann equation can be decomposed into collision and transition process, is shown below.

$g_{\mathrm{\&alpha;}}(x,t)=-\frac{1}{\mathrm{\&tau;}}\left(g_{\mathrm{\&alpha;}}\right(x,t)-g_{\mathrm{\&alpha;}}^{eq}(x,t)),$

g _α(x+ξ _αδ _t,t+δ _t)＝g _α(x,t),

First distribution function will be relaxed to its equilibrium state by collision on node, move thereafter along discrete velocity reversal to adjacent node again.On the whole, if the value of a known upper moment full-field distribution function and macroscopic physical quantity, the step developed to subsequent time is:

1) collision of distribution function, is carried out according to formula.

2), the process of boundary condition.

3), according to formula complete the transition process of distribution function, the distribution function of subsequent time can be obtained.

4), calculated the macroscopic physical quantity of subsequent time according to formula by the distribution function after moving.

Three, the parallel algorithm speed-up computation of graphic based processing unit is adopted.

In the present invention, the CUDA language that NIVIDIA company releases is used to develop the code being adapted at executed in parallel in graphic process unit.When research object can show abundant data parallelism, we just can consider workload partition to become much non-interfering detached process, and these processes are carried out simultaneously thus increased work efficiency, shorten man-hour.In CUDA, the concept relevant to said process is above thread, and thread is responsible for the process of independent a certain item process, and task can be divided into how many processes, just has how many threads corresponding with it.According to the actual demand of problem, thread can be organized into one dimension, two dimension or three-dimensional thread block, and thread block can be organized into the thread grid of one dimension or two dimension.Once thread is organized into orderly structure, so also just provides the mechanism distinguishing thread, thus different work can be distributed to different threads by thread number in a program.

Storage space and the Installed System Memory of modern GPU are separate, so all needed data to transfer to video memory from Installed System Memory before carrying out any calculating with GPU.Access efficiency due to video memory depends on the mode that video memory is accessed significantly, when our access be continuous print a slice data in video memory time, the peak bandwidth of video memory access can be reached, so we need the standardization as far as possible ensureing video memory access when programming.

CUDA paralleling tactic of the present invention has following characteristics specifically:

1), thread structure aspect, thread block is divided into one-dimentional structure, and thread grid is divided into a peacekeeping two-dimensional structure respectively under two and three dimensions condition, and each thread is responsible for the calculating of a net point specially.

2), data structure aspect, adopt one-dimension array storage and distribution function and macroscopic physical quantity, the distribution function along different directions is stored by different arrays.

3), in transition process in order to avoid graphic process unit video memory access non-alignment behavior, utilize shared storage that transition process is divided into two steps, first step distribution function will utilize shared storage along the migration of direction, thread block place, and distribution function edge is moved perpendicular to direction, thread block place by video memory by second step again.Because shared storage size is limited, so need extra kernel function to exchange the distribution function exceeding thread block border in first step migration.

4), under three-dimensional condition, the topological relation storage policy of structure surface triangles grid is: provide all Triangular patch numberings of this node as one of summit according to grid node.Embodiment is in a program definition structure body array, is used for preserving the triangle number information corresponding to node.

Compared with prior art, the present invention has following advantage and beneficial effect:

1, the present invention can process the structure of arbitrary shape, and do not need reconstruct bottoms stream grid in computation process, counting yield is high.

2, in the present invention, structure and fluid interaction implementation procedure are simply direct, and need not arrange boundary condition at interface place, the globality of flow field calculation is better.

3, the Lattice Boltzmann Method algorithm that the present invention is used is simple, can well simulate NS equation when close can not pressure with certain precision.

4, traditional method for numerical simulation energy that needs to cost a lot of money goes to separate Poisson equation and just can obtain pressure, and the pressure state equation in Lattice Boltzmann Method is tried to achieve, more directly and convenient.

5, the present invention has very high concurrency, its evolution mainly local computing that synchronously performs of some classes, and these computings all can show and carry out, and it can directly be realized on parallel machine, and the communication industry between multiprocessing is very easily optimized.Therefore the present invention is less by the restriction of computer capacity and speed, can simulate fairly large high-order flow field.

6, algorithm stability of the present invention is better, to the less-restrictive of structure physical parameter span, is suitable for simulating practical problems.

Accompanying drawing explanation

Fig. 1 is the parametrization method for expressing of structure grid of the invention process.Triangular element when Fig. 1 (a) is structure beinthebalancestate, Fig. 1 (b) is the triangular element after structure distortion, and Fig. 1 (c) utilizes translation to make the above two certain situation be located on the same line.Namely the stress of structure solve by the displacement of distortion front-rear triangular shape unit corresponding vertex.

Fig. 2 is the film center pressure change procedure curve of embodiment 1 Elastic film when reverting to circular state from oval state.

Fig. 3 is the time-varying process that embodiment 2 Elastic film fixes 2 a, b distance thin film center length when petal state reverts to circular state;

Fig. 4 is the time-varying process of round and elastic film spheroid central authorities' pressure and spheroid volume when stretching recovering state is nature circular state outward from entirety in embodiment 3;

The time-varying process of Fig. 5 is embodiment 4 Elastic film when being spheroidal state from spheroid recovering state semi-major axis;

Fig. 6 be in embodiment 5 spheroid elastic film from entirety outward stretching recovering state be nature spheroidal state process, balloon is provided with a cut, in motion process after balloon breaks, the ordinate time-varying process curve of a, b, c that spheroid is arranged 3.

Fig. 7 inflates in embodiment 6 in the spheroid elastic film of initial equilibrium state, make it expand, and in expansion process, balloon is provided with a cut, the time-varying process of a, b, d of structurally arranging 3 distance initial equilibrium state spheroid thin film center point length obtained and the time-varying process of the fluid velocity size at film cut center.

Embodiment

Below by embodiment, illustrate outstanding feature of the present invention and marked improvement further, be only the present invention is described and never limits the present invention.

Step of the present invention is as follows:

1) in professional pre-processing software, set up the immersed boundary build of structure to be measured and internal flow thereof, and divide border grid and read for program;

2) program reading boundary build data, arrange boundary condition and the procedural model parameter of zoning;

3) Finite element arithmetic border is utilized to be applied to acting force on fluid according to border current time configuration and reference configuration;

4) acting force on border is diffused to peripheral fluid via Dirac function;

5) Navier Stokes equation of external force term is utilized Lattice Boltzmann Method to solve to carry;

6) fluid velocity is interpolated into via Dirac function translational speed border obtaining border equally, and then upgrades boundary position;

7) step 3-6 is repeated, until reach calculating terminal.

Described professional pre-processing software is Gambit or ICEM or ProE.

In described step 3:

(1) when immersed boundary build is two dimension, adopt one dimension wire grid discrete topology body, consider that structure opposing stretches and bending ability simultaneously, be defined in tension on node and resistance to bending density is calculated as follows respectively:

${\left(F_s\right)}_l=\frac{K_s}{{\left({\mathrm{\&Delta;\&alpha;}}_l\right)}^2}\underset{m=1}{\overset{n_f-1}{\&Sigma;}}\left(|X_{m-1}-X_m|-{\mathrm{\&Delta;\&alpha;}}_l\right)\frac{X_{m+1}-X_m}{|X_{m+1}-X_m|}\left({\mathrm{\&delta;}}_{ml}-{\mathrm{\&delta;}}_{m+1,l}\right).$

${\left(F_b\right)}_l=\frac{K_b}{{\left({\mathrm{\&Delta;\&alpha;}}_l\right)}^4}\underset{m=2}{\overset{n_f-1}{\&Sigma;}}(X_{m+1}+X_{m-1}-2X_m)(2{\mathrm{\&delta;}}_{ml}-{\mathrm{\&delta;}}_{m+1,l}-{\mathrm{\&delta;}}_{m-1,l})\&CenterDot;$

(F in above formula _s) _l(F _b) _lrepresent respectively to be defined in and be numbered tension on l node and unbending lines of force density, K _sand K _bbe respectively tensile index and bending stiffness, Δ α _lfor the former length of the node of numbering l, X _m+1, X _m-1and X _mbe respectively the node coordinate of numbering m+1, m-1 and m, and δ _mlfor Kronecker symbol, as m=l, its value is 1, and in all the other situations, value is 0; n _ffor being used for total nodes of discrete topology body;

(2) when immersed boundary is three-dimensional, destructuring triangular mesh discrete topology body is adopted.Construction is numbered the acting force density F on m node _mcomputing formula as follows:

$F_m=-\frac{\&part;W\left(X_m\right)}{\&part;X_m}$

In above formula: W is the gross energy of structure, for partial differential symbol, X _mfor the coordinate of node m; W=∫ w ^sdA, w ^sfor drawing of structure can density, and dA is bin area, wherein E is elastic modulus, I ₁and I ₂be respectively the first and second strain invariants.Can prove that strain invariant can be calculated as follows:

${\mathrm{\&lambda;}}_1^2+{\mathrm{\&lambda;}}_2^2=a^2+b^2+c^2,{\mathrm{\&lambda;}}_1^2{\mathrm{\&lambda;}}_2^2=a^2{c}^2$

Wherein λ ₁and λ ₂be the draw ratio along a pair vertical direction in normal plane respectively, and the calculating formula of a, b, c is: $a=\frac{l}{l_0},$

In described step 4 and step 6, acting force diffusion equation is:

$f(x,t)=\underset{m}{\&Sigma;}{F}_m\mathrm{\&delta;}(x-X_m){\mathrm{\&Delta;s}}_k$

The acting force density that f (x, t) in above formula is subject to for t point x place fluid, respectively for two and three dimensions condition, Δ s _kbe defined in be numbered arc length on the structure node of k and area respectively, and δ (x-X _m) be Dirac function, generally take following concrete form:

$\mathrm{\&delta;}(x-X_m)=\left\{\begin{array}{cc}\frac{1}{8}(3-2|x-X_m|+\sqrt{1+4|x-X_m|-4{(x-X_m)}^2}),& 0\&le;|x-X_m|<1\\ \frac{1}{8}(5-2|x-X_m|-\sqrt{-7+12|x-X_m|-4{(x-X_m)}^2}),& 1\&le;|x-X_m|<2,\\ 0,& |x-X_m|\&GreaterEqual;2.\end{array}\right.$

In order to obtain the change procedure of structure shape, we need the configuration upgrading structure according to the speed of structure net point, and under known flow field velocity distribution occasion, structure speed will utilize Dirac function by obtaining from flow field interpolation equally:

$\frac{{\mathrm{dX}}_m}{dt}=U_m=\underset{x}{\&Sigma;}u(x,t)\mathrm{\&delta;}(x-X_m)\mathrm{\&Delta;}h$

U in above formula _mfor the speed of structure net point, u (x, t) is the speed of t point x place fluid, and Δ h is mesh spacing.

In described step 5, Lattice Boltzmann equation is:

$\begin{array}{c}{g}_{\mathrm{\&alpha;}}\left(x+{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t\right)-g_{\mathrm{\&alpha;}}\left(x,t\right)=-\frac{1}{\mathrm{\&tau;}}\left(g_{\mathrm{\&alpha;}}\left(x,t\right)-g_{\mathrm{\&alpha;}}^{eq}\left(x,t\right)\right)+\\ \left(1-\frac{1}{2\mathrm{\&tau;}}\right)\left\{3\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}-u\left(x,t\right)\&rsqb;+9\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\&CenterDot;u\left(x,t\right)\&CenterDot;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\&rsqb;\right\}\&CenterDot;f\left(x,t\right),\end{array}$

In formula, g _α(x, t) represents that t x place is along discrete velocity reversal ξ _αdistribution function, and g _α(x+ ξ _αδ _t, t+ δ _t) be then t+ δ _tmoment x+ ξ _αδ _tplace is along discrete velocity reversal ξ _αdistribution function. that t x place is along discrete velocity reversal ξ _αequilibrium distribution function.τ is nondimensional slack time, relevant with intermolecular viscosity coefficient.δ _tfor time step.

Embodiment 1

Experimental subjects is the film center pressure change procedure of elastic film when reverting to circular state from oval state;

(1) zoning is chosen for the square that the length of side is 1.0m, discrete be 200 × 200 orthogonal cartesian grid.Structure is that semi-major axis is respectively r _a=0.75m and r _bthe ellipse of=0.5m.The circle of the equilibrium state of structure to be radius be 0.5149m, the number of nodes of discrete topology is 1200.

(2) structure elastic modulus is set to 0.55Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

Fig. 2 is the film center pressure change procedure curve of elastic film when reverting to circular state from oval state.As seen from the figure circular membrane to the relaxation of its equilibrium state by the variation of the ripple that causes its central fluid pressure to be decayed gradually with an amplitude, and central pressure convergence one fixed value the most at last;

Embodiment 2

Experimental subjects is the time-varying process that elastic film fixes 2 a, b distance thin film center length when petal state reverts to circular state;

(1) zoning is chosen for the square that the length of side is 2.0m, discrete be 200 × 200 orthogonal cartesian grid.The six lobe films of structure for being described by polar equation ρ=0.5 [1+0.4cos (6 θ)], the number of nodes of discrete topology is 2200.The circle of the equilibrium state of structure to be radius be 0.518m.

(2) structure elastic modulus is set to 0.55Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

Embodiment 3

Experimental subjects is the time-varying process of round and elastic film spheroid central authorities' pressure and spheroid volume when stretching recovering state is nature circular state outward from entirety;

(1) zoning is chosen for the cube that the length of side is 1.0m, discrete be 64 × 64 × 64 orthogonal cartesian grid.In pre-processing software Gambit, set up spheroid build, its radius is 0.2m, and divides gore grid on its surface, exports grid file, reads for master routine.

(2) master routine reads the structure grid file that pre-processing software exports, and structure elastic modulus is set to 0.55Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

Embodiment 4

The time-varying process of experimental subjects is elastic film when being spheroidal state from spheroid recovering state semi-major axis;

(1) zoning is chosen for the cube that the length of side is 1.0m, discrete be 64 × 64 × 64 orthogonal cartesian grid.In pre-processing software Gambit, set up ellipsoid build, its semi-major axis is respectively r _a=0.3m, r _b=0.24m, r _c=0.18m, and divide gore grid on its surface, export grid file, read for master routine.

(2) master routine reads the structure grid file that pre-processing software exports, and structure elastic modulus is set to 0.15Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

As seen from the figure, the semi-major axis of spheroid alternately rises and declines, and the amplitude of its fluctuation is also constantly decaying.Because under surface tension effects, any film all will level off to perfectly spherical, so can see the fixed value that semi-major axis is tending towards identical the most at last.

Embodiment 5

Experimental subjects be spheroid elastic film from entirety outward stretching recovering state be nature spheroidal state process; Wherein measuring point c is arranged near film cut, and a point, b point are respectively 2 end points diametrically, and b point is arranged in balloon far right end, and a point is arranged in the low side of balloon.

(1) zoning is chosen for the cube that the length of side is 1.0m, discrete be 200 × 200 × 200 orthogonal cartesian grid.In pre-processing software Gambit, set up the spheroid build that radius is r=0.156m, and divide gore grid on its surface, export grid file, read for master routine.

(2) master routine reads the structure grid file that pre-processing software exports, and structure elastic modulus is set to 0.15Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.Add prestress to film before calculating starts, make the surface tension of film and film inside and outside differential pressure meet laplace's theorem.Open on the top of film the minor break that Radius is 0.034m afterwards, fluid will be gushed out from spout under the effect of film inside and outside differential pressure.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

As shown in Figure 6, the motion of balloon different piece is nonsynchronous in fact, and its repercussion effect of part the closer to cut embodies more early, and therefore balloon will be pressed into elliptical shape in the early stage.

Embodiment 6

Experimental subjects is inflate in the spheroid elastic film of initial equilibrium state, make it expand, and in expansion process, balloon is provided with a cut, the time-varying process of the time-varying process of a, b, d that spheroid elastic film is arranged 3 distance initial equilibrium state ball centre length and the fluid velocity size at film cut center, wherein a point is arranged near film cut, and b point is arranged in the far right end of balloon, and d point is arranged near balloon nozzle.

(1) zoning is chosen for the cube that the length of side is 1.0m, discrete be 264 × 264 × 264 orthogonal cartesian grid.In pre-processing software Gambit, set up structural style, structure is the assembly of cylindrical nozzle and spheroid balloon, and wherein nozzle is rigid body, and balloon is fexible film.Wherein the radius of balloon film is 0.24m, and the radius of nozzle cross-section is 0.15m, and long is 0.086m.Utilize Gambit to divide gore grid at body structure surface, export grid file, read for master routine.

(2) master routine reads the structure grid file that pre-processing software exports, and structure elastic modulus is set to 0.15Gpa, and fluid density is 1.0g/cm ³, kinematic viscosity is 1.0 × 10 ^-6m ²/ s, fluid original state is static, and border is first kind pressure boundary condition.During beginning from nozzle with 0.02m ³the flow of/s injects liquid to balloon, and removes balloon front end one small pieces of material to a certain extent afterwards at ballooning.

(3) setting exports once result every specific time period, working procedure, till flow field and structure finally settle out, can obtain the interaction process of flow field and structure.

As shown in Figure 7, balloon at cut in a flash by back seat process of short duration for experience one, afterwards in order to be drained via cut by balloon interior fluid as early as possible, balloon longitudinally will be drawn as strip, and final balloon is returned to again original equilibrium configuration after inside and outside pressure tends to balance.

Claims (5)

1., through immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method that GPU accelerates, it is characterized in that: step is as follows:

(1) in professional pre-processing software, set up the immersed boundary build of structure to be measured and internal flow thereof, and divide border grid and read for program;

(2) program reading boundary build data, arrange boundary condition and the procedural model parameter of zoning;

(3) Finite element arithmetic border is utilized to be applied to acting force on fluid according to border current time configuration and reference configuration;

(4) acting force on border is diffused to peripheral fluid via Dirac function;

(5) Navier Stokes equation of external force term is utilized Lattice Boltzmann Method to solve to carry;

(6) fluid velocity is interpolated into via Dirac function translational speed border obtaining border equally, and then upgrades boundary position;

(7) step 3-6 is repeated, until reach calculating terminal.

2., as claimed in claim 1 through immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method that GPU accelerates, it is characterized in that: described professional pre-processing software is Gambit or ICEM or ProE.

3., as claimed in claim 1 through immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method that GPU accelerates, it is characterized in that: in described step 3:

(1) when immersed boundary build is two dimension, adopt one dimension wire grid discrete topology body, consider that structure opposing stretches and bending ability simultaneously, be defined in tension on node and resistance to bending density is calculated as follows respectively:

$\begin{array}{c}{\left(F_s\right)}_l=\frac{K_s}{{\left({\mathrm{\&Delta;\&alpha;}}_l\right)}^2}\underset{m=1}{\overset{n_f-1}{\&Sigma;}}\left(\right|X_{m+1}-X_m|-{\mathrm{\&Delta;\&alpha;}}_l)\frac{X_{m+1}-X_m}{|X_{m+1}-X_m|}({\mathrm{\&delta;}}_{ml}-{\mathrm{\&delta;}}_{m+1,l}).\\ {\left(F_b\right)}_l=\frac{K_b}{{\left({\mathrm{\&Delta;\&alpha;}}_l\right)}^2}\underset{m=2}{\overset{n_f-1}{\&Sigma;}}(X_{m+1}+X_{m-1}-2X_m)(2{\mathrm{\&delta;}}_{ml}-{\mathrm{\&delta;}}_{m+1,l}-{\mathrm{\&delta;}}_{m-1,l}).\end{array}$

(F in above formula _s) _l(F _b) _lrepresent respectively to be defined in and be numbered tension on l node and unbending lines of force density; K _sand K _bbe respectively tensile index and bending stiffness; Δ α _lfor the original arc length of the node of numbering l; X _m+1, X _m-1and X _mbe respectively the node coordinate of numbering m+1, m-1 and m; δ _mlfor Kronecker symbol, as m=l, its value is 1, and in all the other situations, value is 0; n _ffor being used for total nodes of discrete topology body;

(2) when immersed boundary build is three-dimensional, destructuring triangular mesh discrete topology body is adopted,

Construction is numbered the acting force density F on m node _mcomputing formula as follows:

$F_m=-\frac{\&part;W\left(X_m\right)}{\&part;X_m}.;$

Wherein: W is the gross energy of structure, for partial differential symbol, X _mfor the coordinate of node m; The calculating formula of structure gross energy is:

W＝∫w ^sdA

DA in formula is bin area, and w ^sfor drawing of structure can density, its calculating formula is:

wherein E is elastic modulus, I ₁and I ₂be respectively the first and second strain invariants, can prove that strain invariant can be calculated as follows:

${\mathrm{\&lambda;}}_1^2+{\mathrm{\&lambda;}}_2^2=a^2+b^2+c^2,{\mathrm{\&lambda;}}_1^2{\mathrm{\&lambda;}}_2^2=a^2{c}^2$

Wherein λ ₁and λ ₂be the draw ratio along a pair vertical direction in normal plane, and the calculating formula of a, b, c is:

In above formula, geometric parameter is defined as follows: after structure reference configuration and distortion, configuration chooses two pieces of corresponding Triangular patch, in above-mentioned triangle, choose arbitrarily two limits respectively, on reference configuration, its length of side is respectively l ₀with l ' ₀, and after deformation on configuration its length of side be respectively l and l ', with be respectively the angle of above-mentioned two limits after reference configuration and distortion on configuration.

4. as claimed in claim 1 through immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method that GPU accelerates, it is characterized in that: the acting force diffusion equation in described step 4 and step 6 is:

$f(x,t)=\underset{m}{\&Sigma;}{F}_m\mathrm{\&delta;}(x-X_m){\mathrm{\&Delta;s}}_k$

The acting force density that f (x, t) in above formula is subject to for t point x place fluid, F _mfor being numbered the acting force density on m node, respectively for two and three dimensions condition, Δ s _kbe defined in be numbered arc length on the structure node of k and area respectively, and δ (x-X _m) be Dirac function, generally take following concrete form:

$\mathrm{\&delta;}(x-X_m)=\left\{\begin{array}{cc}\frac{1}{8}(3-2|x-X_m|+\sqrt{1+4|x-X_m|-4{(x-X_m)}^2}),& 0\&le;|x-X_m|<1\\ \frac{1}{8}(5-2|x-X_m|-\sqrt{-7+12|x-X_m|-4{(x-X_m)}^2}),& 1\&le;|x-X_m|<2,\\ 0,& |x-X_m|\&GreaterEqual;2.\end{array}\right.$

In order to obtain the change procedure of structure shape, we need the configuration upgrading structure according to the speed of structure net point, and under known flow field velocity distribution occasion, structure speed will utilize Dirac function by obtaining from flow field interpolation equally:

$\frac{{\mathrm{dX}}_m}{dt}=U_m=\underset{x}{\&Sigma;}u(x,t)\mathrm{\&delta;}(x-X_m)\mathrm{\&Delta;}h$

X in above formula _mand U _mbe respectively coordinate and the speed of structure net point, u (x, t) is the speed of t point x place fluid, and Δ h is mesh spacing.

5., as claimed in claim 1 through immersed boundary-Lattice Boltzmann fluid structurecoupling analogy method that GPU accelerates, it is characterized in that: in described step 5, Lattice Boltzmann equation is:

$\begin{array}{c}{g}_{\mathrm{\&alpha;}}(x+{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}{\mathrm{\&delta;}}_t,t+{\mathrm{\&delta;}}_t)-g_{\mathrm{\&alpha;}}(x,t)-\frac{1}{\mathrm{\&tau;}}\left(g_{\mathrm{\&alpha;}}\right(x,t)-g_{\mathrm{\&alpha;}}^{eq}(x,t))+\\ (1-\frac{1}{2\mathrm{\&tau;}})\{3\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}-u(x,t)\&rsqb;+9\&lsqb;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\&CenterDot;u(x,t)\&rsqb;\&CenterDot;{\mathrm{\&xi;}}_{\mathrm{\&alpha;}}\}\&CenterDot;f(x,t),\end{array}$

In formula, g _α(x, t) represents that t x place is along discrete velocity reversal ξ _αdistribution function, and g _α(x+ ξ _αδ _t, t+ δ _t) be then t+ δ _tmoment x+ ξ _αδ _tplace is along discrete velocity reversal ξ _αdistribution function, that t x place is along discrete velocity reversal ξ _αequilibrium distribution function, τ is nondimensional slack time, relevant with intermolecular viscosity coefficient, δ _tfor time step.