CN105808954A - Periodic unsteady flow field prediction method suitable for CFD numerical simulation - Google Patents

Abstract

The invention relates to a periodic unsteady flow field prediction method suitable for CFD (computational fluid mechanics) numerical simulation. According to the method, Fourier transform is adopted to transform a time-dependent periodic unsteady problem in a time domain into N coupled steady equations, and the whole period can be reconstructed through coupling and computing N steady moments, so that the computing efficiency can be greatly improved under the condition of not reducing the computing accuracy. The method can be used for efficiently and accurately obtaining the periodic unsteady flow field distribution of aircrafts, obtaining the periodic unsteady aerodynamic force of the aircrafts according to the obtained flow field distribution and guiding the dynamic stability design of the aircrafts according to the obtained periodic unsteady aerodynamic force.

Description

A kind of Forecasting Methodology of the cycle Unsteady Flow suitable in CFD numerical simulation

Technical field

The present invention relates to the Forecasting Methodology of a kind of cycle Unsteady Flow being applicable to CFD (Fluid Mechanics Computation) numerical simulation, the method can obtain the nonstationary flow field distribution of aircraft cycle on high-efficiency high-accuracy ground, and the cycle unsteady aerodynamic force of aircraft can be obtained according to the Flow Field Distribution obtained, then the dynamic stability of aircraft can be instructed to design according to the cycle unsteady aerodynamic force obtained.

Background technology

Periodically UNSTEADY FLOW has engineering background widely, and the rotor of such as helicopter streams, and comprises the interior stream flowing of stator and the turbine of rotor, streaming of fan blade, and the research of flapping wing power set.In addition numerical value forced oscillation method calculates dynamic derivative, allows object do the forced vibration of given frequency, the sluggish aerodynamic force that key is directly related with frequency of vibration when being also obtain periodic movement.Above-mentioned all flowings have an obvious characteristic frequency, and macroscopically the periodicity of this kind of flowing is all the result of object periodic movement.At present, generally adopting this UNSTEADY FLOW of pseudo-time Method numerical simulation, but pseudo-time Method does not consider the cyclophysis of flowing, therefore computational efficiency is not high, and particularly with hypersonic flowing owing to the characteristic time sharply reduces, assessing the cost can be higher.For reducing amount of calculation, some scholars propose a series of simplification algorithm, but Unsteady flow computation it is crucial that obtain good balance between computational efficiency and computational accuracy, the method of high-accuracy high-resolution is generally all limited by the computing capability of current computer, and low precision arithmetic generally can not react some important physical messages.Therefore, develop a kind of algorithm that can quickly calculate periodically UNSTEADY FLOW, its efficiency and the requirement of precision two aspect will be considered.

Summary of the invention

The technology of the present invention solves problem: overcome the deficiencies in the prior art, it is provided that the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation.

The technical solution of the present invention is:

A kind of Forecasting Methodology of the cycle Unsteady Flow suitable in CFD numerical simulation, step is as follows:

(1) grid is calculated as needed in generating in computational fields；

(2) a calculating cycle is divided into etc. N number of sampled point of time interval, now, the conservation variable W of the n-th sampled point in N-S governing equation_nDerivative term to the timeDiscrete it is:

$\frac{\∂}{\∂t}{W}_n=\underset{j=0}{\overset{N-1}{\Σ}}D\[n,j\]W_j$

Wherein n=0,1,2 ..., N-1, j=0,1,2 ..., N-1, N is natural number, and D [n, j] represents the element that dimension is the line n of N N matrix D, jth row；

(3) grid that calculates obtained according to step (1) adopts rigidity Moving mesh method to obtain the calculating grid of each sampled point；

(4) incoming flow initial condition is set in fluid diagnosis, by the calculating grid in moment corresponding in each sampled point setting steps (3), the time-derivative item in fluid diagnosis is rewritten as in step (2)Form, calculate and obtain N number of sampling instant convergence flow field and aerodynamic force；

(5) aerodynamic force of N number of sampling instant that will obtain in step (4), obtains Fourier expansion coefficient:

${\hat{C}}_k\≈\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}{C}_n{e}^{-ik\mathrm{\ω}n\mathrm{\Δ}t}$

Wherein Δ t=T/N, ω=2 π/T, T is the time in a cycle, and k is harmonic number, k=0 ,+1 ,-1 ,+2 ,-2 ..., C_nIt is the aerodynamic force of the n-th sampled point,For Fourier space kth level number, aerodynamic force C during moment t in the reconstruct cycle_tFor:

$C_t\≈\underset{k=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{\hat{C}}_k{e}^{ik\mathrm{\ω}t}.$

In step (2), W_n=(ρ, ρ u, ρ v, ρ w, ρ E)_n ^T, wherein ρ is fluid density, and (u, v, w) for the velocity component under rectangular coordinate system, E is the gross energy of unit mass gas.

In step (2), the expression formula of D [n, j] is:

N is odd number

N is even number.

In step (3), rigidity Moving mesh method for carrying out integral-rotation or translation to initial mesh.

In step (4), incoming flow initial condition includes reynolds number Re, Mach number Ma and temperature of incoming flow T_∞。

Equal Navier-Stokes (RANS) governing equation solver when fluid diagnosis in step (4) is three dimensional compressible Reynolds, uses two kinds of turbulence model simulation turbulent flow of SA, MenterSST in solver.

Beneficial effect

(1) cycle Non-steady Problem time dependent in time domain is transformed to the permanent equation of N number of coupling by the time spectrum method of the present invention, only need to couple calculating N number of permanent moment, i.e. the restructural whole cycle, thus the raising computational efficiency of larger amplitude.To hypersonic flowing, the time spectrum method of the present invention is not when reducing computational accuracy, and computational efficiency can improve a more than magnitude；

(2) in the method for the present invention, rigidity Moving mesh method obtains the calculating grid of each sampled point, it is not necessary to change mesh topology.

(3) existing solver is had good inheritance by the time spectrum method of the present invention, wherein in solver flux computing module without amendment, it is only necessary to modification time derivative term.

(4) form that when the time spectrum method of the present invention is to cycle unsteady flo w mixing Reynolds, in equal Navier-Stokes equation (RANS) equation, the discrete employing of the time-derivative item of turbulence model equation is unified.

Accompanying drawing explanation

Fig. 1 is the initial calculation grid in embodiment 1；

The grid that Fig. 2 is the initial 0.00s moment rotates the calculating grid obtaining the 0.02s moment；

The grid that Fig. 3 is the initial 0.00s moment rotates the calculating grid obtaining the 0.04s moment；

The grid that Fig. 4 is the initial 0.00s moment rotates the calculating grid obtaining the 0.06s moment；

The grid that Fig. 5 is the initial 0.00s moment rotates the calculating grid obtaining the 0.08s moment；

Fig. 6 is the pitching moment coefficient retardant curve with angle of attack variation in NACA0015 aerofoil profile cycle, wherein, abscissa is the angle of attack, vertical coordinate is pitching moment coefficient, TSM is the retardant curve adopting the method for the present invention to obtain as N=5, the retardant curve that the pseudo-time Method that Dual-Time-Step is traditional obtains, Exp is the result of the test curve of retardant curve；

Fig. 7 is the method comparison with traditional pseudo-time Method CPU calculating time of the present invention；Wherein, TSM adopts the method CPU of the present invention time calculated, the time that the pseudo-time Method CPU that Dual-Time-Step is traditional calculates.

Detailed description of the invention

Set up the effective ways calculating UNSTEADY FLOW, it should take into full account its flow performance.Periodically UNSTEADY FLOW is characterized by that, after several intervals, the flow behavior in flow field can repeat.Traditionally, this cycle unsteady flo w is same with general Non-steady Problem treats.The algorithm major part of past exploitation is all adopt time stepping method.This is owing at time orientation, the solution in any moment all only affects the solution of future time instance.Most of UNSTEADY FLOW are all acceptables by the throwing materilization freatment method of this time orientation, and such as aircraft receives pulse-type disturbance and does maneuvering flight, and the initial condition at this moment determined must accurately, and transient process is also critically important.But in periodicity Non-steady Problem, the hypothesis that this time orientation only affects backward is no longer completely rationally, if because only considering the cycle of a convergence, and being left out influencing each other of different cycles, at this moment it is believed that in the cycle solution in any moment all affect the solution in other moment in this cycle.To this kind of flowing, being generally non-physical by transient process initial condition to steady statue solution, it is very slow that this transient process probably restrains, and substantial amounts of cpu resource all expends in this.To this cycle Non-steady Problem, time-marching method is all too expensive, it is clear that a kind of computational methods direct solution periodic state of exploitation, and avoids solving the transient process of non-physical, by bigger raising computational efficiency.

To time spectrum method, current all flux calculate and all carry out in time domain, and existing solver is had good inheritance.

1D matrix calculus

Adopting Fourier's change, the time-derivative in n moment can be expressed as D matrix and the product of sample sequence in the cycle.

2 sampled point grid computings

Couple the flow field of each sampled point in the calculating cycle, first to obtain the calculating grid of each sampled point.Rigidity Moving mesh method is adopted to obtain the calculating grid of each sampled point.During sampled point N=5, the calculating grid schematic diagram of each sampled point is as shown in the figure.

3 initialize flow field

Initialization for flow field has two ways: the primary data of (1) each sampled point is all given as inlet flow conditions；(2) calculating an initial steady flow field, the primary data of each sampled point is all given as this initial steady flow field.Through comparative study, two kinds of initialization modes are little on constringent impact.

4 flux calculate

Calculate without viscous flux and viscosity flux and all can inherit existing solver, but flux to add, in calculating, the time-derivative comprising D matrix.

5 time stepping method

Permanent virtual time one step process is adopted to carry out time stepping method.The time spectrum method adopting fully implicit solution carries out time stepping method, can improve convergence stability.4th step is repeatable repeatedly to the 5th step, until obtaining satisfied result.

6 reconstruct cycle flow fields

Flow field and aerodynamic force to calculated each sampling instant, utilizes the Fourier inversion can whole cycle flow field.

A kind of efficient period UNSTEADY FLOW Forecasting Methodology, adopt Fourier transformation that cycle Non-steady Problem time dependent in time domain is transformed to the permanent equation of N number of coupling, only need to couple calculating N number of permanent moment, i.e. the restructural whole cycle, thus the raising computational efficiency of larger amplitude.

To SA turbulence model, that MenterSST turbulence model equation also uses time spectrum method is discrete.

The time spectrum method adopting fully implicit solution carries out time stepping method.

Adopt MPI framework, it is achieved that parallel computation, improve calculating scale.

Invention describes the efficient Forecasting Methodology of a kind of cycle UNSTEADY FLOW.Adopt Fourier transformation that cycle Non-steady Problem time dependent in time domain is transformed to the permanent equation of N number of coupling, only need to couple calculating N number of permanent moment, the i.e. restructural whole cycle, when not reducing computational accuracy, can the raising computational efficiency of larger amplitude.

The method that the present invention introduces is compared with existing method, can be applicable to equal Navier-Stokes equation (RANS) equation during cycle unsteady flo w Reynolds, to SA (Spalart-Allmaras) turbulence model, that MenterSST (Shear-StressTransport) turbulence model also uses time spectrum method is discrete, improves simulation precision and Complex Flows adaptability.All flux calculate and all carry out in time domain, and existing solver is had good inheritance.Adopt the fully implicit solution form of time spectrum method, improve the stability problem increased along with sampling number and occur.Parallel computation can be realized, it is adaptable to large-scale engineer applied problem.

A kind of Forecasting Methodology of the cycle Unsteady Flow suitable in CFD numerical simulation, step is as follows:

(1) grid is calculated as needed in generating in computational fields；

(2) a calculating cycle is divided into etc. N number of sampled point of time interval, now, conservation variable W in N-S governing equation_nDerivative term to the timeDiscrete it is: (wherein N is natural number, n=0,1,2 ... N-1, j=0,1,2 ... N-1)

$\frac{\∂}{\∂t}{W}_n=\underset{j=0}{\overset{N-1}{\Σ}}D\[n,j\]W_j$

Wherein D [n, j] is the line n that dimension is N N matrix D, the element of jth row, and the expression formula of D [n, j] is:

N is odd number

N is even number

W=(ρ, ρ u, ρ v, ρ w, ρ E)^TFor conservation variable, wherein ρ is fluid density, and (u, v, w) for the velocity component under rectangular coordinate system, E is the gross energy of unit mass gas；Then W_n=(ρ, ρ u, ρ v, ρ w, ρ E)_n ^TRepresent the conservation variable of the n-th sampled point；

(3) grid that calculates obtained according to step (1) adopts rigidity Moving mesh method to obtain the calculating grid of each sampled point；

(4) incoming flow initial condition is set in fluid diagnosis, including reynolds number Re, Mach number Ma, and temperature of incoming flow T_∞, by the calculating grid in moment corresponding in each sampled point setting steps (3), the time-derivative item in fluid diagnosis is rewritten as in step (2)Form, calculates and obtains convergent current field and aerodynamic force；

(5) aerodynamic force of N number of sampling instant that will obtain in step (4), obtains Fourier expansion coefficient:

${\hat{C}}_k\≈\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}{C}_n{e}^{-ik\mathrm{\ω}n\mathrm{\Δ}t}$

Wherein Δ t=T/N, ω=2 π/T, T is the time in a cycle, and k is harmonic number, k=0 ,+1 ,-1 ,+2 ,-2 ...；C_nIt is the aerodynamic force of the n-th sampled point,For Fourier space kth level number, finally, the aerodynamic force of any time in the restructural cycle:

$C_t\≈\underset{k=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{\hat{C}}_k{e}^{ik\mathrm{\ω}t}.$

Embodiment

(1) computation model adopts NACA0015 aerofoil profile, and aerofoil profile chord length is c=0.3048m, is generate in computational fields to calculate grid within the scope of ten times of chord lengths, as it is shown in figure 1, airfoil surface has 260 points, wall ground floor mesh spacing is 3 × 10^-6M, NACA0015 aerofoil profile makes sinusoidal pitch vibration, namely

α=4.0+4.2sin (20 π t)

Wherein, α is the angle of attack of NACA0015 aerofoil profile, and t is time of vibration；

(2) a calculating cycle is divided into etc. 5 sampled points of time interval, and respectively (0.00s, 0.02s, 0.04s, 0.06s, 0.08s), now, D matrix is:

$D=\left[\begin{array}{ccccc}0& 53.45& -33.03& 33.03& -53.45\\ -53.45& 0& 53.45& -33.03& 33.03\\ 33.03& -53.45& 0& 53.45& -33.03\\ -33.03& 33.03& -53.45& 0& 53.45\\ 53.45& -33.03& 33.03& -53.45& 0\end{array}\right]$

(3) rigidity Moving mesh method is adopted to rotate the calculating grid respectively obtaining 0.02s, 0.04s, 0.06s, the 0.08s moment according to the grid in initial 0.00s moment, successively as shown in Figure 2-5；

(4) incoming flow initial condition, wherein reynolds number Re=1.95 × 10 are set in fluid diagnosis⁶, Mach number Ma=0.29, and temperature of incoming flow T_∞=288.16K；

(5) the time-derivative item in fluid diagnosis adopt D matrix in step (2) carry out discrete, namely

$\left[\begin{array}{c}\∂W_0/\∂t\\ \∂W_1/\∂t\\ \∂W_2/\∂t\\ \∂W_3/\∂t\\ \∂W_4/\∂t\end{array}\right]=\left[\begin{array}{ccccc}0& 53.45& -33.03& 33.03& -53.45\\ -53.45& 0& 53.45& -33.03& 33.03\\ 33.03& -53.45& 0& 53.45& -33.03\\ -33.03& 33.03& -53.45& 0& 53.45\\ 53.45& -33.03& 33.03& -53.45& 0\end{array}\right]\left[\begin{array}{c}{W}_0\\ W_1\\ W_2\\ W_3\\ W_4\end{array}\right]$

(6) the calculating grid of five sampling instants that fluid diagnosis obtains according to step (3), step (4) are arranged incoming flow initial condition and step (5) discrete after time-derivative item calculate the pitching moment coefficient obtaining each sampled point, for

$\left[\begin{array}{c}{\mathrm{Cmz}}_0\\ {\mathrm{Cmz}}_1\\ {\mathrm{Cmz}}_2\\ {\mathrm{Cmz}}_3\\ {\mathrm{Cmz}}_4\end{array}\right]=\left[\begin{array}{c}6.096E-04\\ -1.192E-02\\ -1.919E-02\\ -9.032E-03\\ 2.250E-03\end{array}\right]$

(7) pitching moment coefficient of N number of sampling instant that will obtain in step (6), obtains Fourier expansion coefficient:

${\hat{C}}_k\≈\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}{C}_n{e}^{-ik\mathrm{\ω}n\mathrm{\Δ}t}$

$\left[\begin{array}{c}{\hat{C}}_{-2}\\ {\hat{C}}_{-1}\\ {\hat{C}}_0\\ {\hat{C}}_1\\ {\hat{C}}_2\end{array}\right]=\left[\begin{array}{c}-5.77E-05+i2.66E-04\\ 4.09E-03-i3.89E-03\\ -7.46\\ 4.09E-03+i3.89E-03\\ -5.77E-05-i2.66E-04\end{array}\right]$

(8) according to the result of step (7) the pitching moment coefficient C obtaining any time in the cycle in conjunction with formula (2)_t, and in conjunction with angle of attack rule formula (3) over time, obtain the pitching moment coefficient retardant curve with angle of attack variation in NACA0015 aerofoil profile cycle, as shown in Figure 6；

$C_t\≈\underset{k=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{\hat{C}}_k{e}^{ik\mathrm{\ω}t}---\left(2\right)$

α=4.0+4.2sin (20 π t) (3)

According to Fig. 6 it can be seen that adopt the result that the method for the present invention obtains to coincide better with test value.

According to Fig. 7 it can be seen that adopt the CPU of the method for the present invention Time Calculation calculated to be about the 1/10 of the tradition dual time-stepping method CPU time calculated, computational efficiency can improve a more than magnitude；

The non-detailed description of the present invention is known to the skilled person technology.

Claims (6)

1. the Forecasting Methodology of the cycle Unsteady Flow being applicable to CFD numerical simulation, it is characterised in that step is as follows:

(1) grid is calculated as needed in generating in computational fields；

(2) a calculating cycle is divided into etc. N number of sampled point of time interval, now, the conservation variable W of the n-th sampled point in N-S governing equation_nDerivative term to the timeDiscrete it is:

$\frac{\∂}{\∂t}{W}_n=\underset{j=0}{\overset{N-1}{\Σ}}D\[n,j\]W_j$

Wherein n=0,1,2 ..., N-1, j=0,1,2 ..., N-1, N is natural number, and D [n, j] represents the element that dimension is the line n of N N matrix D, jth row；

(3) grid that calculates obtained according to step (1) adopts rigidity Moving mesh method to obtain the calculating grid of each sampled point；

(4) incoming flow initial condition is set in fluid diagnosis, by the calculating grid in moment corresponding in each sampled point setting steps (3), the time-derivative item in fluid diagnosis is rewritten as in step (2)Form, calculate and obtain N number of sampling instant convergence flow field and aerodynamic force；

(5) aerodynamic force of N number of sampling instant that will obtain in step (4), obtains Fourier expansion coefficient:

${\hat{C}}_k\≈\frac{1}{N}\underset{n=0}{\overset{N-1}{\Σ}}{C}_n{e}^{-ik\mathrm{\ω}n\mathrm{\Δ}t}$

Wherein Δ t=T/N, ω=2 π/T, T is the time in a cycle, and k is harmonic number, k=0 ,+1 ,-1 ,+2 ,-2 ..., C_nIt is the aerodynamic force of the n-th sampled point,For Fourier space kth level number, aerodynamic force C during moment t in the reconstruct cycle_tFor:

$C_t\≈\underset{k=-(N-1)/2}{\overset{(N-1)/2}{\Σ}}{\hat{C}}_k{e}^{ik\mathrm{\ω}t}.$

2. the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation according to claim 1, it is characterised in that: in step (2), W_n=(ρ, ρ u, ρ v, ρ w, ρ E)_n ^T, wherein ρ is fluid density, and (u, v, w) for the velocity component under rectangular coordinate system, E is the gross energy of unit mass gas.

3. the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation according to claim 1, it is characterised in that: in step (2), the expression formula of D [n, j] is:

4. the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation according to claim 1, it is characterised in that: in step (3), rigidity Moving mesh method for carrying out integral-rotation or translation to initial mesh.

5. the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation according to claim 1, it is characterised in that: in step (4), incoming flow initial condition includes reynolds number Re, Mach number Ma and temperature of incoming flow T_∞。

6. the Forecasting Methodology of a kind of cycle Unsteady Flow suitable in CFD numerical simulation according to claim 1, it is characterized in that: equal Navier-Stokes (RANS) governing equation solver when the fluid diagnosis in step (4) is three dimensional compressible Reynolds, solver uses two kinds of turbulence model simulation turbulent flow of SA, MenterSST.