CN104636566A - Dynamic mesh numerical solution method based on modified fluid motion equation - Google Patents

Abstract

The invention relates to a dynamic mesh numerical solution method based on a modified fluid motion equation, and belongs to the technical field of optimization methods. The method comprises the following steps: introducing an original non-linear convective term inside a concept deconstruction equation of convection speed into a speed equation on the basis of a classic cauchy equation of motion; deriving a linear modified fluid motion equation so as to establish a dynamic mesh numerical iteration solver special for the modified fluid motion equation under a finite volume discretization method. The linear characteristic of the dynamic mesh numerical iteration solver achieves obvious advantage on solving the problem of an internal flow field of a rotary machine, namely the speed item coefficient a of an equation to be solved does not participate in the cyclical iteration process of a PISO algorithm any longer, the equation can be solved by only utilizing a dynamic mesh equation based on the convection speed w, flux change caused by mesh motion can be directly obtained through the convection speed, and thus a numerical calculation process becomes an iteration solution process of a simple speed item U and a pressure item p. The dynamic mesh numerical solution method disclosed by the invention has the characteristics of high operation efficiency, short calculating time and low operation cost aiming at the design calculation of the internal flow field of the rotary machine.

Description

A kind of dynamic mesh method of value solving based on revising fluid motion equation

Technical field

The present invention relates to a kind of dynamic mesh method of value solving based on revising fluid motion equation, belonging to optimization method technical field.

Background technology

Due to the active demand in the industrial research fields such as fluid machinery, the developing rapidly of computer technology over especially past more than 30 years, makes Fluid Mechanics Computation (CFD) achieve very large achievement.Today is, in the research and development that have been widely applied to rotating machinery of the CFD technology of representative and design process, achieve the achievement attracted people's attention with solving NS equations.But along with going deep into gradually of technical development, CFD is also faced with increasing puzzlement, especially on the Solve problems of NS equation.This partial differential equation, after finite volume method is discrete, often generates a numerous and jumbled Groebner Basis.Mathematically, the process of this system of equations, not only increases the difficulty solving CFD problem, also seriously hinders the research and development of relevant numerical simulation method simultaneously.In addition, NS equation also has harsh requirement to the quality of grid, crosses the load that high-quality grid seriously can increase computing machine undoubtedly.These superelevation that all result in CFD analog computation are consuming time, thus have impact on the design efficiency of fluid machinery, add its R&D costs.

Classical NS fluid motion equation can simply be described to: it solves nonlinear iteration steady-state algorithm and often carries out a speed iteration (U _piteration) matrix H (U) all will upgrade, until speed and pressure reach current time step or time step convergence.Use these convergence numerical value, substitute into and recalculate coefficient a _pand a _ncarry out the next time step of initialization.As speed U _presidual error reach the calculating condition of convergence, whole calculating convergence.In this approach, the non-linear and rate-pressure coupled problem of NS equation is all solved by iterative technique.

The nineties in last century, Scale invariant type mechanics theory applies in the middle of fluid mechanics by professor Sohrab of Northwestern Univ USA, and the relationship of system of linear equations conducts in-depth research, and proposes Scale invariant type fluid motion equation concept.A large amount of hydrodynamic calculations and fluid model analytic solution all demonstrate science, accuracy and the linearity benefits that this equation Fluid Mechanics problem carries out solving.Although the appearance of this linear equation concept, the CFD technical barrier that some are caused by NS equation can be solved in theory, but be subject to the restriction of conventional analytic method, the research of Sohrab professor, substantially rest on theoretical research stage, be difficult to expand in more complicated industrial machine flow field model.How finding one and accurately can either solve update equation, can the numerical operation method of this equation linear feature of high efficiency performance will be the key improving this equation range of application further again.

Summary of the invention

Technical matters: the object of the invention is to overcome in prior art the deficiency processing rotating machinery flow field CFD and calculate, improve iterative algorithm and the mesh motion technology of solving equation in existing CFD technology, a kind of dynamic mesh method of value solving based on revising fluid motion equation is provided, can meet solve the prerequisite of accuracy under greatly optimize the iterative algorithm of CFD numerical simulation, the non-linear iterative that NS equation transient state classical under rotating coordinate system of comparing solves, for rotating machinery flow field designing and calculating, there is operation efficiency high, short and the operating cost computing time advantage such as low.

Technical scheme: the present invention is based on the dynamic mesh method of value solving revising fluid motion equation, comprise the following steps:

(1) built the physical arrangement model of rotating machinery design by computing machine, complete dynamic remeshing, based on Scale invariant type statistical parsing, conditions setting convection velocity w, Grid Mobile speed U _b, kinematic viscosity ν, model cootrol body volume V, model cootrol bulk area s initiation parameter, according to rotating machinery designing requirement setting-up time step-length, carry out loop iteration calculating by setting-up time step-length:

According to dynamic flow correction fluid motion equation:

First absolutization process is carried out to flux φ, then dynamic mesh displacement x is determined, and then flux φ is revised, complete the relative processing of the relative Grid Mobile speed of flux φ; In formula: q _φfor the pressure p solved;

(2) spot speed term coefficient a to be solved is calculated _pwith solve neighbor point speed term coefficient a _n, wherein, w _pfor the convection velocity at P point place, U _bPfor the Grid Mobile speed at P point place, w _nfor the convection velocity of P point adjacent place each point, U _bNfor the Grid Mobile speed of adjacent place each point;

(3) iterative process of PISO algorithm circulation is carried out, first compute matrix H (U): wherein U _nrepresent the speed of P point adjacent place each point;

Then calculating pressure p: the equation of momentum at corresponding P point place is: rate equation is: $U_P=\frac{H\left(U\right)}{a_P}-\frac{\▿p}{a_P}$ , pressure equation is: $\▿\·(\frac{1}{a_P}\▿p)=\underset{f}{\mathrm{\Σ}}A\·{\left(\frac{H\left(U\right)}{a_P}\right)}_f,$ Pressure p is the q in equation _φ;

Further relative processing is carried out to flux φ;

To speed U _psolve, wherein, U _pds=φ;

(4) speed U is judged _pwhether restrain,

When PISO algorithm loop iteration evaluation tends towards stability, then judge speed U _pmeet the condition of convergence, then enter next step operation program;

When loop iteration evaluation still changes instability, then judge speed U _pdo not meet the condition of convergence, then repeat step (3), then carry out the iterative computation of PISO algorithm circulation, until speed U _ptill meeting the condition of convergence;

(5) judge to calculate whether reach time step,

As speed U _pmeet the condition of convergence, then carry out the judgement that time step requires whether to reach, if reached, computing terminates, and extracts operation result and provides design parameter for rotating machinery design;

If do not reach set time step, then repeat step (1), then carry out the calculating of iterative loop, till the requirement time step of computing reaches setting value.

Described setting-up time step-length is determined according to actual demands such as the rotating speeds of different rotary machinery.

Beneficial effect: revise the application of fluid motion equation at field of fluid mechanics to explain, the present invention is first from low speed boundary layer flow, compare by analytic solution and well-known Nikuradse Classic Experiments result, and be aided with Modern Laser experimental result, fully demonstrate the accuracy revising fluid motion equation processing layer flow model.Then, in turn give other several flow field model correction fluid motion equation analytic solution with hyperspin feature, and compare with experimental result, be obtained for and well coincide.Based on the proposition of above-mentioned linear revise fluid motion equation, use its update equation to set up the mathematical model of interior flow field in rotating machinery, and adopt the method for numerical simulation to solve.Specifically, under the environment of open source software OpenFOAM, one is specific to the dynamic mesh iterative numerical solver revising fluid motion equation and is established.Convective term in equation is by after discrete, the linear feature of this equation fully shows it and is processing the advantage in rotating machinery flow field problem, namely obtained system of equations is a system of linear equations, wherein speed term coefficient a does not need to participate in the circulation of PISO algorithm, can directly be solved with being coupled of Dynamic mesh by convection velocity w, and the variations of flux that mesh motion causes can directly be obtained by convection velocity, therefore the solution procedure of system of equations is speed term U and the iterative process that is coupled of pressure item p by abbreviation.

And be roughly based on the calculating process of non-linear NS equation numerical value iterative: at moment t, first carry out flow absolutization, dynamic mesh displacement determined, the process of the relativization of flow correction, flow, then speed term coefficient a is solved via the interative computation in a upper moment t-1, then the implicit expression computing of solver commencing speed U, and based on this result, solve pressure p, and then solve U _p, work as U _pduring convergence, judge to calculate whether reach time step, the iterative loop reaching rear whole time step calculates and terminates.The numerical solution calculating process of contrast two equations, apparent, the operational method revising fluid motion equation iterative numerical saves the iteration cycle process of coefficient a in time step.Fluid motion equation dynamic mesh method of value solving is revised in realization of the present invention, can save operation time in essence from computing, promotes the operation efficiency of iterative numerical solver.Convection velocity w in described correction fluid motion equation _βbe implemented in upper level Scale Model β+1, under current scale β, can be counted as the constant term being similar to coefficient of viscosity, therefore update equation is a linear equation, the nonlinear NS equation compared to classics: the solution procedure of update equation will be greatly simplified.Utilize high-performance cloud calculation server, take OpenFOAM as implementation platform, the program of increasing income is write under limited bulk discrete method environment, set up an exclusive iterative numerical solver to carry out dynamic mesh and solve, the non-linear iterative that NS equation transient state classical under rotating coordinate system of comparing solves, do not need to participate in the circulation of PISO algorithm at this iterative process medium velocity term coefficient a, can directly be solved with being coupled of Dynamic mesh by convection velocity w, and the variations of flux that mesh motion causes can directly be obtained by convection velocity, therefore the solution procedure of system of equations is speed term U and the iterative process that is coupled of pressure item p by abbreviation.

This dynamic mesh numerical solution device can be applied to simulation centrifugal pump rotor stator flow coupling phenomenon, and simulation and the relevant modern PIV Laser Experiments of its analog result and traditional NS solver have desirable result identical.Added up by a large amount of simulated experiment, for above CFD model, revise fluid motion equation solver has only accounted for NS solver operation time about 2/3, show the clear superiority of its linear feature on Dynamic mesh uses.

Method for solving of the present invention same meet solve the prerequisite of accuracy under, in contrast to classical Navier Stokes equation iterative numerical method for solving, have operation efficiency high, computing time is short, the feature that operating cost is low, can promote and be applied to the industrial research fields such as rotating machinery.

Accompanying drawing explanation

Fig. 1 is correction fluid motion equation dynamic mesh numerical solution algorithm flow chart of the present invention;

Fig. 2 is centrifugal pump flow field stress and strain model schematic diagram;

Fig. 3 (a), (b) be respectively NS solver, revise the ratio curve comparison diagram of solver, experimental measurements relative velocity under blade radial gap ratio condition vertical and horizontal component and peripheral speed;

Fig. 4 is the iteration time comparison diagram of existing NS equation solver and update equation solver of the present invention.

Embodiment

Below in conjunction with accompanying drawing, one embodiment of the present of invention are further described:

Dynamic solution method based on revising fluid motion equation linear iteration of the present invention, it is the Scale invariant pattern type based on Scale invariant type mechanics theory and broad sense, progressively derive the correction fluid motion equation solved for potential flows grid values, fully using on FVM discrete method basis, complete the correction fluid motion equation analytic solution modeling work under dynamic grid environment.With reference to the technical characterstic of PISO algorithm typical algorithm, and consider the solution such as explicit difference, implicit difference skill and art of programming, under LINUX environment, develop software to increase income as operating platform, the advantage of linear feature in numerical simulation calculation on algorithm of fluid motion equation is revised in exploitation, simplifies the calculating process of solver.On this basis, theoretical validation is carried out by the flow model model analyzing method research of MATLAB, and carry out case verification further by large amount of complex three-dimensional blades machinery flow field, and compare computing time with NS equation numerical value solver iterative algorithm, demonstrate the characteristic and advantage of this linear iterative algorithm.

As shown in Figure 1, realize the dynamic mesh numerical solution algorithm flow chart revising fluid motion equation, its algorithm steps is:

(1) built the physical arrangement model of rotating machinery design by computing machine, complete dynamic remeshing, based on Scale invariant type statistical parsing, conditions setting convection velocity w, Grid Mobile speed U _b, kinematic viscosity ν, model cootrol body volume V, model cootrol bulk area s initiation parameter, according to rotating machinery designing requirement setting-up time step-length, carry out loop iteration calculating by setting-up time step-length:

According to dynamic flow correction fluid motion equation:

Q _φfor the pressure p solved, first absolutization process is carried out to flux φ, then dynamic mesh displacement x is determined, and then flux φ is revised, complete the relative processing of the relative Grid Mobile speed of flux φ;

(2) spot speed term coefficient a to be solved is calculated _pwith solve neighbor point speed term coefficient a _n, wherein, $a_P=2v\frac{1}{\mathrm{\Δx}}+(w_P-U_{\mathrm{bP}}),a_N=v\frac{1}{\mathrm{\Δx}}+(w_N-U_{\mathrm{bN}}),$ W _pfor the convection velocity at P point place, U _bPfor the Grid Mobile speed at P point place, w _nfor the convection velocity of P point adjacent place each point, U _bNfor the Grid Mobile speed of adjacent place each point;

(3) iterative process of PISO algorithm circulation is carried out, first compute matrix H (U): wherein U _nrepresent the speed of P point adjacent place each point;

Then calculating pressure p: the equation of momentum at corresponding P point place is: rate equation is: $U_P=\frac{H\left(U\right)}{a_P}-\frac{\▿p}{a_P}$ , pressure equation is: $\▿\·(\frac{1}{a_P}\▿p)=\underset{f}{\mathrm{\Σ}}A\·{\left(\frac{H\left(U\right)}{a_P}\right)}_f,$ Pressure p is the q in equation _φ;

Further relative processing is carried out to flux φ;

To speed U _psolve, wherein, U _pds=φ;

(4) speed U is judged _pwhether restrain,

When PISO algorithm loop iteration evaluation tends towards stability, then judge speed U _pmeet the condition of convergence, then enter next step operation program;

When loop iteration evaluation still changes instability, then judge speed U _pdo not meet the condition of convergence, then repeat step (3), then carry out the iterative computation of PISO algorithm circulation, until speed U _ptill meeting the condition of convergence;

(5) judge to calculate whether reach time step,

As speed U _pmeet the condition of convergence, then carry out the judgement that time step requires whether to reach, if reached, computing terminates, and extracts operation result and provides design parameter for rotating machinery design;

If do not reach set time step, then repeat step (1), then carry out the calculating of iterative loop, till the requirement time step of computing reaches setting value.

Correction fluid motion equation is the Scale invariant pattern type based on Scale invariant type mechanics theory and broad sense, proposes the relationship between adjacent yardstick: w _β=U _β+1, the Scale invariant type convection velocity w namely under yardstick β equals the local speed U of Scale invariant type under yardstick β+1; Simultaneously in conjunction with fluid convection and diffusion principle: U _β=w _β+ v _β, i.e. local speed U _βequal convection velocity w _βwith rate of propagation v _βsuperposition, derive the Boltzmann equation of Scale invariant type: and then draw Scale invariant type fluid motion equation or revise fluid motion equation to the expression-form of incompressible flow: $\frac{\∂}{\∂t}{U}_{\mathrm{\β}}+w_{\mathrm{\β}}\·\▿U_{\mathrm{\β}}-{\mathrm{\ν}}_{\mathrm{\β}}{\▿}^2{U}_{\mathrm{\β}}=-\frac{\▿p_{\mathrm{\β}}}{{\mathrm{\ρ}}_{\mathrm{\β}}}.$

Dynamic solution method based on correction fluid motion equation linear iteration is first started with from the angle of math equation, provides another kind of a _pand a _nmathematical method.These two coefficients constitute transmission speed, and independent of soluble variable local velocity, can find a _pand a _ndirectly can be fixed on the initial step of limited bulk computing method.Therefore, these two coefficients need not participate in the iterative calculating of PISO algorithm algorithm circulation.In this approach, it is linear for revising fluid motion equation, and that is a rate-pressure coupling difficult problem is solved by solving equation and iterative algorithm double optimization technology.

Figure 2 shows that centrifugal pump flow field stress and strain model schematic diagram, this figure for physical model, carries out stress and strain model to this physical model with a centrifugal water pump impeller interior flow field simultaneously.Under FVM method, when adopting update equation solver and NS equation solver to carry out numerical simulation calculation, all adopt the stress and strain model pattern of this physical model.

Fig. 3 (a) is depicted as NS solver, revises solver, experimental measurements at moving vane and stator blade coupling regime, the ratio curve comparison diagram of relative velocity vertical component and peripheral speed under its blade radial gap ratio condition; Fig. 3 (b) is NS solver, revises the ratio curve comparison diagram of solver, experimental measurements relative velocity horizontal component and peripheral speed under same radial gap ratio condition.In figure, the circumferential position of impeller moving vane and diffuser stator blade has been illustrated respectively in the top and bottom of form, to illustrate the circumferential distance of two kinds of situations relative to transverse axis.Understand the high consistency between numerical result and measurement result by contrast table, except peak-to-valley value region, two numerical results and experimental measurements curve substantially identical, this illustrate two kinds of numerical solution utensils have higher accuracy.Also can find out NS solver by comparison diagram and revise the matching that solver has height in numerical solution process, two kinds of numerical solution devices, under result of calculation accurately prerequisite, calculate effect basically identical.By the quantitative comparison of relative velocity component, we are understood the effect that stator-rotor interaction causes better: in relative velocity distribution, we can observe and organize peak-to-valley value more, and these values all produce at impeller moving vane trailing edge and diffuser stator blade leading edge; Between the flow field area of moving vane and stator blade, velocity gradient is relatively little.The more important thing is, diffuser vane also creates very large impact to the output stream of impeller, and cause flow velocity creating speed minimum value near its leading edge place, especially when directly impacting blade when the air-flow along blade, speed experienced by very large graded.

Figure 4 shows that the iteration time comparison diagram of NS equation solver and update equation solver of the present invention.Shown in figure, NS fluid motion equation solver iterative convergent process needs computing time ratio to be 2.3, and revising fluid motion equation solver iterative convergent process needs be 1.0 than then computing time, and computing time is saved over half.This is mainly because revise fluid motion equation algorithm to save coefficient a _pand a _nin the process of each time step double counting, the iterative mode of being programmed by this algorithmic code decreases the convergence oscillatory occurences that NS fluid motion equation solver obviously occurs in an iterative process.This instability concussion needs extra time step to calculate, thus adds computing time.This is also revise fluid motion equation solver to decrease total run time main cause.

From Fig. 3 and Fig. 4 in general, correction fluid motion equation dynamic mesh method of value solving of the present invention solves consistance in guarantee and experimental measurements and NS solver, under guaranteeing the prerequisite of accuracy, improve the efficiency of rotating machinery design value analog computation, compared with prior art there is obvious advantage.

The above is only the preferred embodiment of the present invention, and the present invention should not be confined to the content disclosed in this embodiment and accompanying drawing.The amendment of any unsubstantiality done within every basic ideas at algorithm proposed by the invention and framework, conversion and improvement all should be included within protection scope of the present invention.

Claims (2)

1., based on the dynamic mesh method of value solving revising fluid motion equation, it is characterized in that comprising the following steps:

(1) built the physical arrangement model of rotating machinery design by computing machine, complete dynamic remeshing, based on Scale invariant type statistical parsing, conditions setting convection velocity w, Grid Mobile speed U _b, kinematic viscosity ν, model cootrol body volume V, model cootrol bulk area s initiation parameter, according to rotating machinery designing requirement setting-up time step-length, carry out loop iteration calculating by setting-up time step-length:

According to dynamic flow correction fluid motion equation:

First absolutization process is carried out to flux φ, then dynamic mesh displacement x is determined, and then flux φ is revised, complete the relative processing of the relative Grid Mobile speed of flux φ; In formula: q _φfor the pressure p solved;

(2) spot speed term coefficient a to be solved is calculated _pwith solve neighbor point speed term coefficient a _n, wherein, $a_P=2v\frac{1}{\mathrm{\Δx}}+(w_P-U_{\mathrm{bP}}),a_N=v\frac{1}{\mathrm{\Δx}}+(w_N-U_{\mathrm{bN}}),$ W _pfor the convection velocity at P point place, U _bPfor the Grid Mobile speed at P point place, w _nfor the convection velocity of P point adjacent place each point, U _bNfor the Grid Mobile speed of adjacent place each point;

(3) iterative process of PISO algorithm circulation is carried out, first compute matrix H (U): u in formula _nrepresent the speed of P point adjacent place each point;

Then calculating pressure p: the equation of momentum at corresponding P point place is: rate equation is: $U_P=\frac{H\left(U\right)}{a_P}-\frac{\▿p}{a_P},$ Pressure equation is: $\▿\·(\frac{1}{a_P}\▿p)=\underset{f}{\mathrm{\Σ}}A\·{\left(\frac{H\left(U\right)}{a_P}\right)}_f,$ Pressure p is the q in equation _φ;

Further relative processing is carried out to flux φ;

To speed U _psolve, wherein, U _pds=φ;

(4) speed U is judged _pwhether restrain,

When PISO algorithm loop iteration evaluation tends towards stability, then judge speed U _pmeet the condition of convergence, then enter next step operation program;

When PISO algorithm loop iteration evaluation still changes instability, then judge speed U _pdo not meet the condition of convergence, then repeat step (3), then carry out the iterative computation of PISO algorithm circulation, until speed U _ptill meeting the condition of convergence;

(5) judge to calculate whether reach time step,

As speed U _pmeet the condition of convergence, then carry out the judgement that time step requires whether to reach, if reached, computing terminates, and extracts operation result and provides design parameter for rotating machinery design;

If do not reach set time step, then repeat step (1), then carry out the calculating of iterative loop, till the requirement time step of computing reaches setting value.

2. a kind of dynamic mesh method of value solving based on revising fluid motion equation according to claim 1, is characterized in that: described setting-up time step-length is determined according to the rotating speed actual demand of different rotary machinery.