CN104732109A - Trailing vortex prediction method for vertical shaft rotating power generation device - Google Patents

Abstract

The invention provides a trailing vortex prediction method for a vertical shaft rotating power generation device. The method comprises the concrete steps that firstly, airfoil contour coordinates etai and epsiloni are input, and potential flow solution and viscous boundary layer solution are performed; secondly, trailing vortex shedding is performed, wherein according to the circulation theorem, the position of a trailing vortex is determined, the strength omegai of the trailing vortex is solved, and airfoil stress and torque are solved; thirdly, after the airfoil stress and the torque are obtained, the next time step movement of an airfoil is determined, and a trailing vortex graph is drawn according to relevant data. The trailing vortex prediction method for the vertical shaft rotating power generation device is combined with a new Kutta condition on the basis of a vortex panel method, a limited vortex method is developed, and the trailing vortex prediction method has better adaptability to the dynamic stall phenomenon.

Description

A kind of trailing vortex forecasting procedure of vertical pivot rotary generating device

Technical field

The present invention relates to trailing vortex forecasting procedure field, particularly relate to a kind of trailing vortex forecasting procedure of vertical pivot rotary generating device.

Background technology

The demand of the whole world to the energy increases day by day.Wind energy represents as regenerative resource, has become the important component part of many National Sustainable Development Strategies, by 2011 the end of the year China blower fan installation amount accumulative reach 62GW, occupy the first in the world.But expansion blindly result in over capacity, according to the recent statistics data display of global investing bank, within 2011, China Power Grids can only receive 15GW-16GW, has the wind-powered electricity generation capacity of 20GW cannot be digested at present, and therefore Chinese feature electric industry needs to think deeply expansion plans further.A kind of strategic choice is: the resource and the regional superiority that play coastal regenerative resource, extensive development regenerable marine energy, and this energy development also having fully demonstrated economical rationality is seen.Therefore, blower fan enters shallow sea from land: wind energy and marine tidal-current energy core conversion equipment are impellers, and be mainly divided into transverse axis and vertical pivot two kinds according to version, wherein transverse axis form is business main flow.Experience and research show: be 1/3 of transverse axis turbine cost with installed capacity vertical turbine cost.It is cheap with it, easy to maintenance, flow to self-adaptation to various, is applicable to the main source being suitable as coastal economy group particularly island energy resource supply from the feature of net generating.

The core apparatus of vertical turbine is impeller, and the capacity usage ratio of flow dynamic characteristic to turbine of impeller has vital effect.At present, the hydrodynamic force forecasting procedure of impeller mainly contains three kinds: the stream tube model based on blade momentum theory, the vortex method based on potential vortex theory and based on the viscosity CFD method solving N-S equation.

Stream tube model can provide the impeller bulk fluid characteristic such as average load of the capacity usage ratio-speed ratio characteristic of impeller, power-flow velocity (rotating speed) characteristic, blade, can also provide the details such as braking surface flow field velocity flow profile.But, because stream tube model is on average supposed based on infinite blade, thus cannot the unsteady characteristic of accurate forecast blade and impulsive load; Secondly, the equation of momentum solves under large speed ratio, high solidity easily to be dispersed, and is therefore not too applicable to fluid behaviour when calculating turbine heavy duty.Vortex method can describe Flow details effectively, but also there is intrinsic defect: one is set up under potential barrier framework based on the method for vortex theory, therefore when turbine is operated in little speed ratio, attack angle of blade amplitude is larger, easily there is separation flow and dynamic stall phenomenon, viscous effects effect is comparatively large, calculates and easily disperses; Another is that vortex method is consuming time longer, can not meet the requirement of engineering Fast Prediction impeller fluid property.Viscosity CFD method calculates then consuming time, and result of calculation requires that software user has rich experience, uses when being suitable as the checking after turbine type selecting and Fine design.

Because the energy conversion device of vertical pivot blower fan and tidal current energy water turbine is impeller, the fluid power component of impeller normally one group of cross section is the blade of aerofoil profile, if have ignored the factors such as the interference of free surface effect, pylon effect, impeller and supporting construction, single hydrodynamic characteristics and mechanism analyzing vertical turbine from the angle of desirable impeller, it inherently can be summed up as multiple-blade circumferential motion problem.The Kutta's condition that blade is traditional when doing large attack angle high frequency motion is also inapplicable, needs in conjunction with suitable unsteady flo w Kutta's condition.

Summary of the invention

For above problem, the invention provides a kind of trailing vortex forecasting procedure of vertical pivot rotary generating device, it is on the basis of whirlpool panel method, in conjunction with a kind of new Kutta's condition, developed limited vortex method, the method has better adaptability to dynamic stall phenomenon, for reaching this object, the invention provides a kind of trailing vortex forecasting procedure of vertical pivot rotary generating device, concrete steps are as follows: 1) input aerofoil profile week line coordinates η _i, ε _i, carry out potential barrier and solve and solve with viscous boundary layer; Described potential barrier solve first arrange source converge with the linear whirlpool of mean camber line, rear in conjunction with Kutta's condition and normal direction can not penetrate equation of condition group obtain source converge intensity σ _i, γ;

Described system of equations is as follows:

Set up earth coordinates (O:XY), wing local coordinate is (o:xy).The angular velocity Ω that the motion of wing can be rotated around reference point o by the translational velocity U of local coordinate reference point o and wing describes; If incoming flow is V _a(t), at infinity speed is:

V _A＝(V _x,V _y) (1)

Fluid domain τ _ethe velocity potential Φ (p, t) of interior field point p meets Laplace equation:

${\▿}^2\mathrm{\Φ}=0---\left(2\right)$

At aerofoil surface S _bon meet normal direction and can not penetrate condition:

$\frac{\∂\mathrm{\Φ}}{\∂n}{|}_{S_b}=(U+\mathrm{\Ω}\×r)\·n_b---\left(3\right)$

Wherein r is the radius vector of point under local coordinate system (o:x, y) in aerofoil surface; n _bunit normal vector for this directed towards object inside:

Velocity potential at infinity meets: Φ (p, t) ≈ xV _x+ yV _yp → ∞ (4)

If perturbation velocity potential is then meet following Definite problem:

Governing equation:

Normal direction can not penetrate condition:

Distant place condition:

Starting condition:

At wing profile mean camber line S _mdistribute linear whirlpool to simulate lift, mean camber line coordinate starting point is based upon the first edge of wing and leads round circle centre position, if mean camber line arc length is C _f, line taking vortices breakdown distribution function is:

γ (s)=2 γ _f(C _f-s)/C _f(9) the total intensity Γ in mean camber line distribution whirlpool _ffor:

${\mathrm{\Γ}}_f=\underset{S_m}{\∫}\mathrm{\γ}\left(s\right)\mathrm{ds}={\mathrm{\γ}}_f{C}_f---\left(10\right)$

At aerofoil surface S _bupper distribution constant value source in order to simulated thickness effect, from lower surface trailing edge place along clockwise direction by discrete for aerofoil surface be N number of unit, be designated as S _b1..., S _bNif the constant value source strength of i-th unit is σ _i; Tailwater system S _wreplace with discrete point unit vortex is approximate, suppose that each time step is de-and let out single trailing vortex, the newborn trailing vortex vortex strength of note kth time step is the velocity potential of any field point p of such singular point in flow field is:

To field point p differentiate, obtain induced velocity:

$V(p,t)={\▿}_p\underset{S_b}{\∫}\mathrm{\σ}(q,t)G(p,q){\mathrm{ds}}_q+\underset{S_m+S_w}{\∫}\mathrm{\Γ}(q,t)K(p,q){\mathrm{ds}}_q---\left(12\right)$

Formula (11), in formula (12), σ and Γ represents distribution singular point source, vortices breakdown respectively, and all the other are expressed as follows:

$G(p,q)=\frac{1}{2\mathrm{\π}}1n{r}_{p,q}---\left(13\right)$

$\mathrm{\θ}(p,q)=\frac{1}{2\mathrm{\π}}\arctan \frac{y_p-y_q}{x_p-x_q}---\left(14\right)$

$K(p,q)=\frac{1}{2\mathrm{\π}}(\frac{y_p-y_q}{r_{\mathrm{pq}}^2},-\frac{x_p-x_q}{r_{\mathrm{pq}}^2})---\left(15\right)$

$r_{p,q}=\sqrt{{(x_p-x_q)}^2+{(y_p-y_q)}^2}---\left(16\right)$

If the total intensity of trailing vortex is Γ _w, get time series:

t _k＝kΔt(k＝1,2,…) (17)

T _kmoment velocity potential the Algebraic Equation set meeting Definite problem Chinese style (6) at aerofoil surface reference mark place is:

$\begin{array}{c}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_{i,j}{\mathrm{\σ}}_j^{\left(k\right)}+A_{i,N+1}\·{\mathrm{\γ}}_f^{\left(k\right)}+\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{C}_{i,j}\·{\mathrm{\γ}}_w^{\left(j\right)}\\ =n_i\·(U-V_A+\mathrm{\Ω}\×r_i),(i=\mathrm{1,2},...,N)\end{array}---\left(18\right)$

Wherein:

the intensity of-jth time step trailing vortex;

R _ithe radius vector of i-th unit controls point under wing coordinate on-object plane;

N _ithe unit normal vector at i-th unit controls point place on-object plane;

Influence coefficient is expressed as follows:

$A_{i,j}=\underset{S_{\mathrm{bj}}}{\∫}\frac{\∂}{\∂n_i}G(p_i,q){\mathrm{ds}}_q,j=\mathrm{1,2},...,N---\left(19\right)$

$A_{i,N+1}=n_i\·\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{S_{\mathrm{mj}}}{\∫}2(C_f-s)/C_f\·K(p_i,q){\mathrm{ds}}_q,i=\mathrm{1,2},...,N---\left(20\right)$

C _i,j＝n _i·Κ(p _i,q _w,j)i＝1,2,…,N；j＝1,2,…,k (21)

System of equations (18) has N number of equation, and solution is current step distribution singular point intensity (σ ₁..., σ _n, γ _f, γ _w) ^twith implicit newborn vortex location n+3 unknown quantity altogether, so need supplementary 3 equations to make full scale equation group close;

2) trailing vortex takes off bryonia and solves trailing vortex vortex strength ω according to circulation theorem determination vortex location _i, and solve the stressed and moment of aerofoil profile;

3), after solving the stressed and moment of aerofoil profile, determine the motion of aerofoil profile future time step, draw trailing vortex figure according to related data.

Further improvement of the present invention, described step 1 Kutta's condition is as follows:

p _u-p _d＝Δp (22)

Kelvin theorem

$\frac{d{\mathrm{\Γ}}_{\mathrm{Total}}}{\mathrm{dt}}=\frac{d({\mathrm{\Γ}}_f+{\mathrm{\Γ}}_w)}{\mathrm{dt}}=0---\left(23\right)$

The position of newborn trailing vortex determine according to following empirical formula ^[9]

$r_{w,k}^{\left(k\right)}={\stackrel{\‾}{r}}_{\mathrm{TE}}-\mathrm{\β}(U-V_A+\mathrm{\Ω}\×{\stackrel{\‾}{r}}_{\mathrm{TE}}-{\stackrel{\‾}{V}}_{\mathrm{TE}})\mathrm{\Δt}---\left(24\right)$

Wherein β is a constant, gets 0.4 ~ 0.6.

The present invention is by vertical pivot impeller fluid property forecasting problem in developing for wind energy and marine tidal-current energy, propose a kind of new computing method, by the limited vortex method that whirlpool panel method combines with a kind of new Kutta's condition, the example that steady motion and multiple-blade fly steady motion is done by individual blade, confirm the validity of the method, it has developed limited vortex method, and the method has better adaptability to dynamic stall phenomenon.

Accompanying drawing explanation

Fig. 1 is algorithm flow chart of the present invention;

Fig. 2 is reference frame figure of the present invention;

Fig. 3 is singular point distribution mode schematic diagram of the present invention;

Fig. 4 is the blade surface pressure distribution under the different angle of attack of the embodiment of the present invention limited whirlpool model NACA0018;

Fig. 5 is that embodiment of the present invention Fluent calculates NACA0012 different angle of attack lower blade surface pressure distribution;

Fig. 6 is the vortex street that embodiment of the present invention NACA0012 aerofoil profile does heave movement;

Fig. 7 is that blade normal force coefficient of the present invention is with position angle change (Turbine A, λ=2.5);

Fig. 8 is that blade tangential force coefficient of the present invention is with position angle change (Turbine A, λ=2.5);

Embodiment

Below in conjunction with drawings and Examples, invention is described in detail:

The invention provides a kind of trailing vortex forecasting procedure of vertical pivot rotary generating device, it is on the basis of whirlpool panel method, and in conjunction with a kind of new Kutta's condition, developed limited vortex method, the method has better adaptability to dynamic stall phenomenon.

As an embodiment of the present invention, the invention provides the trailing vortex forecasting procedure of a kind of vertical pivot rotary generating device as shown in Figure 1, concrete steps are as follows:

1) aerofoil profile week line coordinates η is inputted _i, ε _i, carry out potential barrier and solve and solve with viscous boundary layer;

Described potential barrier solve first arrange source converge with the linear whirlpool of mean camber line, rear in conjunction with Kutta's condition and normal direction can not penetrate equation of condition group obtain source converge intensity σ _i, γ;

Described system of equations is as follows:

Set up earth coordinates (O:XY) as shown in Fig. 2 Fig. 3, wing local coordinate is (o:xy).The angular velocity Ω that the motion of wing can be rotated around reference point o by the translational velocity U of local coordinate reference point o and wing describes; If incoming flow is V _a(t), at infinity speed is:

V _A＝(V _x,V _y) (1)

Fluid domain τ _ethe velocity potential Φ (p, t) of interior field point p meets Laplace equation:

${\▿}^2\mathrm{\Φ}=0---\left(2\right)$

At aerofoil surface S _bon meet normal direction and can not penetrate condition:

$\frac{\∂\mathrm{\Φ}}{\∂n}{|}_{S_b}=(U+\mathrm{\Ω}\×r)\·n_b---\left(3\right)$

Wherein r is the radius vector of point under local coordinate system (o:x, y) in aerofoil surface; n _bunit normal vector for this directed towards object inside:

Velocity potential at infinity meets: Φ (p, t) ≈ xV _x+ yV _yp → ∞ (4)

If perturbation velocity potential is then meet following Definite problem:

Governing equation:

Normal direction can not penetrate condition:

Distant place condition:

Starting condition:

At wing profile mean camber line S _mdistribute linear whirlpool to simulate lift, mean camber line coordinate starting point is based upon the first edge of wing and leads round circle centre position, if mean camber line arc length is C _f, line taking vortices breakdown distribution function is:

γ (s)=2 γ _f(C _f-s)/C _f(9) the total intensity Γ in mean camber line distribution whirlpool _ffor:

${\mathrm{\Γ}}_f=\underset{S_m}{\∫}\mathrm{\γ}\left(s\right)\mathrm{ds}={\mathrm{\γ}}_f{C}_f---\left(10\right)$

At aerofoil surface S _bupper distribution constant value source in order to simulated thickness effect, from lower surface trailing edge place along clockwise direction by discrete for aerofoil surface be N number of unit, be designated as S _b1..., S _bNif the constant value source strength of i-th unit is σ _i; Tailwater system S _wreplace with discrete point unit vortex is approximate, suppose that each time step is de-and let out single trailing vortex, the newborn trailing vortex vortex strength of note kth time step is the velocity potential of any field point p of such singular point in flow field is:

To field point p differentiate, obtain induced velocity:

$V(p,t)={\▿}_p\underset{S_b}{\∫}\mathrm{\σ}(q,t)G(p,q){\mathrm{ds}}_q+\underset{S_m+S_w}{\∫}\mathrm{\Γ}(q,t)K(p,q){\mathrm{ds}}_q---\left(12\right)$

Formula (11), in formula (12), σ and Γ represents distribution singular point source, vortices breakdown respectively, and all the other are expressed as follows:

$G(p,q)=\frac{1}{2\mathrm{\π}}1n{r}_{p,q}---\left(13\right)$

$\mathrm{\θ}(p,q)=\frac{1}{2\mathrm{\π}}\arctan \frac{y_p-y_q}{x_p-x_q}---\left(14\right)$

$K(p,q)=\frac{1}{2\mathrm{\π}}(\frac{y_p-y_q}{r_{\mathrm{pq}}^2},-\frac{x_p-x_q}{r_{\mathrm{pq}}^2})---\left(15\right)$

$r_{p,q}=\sqrt{{(x_p-x_q)}^2+{(y_p-y_q)}^2}---\left(16\right)$

If the total intensity of trailing vortex is Γ _w, get time series:

t _k＝kΔt(k＝1,2,…) (17)

T _kmoment velocity potential the Algebraic Equation set meeting Definite problem Chinese style (6) at aerofoil surface reference mark place is:

$\begin{array}{c}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_{i,j}{\mathrm{\σ}}_j^{\left(k\right)}+A_{i,N+1}\·{\mathrm{\γ}}_f^{\left(k\right)}+\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{C}_{i,j}\·{\mathrm{\γ}}_w^{\left(j\right)}\\ =n_i\·(U-V_A+\mathrm{\Ω}\×r_i),(i=\mathrm{1,2},...,N)\end{array}---\left(18\right)$

Wherein:

the intensity of-jth time step trailing vortex;

R _ithe radius vector of i-th unit controls point under wing coordinate on-object plane;

N _ithe unit normal vector at i-th unit controls point place on-object plane;

Influence coefficient is expressed as follows:

$A_{i,j}=\underset{S_{\mathrm{bj}}}{\∫}\frac{\∂}{\∂n_i}G(p_i,q){\mathrm{ds}}_q,j=\mathrm{1,2},...,N---\left(19\right)$

$A_{i,N+1}=n_i\·\underset{j=1}{\overset{M}{\mathrm{\Σ}}}\underset{S_{\mathrm{mj}}}{\∫}2(C_f-s)/C_f\·K(p_i,q){\mathrm{ds}}_q,i=\mathrm{1,2},...,N---\left(20\right)$

C _i,j＝n _i·Κ(p _i,q _w,j)i＝1,2,…,N；j＝1,2,…,k (21)

System of equations (18) has N number of equation, and solution is current step distribution singular point intensity (σ ₁..., σ _n, γ _f, γ _w) ^twith implicit newborn vortex location n+3 unknown quantity altogether, so need supplementary 3 equations to make full scale equation group close;

2) trailing vortex takes off bryonia and solves trailing vortex vortex strength ω according to circulation theorem determination vortex location _i, and solve the stressed and moment of aerofoil profile;

3), after solving the stressed and moment of aerofoil profile, determine the motion of aerofoil profile future time step, draw trailing vortex figure according to related data.

Wherein Kutta's condition is as follows:

p _u-p _d＝Δp (22)

Kelvin theorem

$\frac{d{\mathrm{\Γ}}_{\mathrm{Total}}}{\mathrm{dt}}=\frac{d({\mathrm{\Γ}}_f+{\mathrm{\Γ}}_w)}{\mathrm{dt}}=0---\left(23\right)$

The position of newborn trailing vortex determine according to following empirical formula ^[9]

$r_{w,k}^{\left(k\right)}={\stackrel{\‾}{r}}_{\mathrm{TE}}-\mathrm{\β}(U-V_A+\mathrm{\Ω}\×{\stackrel{\‾}{r}}_{\mathrm{TE}}-{\stackrel{\‾}{V}}_{\mathrm{TE}})\mathrm{\Δt}---\left(24\right)$

Wherein β is a constant, gets 0.4 ~ 0.6.

Concrete example 1 is as follows: the static lift coefficient of the different angle of attack of NACA0018 airfoil fan, in order to compare trailing edge upper and lower surface pressure differential equivalent parameters l to calculating the impact brought.

Parameter l is on the impact of lift coefficient for table 1 trailing edge pressure differential, NACA0018, α=5 °, 13 °, 30 °;

Lift coefficient increases with angle of attack approximately linear before corresponding NACA0018 airfoil stall respectively from small to large, lift coefficient reaches maximum to static stall when occurring, lift coefficient recovers three kinds of physical context gradually after stall to calculate the angle of attack chosen.Cls represents the static experiment value of lift coefficient, and Cl-Cal represents calculated value.As can be seen from Table 1, make the l < 0 that calculated value and experiment value are comparatively close, this means i.e. trailing edge place p _u> p _d, along with the angle of attack increases, making lift coefficient close to experiment value | l| is also larger.Fig. 4 gives the limited whirlpool the model calculation of aerofoil surface pressure distribution when angle of attack is respectively 5 ° and 13 °, significantly can find out the pressure differential of trailing edge.The result of calculation of some viscosity CFD technology also captures the existence of blade trailing edge upper and lower surface pressure differential, and Fig. 5 is the NACA0012 airfoil fan adopting Fluent to simulate is 6 ° and 12 ° at the angle of attack, and Reynolds number is 3 × 10 ⁶time blade surface pressure distribution.P when conclusion is Low Angle Of Attack _u> p _d, but during along with angle of attack increase stall generation, p _u< p _d.

Concrete example 2 is as follows: NACA0012 airfoil fan does different frequency and compares with experiment value with separated vorticcs track display, lift coefficient with the heave movement trailing vortex of amplitude, angle of attack=0 °, blade chord length C=1.0m, heave displacement y=h sin (ω t), k is reduced frequency (reduced frequency, ω C/U), h is the zero dimension amplitude relative to chord length.

Calculating the parameter adopted is ^[10]

Fig. 6 (a) k=3.0, h=0.20, St=kh=0.6

Fig. 6 (b) k=10.1, h=0.2, St=2.2

When 6 (a) figure Leaf does intermediate frequency heave movement, thrust can be obtained; When St number is larger, wing motion frequency is higher, because the asymmetry of trailing vortex can produce skew, occur two whirlpool pairing phenomenon (Dual-Mode), produce uplift, this phenomenon is called Knoller-Betz effect simultaneously.This example shows this structure adopting separated vorticcs model can catch trailing vortex preferably, and this double vortex structure is formed by trailing vortex and separated vorticcs acting in conjunction in essence, have employed higher-order spectrum method and has simulated this phenomenon.

Concrete example 3 is as follows: the instantaneous stressed comparison and analysis of Strickland A turbine;

People's vertical pivot blower fan experiments such as nineteen eighty-three Stickland, conveniently measurement and display flow field, this series of experiments room completes in pond.The data measured and record include the instantaneous tangential force of blade, instantaneous normal force, blade surface pressure distribution, flow field etc.Experiment maintains rotor speed and fixes, and changes speed of incoming flow, measures three speed ratios respectively as shown in Figure 7 and Figure 8.

As can be seen from analog result, this method is to the forecast of blade normal force coefficient peak value, valley a little less than experiment value, and in upstream, card result of calculation is better than downstream card, and this is because the interference between trailing vortex, blade causes.

The above is only preferred embodiment of the present invention, is not restriction the present invention being made to any other form, and any amendment done according to technical spirit of the present invention or equivalent variations, still belong to the present invention's scope required for protection.

Claims (2)

1. a trailing vortex forecasting procedure for vertical pivot rotary generating device, concrete steps are as follows, it is characterized in that:

1) aerofoil profile week line coordinates η is inputted _i, ε _i, carry out potential barrier and solve and solve with viscous boundary layer;

Described potential barrier solve first arrange source converge with the linear whirlpool of mean camber line, rear in conjunction with Kutta's condition and normal direction can not penetrate equation of condition group obtain source converge intensity σ _i, γ;

Described system of equations is as follows:

Set up earth coordinates (O:XY), wing local coordinate is (o:xy), and the angular velocity Ω that the motion of wing can be rotated around reference point o by the translational velocity U of local coordinate reference point o and wing describes; If incoming flow is V _a(t), at infinity speed is:

V _A＝(V _x,V _y) (1)

Fluid domain τ _ethe velocity potential Φ (p, t) of interior field point p meets Laplace equation:

▽ ²Φ＝0 (2)

At aerofoil surface S _bon meet normal direction and can not penetrate condition:

$\frac{\&PartialD;\mathrm{\&Phi;}}{\&PartialD;n}{|}_{S_b}=(U+\mathrm{\&Omega;}\&times;r)\&CenterDot;n_b---\left(3\right)$

Wherein r is the radius vector of point under local coordinate system (o:x, y) in aerofoil surface; n _bunit normal vector for this directed towards object inside:

Velocity potential at infinity meets: Φ (p, t) ≈ xV _x+ yV _yp → ∞ (4)

If perturbation velocity potential is then meet following Definite problem:

Governing equation:

Normal direction can not penetrate condition:

Distant place condition:

Starting condition:

At wing profile mean camber line S _mdistribute linear whirlpool to simulate lift, mean camber line coordinate starting point is based upon the first edge of wing and leads round circle centre position, if mean camber line arc length is C _f, line taking vortices breakdown distribution function is:

γ(s)＝2γ _f(C _f-s)/C _f(9)

The total intensity Γ in mean camber line distribution whirlpool _ffor:

${\mathrm{\&Gamma;}}_f=\underset{S_m}{\&Integral;}\mathrm{\&gamma;}\left(s\right)\mathrm{ds}={\mathrm{\&gamma;}}_f{C}_f----\left(10\right)$

At aerofoil surface S _bupper distribution constant value source in order to simulated thickness effect, from lower surface trailing edge place along clockwise direction by discrete for aerofoil surface be N number of unit, be designated as S _b1..., S _bNif the constant value source strength of i-th unit is σ _i; Tailwater system S _wreplace with discrete point unit vortex is approximate, suppose that each time step is de-and let out single trailing vortex, the newborn trailing vortex vortex strength of note kth time step is the velocity potential of any field point p of such singular point in flow field is:

To field point p differentiate, obtain induced velocity:

$V(p,t)={\&dtri;}_p\underset{S_b}{\&Integral;}\mathrm{\&sigma;}(q,t)G(p,q){\mathrm{ds}}_q+\underset{S_m+S_w}{\&Integral;}\mathrm{\&Gamma;}(q,t)K(p,q){\mathrm{ds}}_q---\left(12\right)$

Formula (11), in formula (12), σ and Γ represents distribution singular point source, vortices breakdown respectively, and all the other are expressed as follows:

$G(p,q)=\frac{1}{2\mathrm{\&pi;}}\ln r_{p,q}---\left(13\right)$

$\mathrm{\&theta;}(p,q)=\frac{1}{2\mathrm{\&pi;}}\arctan \frac{y_p-y_q}{x_p-x_q}---\left(14\right)$

$K(p,q)=\frac{1}{2\mathrm{\&pi;}}(\frac{y_p-y_q}{r_{\mathrm{pq}}^2},\frac{x_p-x_q}{r_{\mathrm{pq}}^2})---\left(15\right)$

$r_{p,q}=\sqrt{{(x_p-x_q)}^2+{(y_p-y_q)}^2}---\left(16\right)$

If the total intensity of trailing vortex is Γ _w, get time series:

t _k＝kΔt (k＝1,2,…) (17)

T _kmoment velocity potential the Algebraic Equation set meeting Definite problem Chinese style (6) at aerofoil surface reference mark place is:

$\begin{array}{c}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{A}_{i,j}{\mathrm{\&sigma;}}_j^{\left(k\right)}+A_{i,N+1}\&CenterDot;{\mathrm{\&gamma;}}_f^{\left(k\right)}+\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{C}_{i,j}\&CenterDot;{\mathrm{\&gamma;}}_w^{\left(j\right)}\\ =n_i\&CenterDot;(U-V_A+\mathrm{\&Omega;}\&times;r_i),(i=\mathrm{1,2},...,N)\end{array}---\left(18\right)$

Wherein:

the intensity of-jth time step trailing vortex;

R _ithe radius vector of i-th unit controls point under wing coordinate on-object plane;

N _ithe unit normal vector at i-th unit controls point place on-object plane;

Influence coefficient is expressed as follows:

$A_{i,j}=\underset{S_{\mathrm{bj}}}{\&Integral;}\frac{\&PartialD;}{{\&PartialD;n}_i}G(p_i,q){\mathrm{ds}}_q,j=1,2,...,N---\left(19\right)$

$A_{i,N+1}=n_i\&CenterDot;\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{S_{\mathrm{mj}}}{\&Integral;}2(C_f-s)/C_f\&CenterDot;K(p_i,q){\mathrm{ds}}_q,i=\mathrm{1,2},...,N---\left(20\right)$

C _i,j＝n _i·Κ(p _i,q _w,j) i＝1,2,…,N；j＝1,2,…,k (21)

System of equations (18) has N number of equation, and solution is current step distribution singular point intensity (σ ₁..., σ _n, γ _f, γ _w) ^twith implicit newborn vortex location n+3 unknown quantity altogether, so need supplementary 3 equations to make full scale equation group close;

2) trailing vortex takes off bryonia and solves trailing vortex vortex strength ω according to circulation theorem determination vortex location _i, and solve the stressed and moment of aerofoil profile;

3), after solving the stressed and moment of aerofoil profile, determine the motion of aerofoil profile future time step, draw trailing vortex figure according to related data.

2. the trailing vortex forecasting procedure of a kind of vertical pivot rotary generating device according to claim 1, is characterized in that: described step 1 Kutta's condition is as follows:

p _u-p _d＝Δp (22)

Kelvin theorem

$\frac{{\mathrm{d\&Gamma;}}_{\mathrm{Total}}}{\mathrm{dt}}=\frac{d({\mathrm{\&Gamma;}}_f+{\mathrm{\&Gamma;}}_w)}{\mathrm{dt}}=0---\left(23\right)$

The position of newborn trailing vortex determine according to following empirical formula ^[9]

$r_{w,k}^{\left(k\right)}={\stackrel{\&OverBar;}{r}}_{\mathrm{TE}}-\mathrm{\&beta;}(U-V_A+\mathrm{\&Omega;}\&times;{\stackrel{\&OverBar;}{r}}_{\mathrm{TE}}-{\stackrel{\&OverBar;}{V}}_{\mathrm{TE}})\mathrm{\&Delta;t}---\left(24\right)$

Wherein β is a constant, gets 0.4 ~ 0.6.