METHOD AND SYSTEM FOR DETERMINING HELICOPTER ROTOR AIRFOIL

Information

  • Patent Application
  • 20220033062
  • Publication Number
    20220033062
  • Date Filed
    August 07, 2020
    4 years ago
  • Date Published
    February 03, 2022
    2 years ago
Abstract
The present disclosure provides a method and system for determining a helicopter rotor airfoil. The method includes: randomly generating a sample point by using a Latin hypercube sampling (LHS) method (S1); determining characterization equations of upper and lower airfoil surfaces of an airfoil based on the airfoil sample point by using a class shape transformation (CST) method (S2); performing dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a computational fluid dynamics (CFD) method, to obtain a flow field characteristic of the airfoil (S3); establishing a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and training the mapping relationship by using a maximum likelihood estimation method and an expected improvement (EI) criterion, to obtain a trained mapping relationship (S4); determining an optimal sample point based on the trained mapping relationship by using Non-dominated Sorting Genetic Algorithm II (NSGA-II) (S5); and determining a rotor airfoil based on the optimal sample point (S6). The method performs optimized design on aerodynamic characteristics of the airfoil in a state with a changing incoming flow and a changing angle of attack, and can effectively alleviate dynamic stall in this state.
Description

This application claims priority to the Chinese Patent Application No. 201911017404.0, filed with the China National Intellectual Property Administration (CNIPA) on Oct. 24, 2019, and entitled “METHOD AND SYSTEM FOR DETERMINING HELICOPTER ROTOR AIRFOIL”, which is incorporated herein by reference in its entirety.


TECHNICAL FIELD

The present disclosure relates to the technical field of rotor airfoil design, and in particular, to a method and system for determining a helicopter rotor airfoil.


BACKGROUND

A rotor is a key lift and control component of a helicopter. Aerodynamic characteristics of an airfoil have important influence on aerodynamic performance and control characteristics of the rotor and the helicopter. When the helicopter flies forward, rotor blades perform rotation, flapping, and pitching simultaneously. In this way, blade sections are under periodically changing incoming flows and angles of attack, and dynamic stall may occur on a trailing side, resulting in a sharp drop in lift and a sudden increase in drag and moments. This is not conducive to flight. Therefore, reducing or even avoiding the occurrence of dynamic airfoil stall in a rotor environment plays an important role in improving the aerodynamic performance and control performance of the rotor.


A lot of research work has been done on design of a helicopter rotor airfoil at home and abroad. Some aviation-developed countries have established corresponding helicopter-specific airfoil libraries, including Boeing's VR series, OA series in France, TsAGI series in Russia, and the like. For the rotor airfoil, NASA, ONERA and Boeing have proposed their own airfoil design indicators. Researchers at home and abroad have also been devoted to designing the rotor airfoil. However, all these mainly consider static characteristics of the airfoil and dynamic characteristics of the airfoil in a condition with a constant velocity and a changing angle of attack. Currently, for an actual rotor environment, rotor airfoil design is still unavailable for a condition with a changing incoming flow and a changing angle of attack.


SUMMARY

In view of this, it is necessary to provide a method and system for determining a helicopter rotor airfoil, to perform, based on an actual working environment of a rotor airfoil, optimization design on aerodynamic characteristics of the airfoil in a condition with a changing incoming flow and a changing angle of attack, thereby effectively alleviating dynamic stall in this condition.


To achieve the above objective, the present disclosure provides the following solutions:


The present disclosure provides a method for determining a helicopter rotor airfoil, including:


randomly generating a sample point by using a Latin hypercube sampling (LHS) method;


determining characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a class shape transformation (CST) method;


performing dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a computational fluid dynamics (CFD) method, to obtain a flow field characteristic of the airfoil, where the flow field characteristic includes a lift coefficient, a drag coefficient, and a moment coefficient;


establishing a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and training the mapping relationship by using a maximum likelihood estimation method and an expected improvement (EI) criterion, to obtain a trained mapping relationship;


determining an optimal sample point based on the trained mapping relationship by using Non-dominated Sorting Genetic Algorithm II (NSGA-II), where the optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic includes an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient; and


determining a rotor airfoil based on the optimal sample point.


The present disclosure further provides a system for determining a helicopter rotor airfoil, including:


a sample point generation module, configured to randomly generate a sample point by using an LHS method;


an airfoil surface equation determining module, configured to determine characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a CST method;


a flow field characteristic calculation module, configured to perform dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a CFD method, to obtain a flow field characteristic of the airfoil, where the flow field characteristic includes a lift coefficient, a drag coefficient, and a moment coefficient;


a mapping relationship calculation module, configured to establish a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and train the mapping relationship by using a maximum likelihood estimation method and an EI criterion, to obtain a trained mapping relationship;


a multi-objective optimization module, configured to determine an optimal sample point based on the trained mapping relationship by using NSGA-II, where the optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic includes an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient; and


an airfoil determining module, configured to determine a rotor airfoil based on the optimal sample point.


Compared with the prior art, the present disclosure has the following beneficial effects:


The present disclosure provides the method and system for determining a helicopter rotor airfoil. An optimization method for determining dynamic characteristics of the airfoil in the condition with the changing incoming flow and the changing angle of attack is constructed based on the CST method, the LHS method, the CFD method, and a multi-objective genetic algorithm NSGA-II. An optimized airfoil determined by using the method or the system can effectively suppress a leading edge vortex of the airfoil in the condition with the changing incoming flow and the changing angle of attack, so that lift, drag, and moments of the airfoil are not affected when the leading edge vortex moves on or falls off a surface of the airfoil, thereby avoiding a sharp drop in the lift and divergence of the drag and the moment, and greatly improving dynamic aerodynamic characteristics of the airfoil.





BRIEF DESCRIPTION OF THE DRAWINGS

To describe the technical solutions in the embodiments of the present disclosure or in the prior art more clearly, the following briefly describes the accompanying drawings required for describing the embodiments. Apparently, the accompanying drawings in the following description show merely some embodiments of the present disclosure, and a person of ordinary skill in the art may still derive other drawings from these accompanying drawings without creative efforts.



FIG. 1 is a flowchart of a method for determining a helicopter rotor airfoil according to an embodiment of the present disclosure;



FIG. 2 is a schematic diagram of a moving embedded grid according to an embodiment of the present disclosure;



FIG. 3 is a flowchart of NSGA-II—based calculation according to an embodiment of the present disclosure;



FIG. 4 is a schematic diagram of a coordinate system and a profile of an airfoil according to an embodiment of the present disclosure;



FIG. 5 compares an SC1095 airfoil and an optimized airfoil;



FIG. 6 compares lift coefficients of an SC1095 airfoil and an optimized airfoil;



FIG. 7 compares drag coefficients of an SC1095 airfoil and an optimized airfoil;



FIG. 8 compares moment coefficients of an SC1095 airfoil and an optimized airfoil;



FIG. 9 compares streamline charts of an SC1095 airfoil and an optimized airfoil;



FIG. 10 compares pressure coefficient distribution curves of an SC1095 airfoil and an optimized airfoil; and



FIG. 11 is a schematic structural diagram of a system for determining a helicopter rotor airfoil according to an embodiment of the present disclosure.





DETAILED DESCRIPTION

The technical solutions in the embodiments of the present disclosure are clearly and completely described below with reference to the accompanying drawings in the embodiments of the present disclosure. Apparently, the described embodiments are only a part rather than all of the embodiments of the present disclosure. All other embodiments obtained by a person of ordinary skill in the art based on the embodiments of the present disclosure without creative efforts shall fall within the protection scope of the present disclosure.


To make objective, features, and advantages of the present disclosure more comprehensible, the following describes the present disclosure in more detail with reference to accompanying drawings and specific implementations.


Refer to FIG. 1. A method for determining a helicopter rotor airfoil in an embodiment includes the following steps.


Step S1: Randomly generate a sample point by using an LHS method.


The step S1 is specifically: generating a sample point xji by using the LHS method, where j represents a number of the sample point, and i represents a number of a variable. It is assumed that m samples are to be extracted from n-dimensional vector space. LHS includes the following steps: (1) dividing each dimension into m non-overlapping intervals, so that all the intervals have a same probability; (2) randomly extracting a point from each interval of each dimension; and (3) extracting, from each dimension, the point selected in the step (2), to constitute a vector.


Step S2: Determine characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a CST method. In this step, the CST method is used to fit a baseline airfoil and generate an airfoil in sample space. Based on the CST method, the upper and lower airfoil surfaces of the airfoil may be represented as follows:






{







ψ
u



(
ζ
)


=



C


(
ζ
)


·




k
=
0

n




A
u
k

·


S
k



(
ζ
)





+


ζ
·
Δ







z
u











ψ
l



(
ζ
)


=



C


(
ζ
)


·




k
=
0

n




A
l
k

·


S
k



(
ζ
)





+


ζ
·
Δ







z
l













where ψu=yu|c,ψl=yl|c; ψu represents a dimensionless vertical coordinate of a point on the upper airfoil surface of the airfoil; ψl represents a dimensionless vertical coordinate of a point on the lower airfoil surface of the airfoil; yu represents a vertical coordinate of the point on the upper airfoil surface of the airfoil; yl represents a vertical coordinate of the point on the lower airfoil surface of the airfoil; c represents a chord length of the airfoil; ζ=x|c; ζ represents a dimensionless horizontal coordinate of a point on the airfoil; x represents a horizontal coordinate of the point on the airfoil; Δzu=yuTE|c, Δzl=ylTE|c, Δzu represents a dimensionless trailing edge thickness of the upper airfoil surface of the airfoil; Δzl represents a dimensionless trailing edge thickness of the lower surface of the airfoil; yuTE represents a vertical coordinate of a trailing edge point on the upper airfoil surface of the airfoil; ylTE represents a vertical coordinate of a trailing edge point of the lower airfoil surface of the airfoil; C(ζ)=ζN1·(1−ζ)N2; C(ζ) represents a class function; N1 and N2 are both constants; for a conventional airfoil with a rounded head and a pointed trailing edge, N1=0.5, and N2=1; Auk represents a coefficient for controlling a profile of upper airfoil surface of the airfoil; Alk represents a coefficient for controlling a profile of the lower airfoil surface of the airfoil; k represents a coefficient number; n represents a total quantity of coefficients; Sk(ζ) represents a shape function, and a Bernstein polynomial is used as the shape function in this embodiment;









S
k



(
ζ
)


-



b
!



a
!




(

b
-
k

)

!



·

ζ
k

·


(

1
-
ζ

)


b
-
k




;




a represents a number of a polynomial; b represents an order of the polynomial; and a=0, 1, . . . , b. Different geometric profiles of the airfoil can be obtained by adjusting values of the coefficients Auk and Alk. Therefore, the coefficients Auk and Alk are used as design variables for optimization design of the airfoil, and are optimized to obtain an optimized airfoil profile satisfying a requirement. In this embodiment, a six-order Bernstein polynomial is used as the shape function. A coefficient of the polynomial is a design variable, and is optimized to obtain a new optimized airfoil.


The baseline airfoil, namely, an SC1095 airfoil, is fitted to obtain polynomial coefficients Auk and Alk of the airfoil, where k=0, 1, . . . , 6. The sample point generated by the LHS method is used to obtain polynomial coefficients of a sample airfoil numbered j (the jth sample point xji): Aujk=ƒ(Auik,xji) and Aljk=ƒ(Alik,xji), where i represents the number of the variable, and ƒ(g) represents a transformation function. The foregoing polynomial coefficients are substituted into the characterization equations of the upper and lower airfoil surfaces of the airfoil to obtain coordinates of all sample points.


Step S3: Perform dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a CFD method, to obtain a flow field characteristic of the airfoil, where the flow field characteristic includes a lift coefficient, a drag coefficient, and a moment coefficient.


In this step, the CFD method is used to perform dynamic characteristic simulation on the airfoil in the sample space. The following typical working condition of the airfoil is selected as a simulation condition: Ma=0.3+0.12sin(ωt), α=10°−8° sin(ωt ), k=0.069, and Re=3.75×106, where Ma represents the Mach number, ω represents a circular frequency of airfoil oscillation, t represents physical time, k represents a reduced frequency, and Re represents the Reynolds number. Specifically:


31) Determine coordinates (ζ,ψ) of the sample point according to the characterization equations of the upper and lower airfoil surfaces.


32) Determine an initial grid, where the initial grid is generated based on the coordinates of the sample point.


33) Generate, based on the initial grid and by using a Poisson equation—based elliptic equation grid generation method, a C-shaped body-fitted grid around the airfoil, where the Poisson equation is as follows:






{







ξ
ζζ

+

ξ
ψψ


=

P


(

ξ
,
η

)










η
ζζ

+

η
ψψ


=

Q


(

ξ
,
η

)












where (ξ,η) represents a point in a curvilinear coordinate system under a calculation plane, there is a mapping relationship between (ξ,η) and a point (ζ,ψ)in a linear coordinate system under a physical plane,








ξ
ζζ

=




2


ξ




ζ
2




,


ξ
ψψ

=




2


ξ




ψ
2




,


η
ζζ

=




2


η




ζ
2




,


η
ψψ

=




2


η




ψ
2




,




P(ξ,η) represents an orthogonality control function, and Q(ξ,η) represents a spacing control function. Grid orthogonality can be controlled by changing P(ξ,η), and grid spacing can be controlled by changing Q(ξ,η). The Poisson equation can be used to finally obtain the C-shaped body-fitted grid that is of the airfoil and whose orthogonality and spacing satisfy calculation requirements. In this embodiment, a moving embedded grid method is used. Rotation and translational motion are performed on the airfoil grid in a background grid to periodically change a velocity of an incoming flow and an angle of attack, as shown in FIG. 2.


34) In this embodiment, use a Reynolds-averaged Navier-Stokes (RANS) equation in integral form as a rotor flow field solving control equation, and when solving the control equation, read coordinates of a grid point of the airfoil first, and then calculate an area of a grid cell, a volume of the grid cell, and a normal vector of a grid cell surface. Therefore, in this step, an area and a volume of each grid cell in the C-shaped body-fitted grid are calculated.


35) Determine a fluid control equation (the RANS equation) in integral form in an inertial coordinate system based on the area, the volume, a conservative variable, a convective flux term, and a viscous flux term. The fluid control equation is as follows:












t











Ω


Wd





Ω

+








Ω




(


F
c

-

F
v


)


dS


=
0




where t represents the physical time, Fc represents the convective flux term, Fv represents the viscous flux term, Ω represents the volume of the grid cell, S represents the area of the grid cell, W represents the conservative variable,







W
=

[



ρ





ρ





u






ρ





v






ρ





E




]


,






F
c

=

[




ρ






V
r








ρ






uV
r


+


n
x


p








ρ






vV
r


+


n
y


p








ρ






HV
r


+


V
t


p





]


,






F
v

=

[



0







n
x



τ
xx


+


n
y



τ
xy










n
x



τ
yx


+


n
y



τ
yy










n
x



Θ
x


+


n
y



Θ
y






]


,




ρ represents air density, E represents total energy, H represents total enthalpy, p represents pressure, u represents a velocity component in an x-axis direction, v represents a velocity component in a Y-axis direction, Vr represents a velocity of relative motion, Vt represents a grid moving velocity, (nx,ny) represents a unit normal vector of the grid cell surface, τxx represents a viscous stress component acting on a plane perpendicular to an x axis and along the x-axis direction, τxy represents a viscous stress component acting on the plane perpendicular to the x axis and along the y-axis direction, τyy represents a viscous stress component acting on a plane perpendicular to a y axis and along the y-axis direction, Θx represents a viscous stress term, and Θy represents a heat conduction term.


36) Solve the fluid control equation, where a solution of the fluid control equation is






W
=


[



ρ





ρ





u






ρ





v






ρ

E




]

.





The flow field characteristic of the airfoil is calculated based on the solution of the fluid control equation.


Step S4: Establish a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and train the mapping relationship by using a maximum likelihood estimation method and an EI criterion, to obtain a trained mapping relationship. The step S4 specifically includes the following substeps.


41) Use the Kriging model to establish the mapping relationship between the sample point and the flow field characteristic, in other words, a mapping relationship between xji and each of the lift coefficient Cl, the drag coefficient Cd, and the moment coefficient Cm. Based on the Kriging model, predicted values Ĉl, Ĉd, Ĉm of the lift coefficient, the drag coefficient, and the moment coefficient at an unknown point x may be represented as a linear sample point combination:






Ĉl(x)=rlTRl−1Yl−(FTRl−1rl−f)T(FTRl−1F)−1FTRl−1Yl,






Ĉd(x)=rdTRd−1Yd−(FTRd−1rd−f)T(FTRd−1F)−1FTRd−1Yd,






Ĉm(x)=rmTRm−1Ym−(FTRm−1rm−f)T(FTRm−1F)−1FTRm−1Ym,


where Ĉl represents the predicted value of the lift coefficient, Ĉd represents the predicted value of the drag coefficient, Ĉm represents the predicted value of the moment coefficient, Yl represents a vector constituted by lift coefficients Cl corresponding to sample points, Yd represents a vector constituted by drag coefficients Cd corresponding to the sample points, Ym represents a vector constituted by moment coefficients Cm corresponding to the sample points, Yl=[Cl1 Cl2 L Clns]T, Yd=[Cd1 Cd2 L Cdns]T, Ym=[Cm1 Cm2 L Cmns]T, ns represents a quantity of the sample points,







f
=

[





f
1



(
x
)






M






f
p



(
x
)





]


,

F
=

[





f
1



(

x
1

)




L




f
p



(

x
1

)






M


O


M






f
1



(

x
ns

)




L




f
p



(

x
ns

)





]


,




ƒl(g), K, ƒp(g) represents a known polynomial function, x1, x2,and xns are all vectors constituted by sample points xji, x1=[x11 x12 K x1e], x2=[x21 x22 K x2e], xns=[xns1 xns2 K xnse], e represents a vector dimension, i represents the number of the variable, j represents a number of the sample point,








R
l

=

[





R
l



(


x
1

,

x
1


)




L




R
l



(


x
1

,

x
ns


)






M


O


M






R
l



(


x
ns

,

x
1


)




L




R
l



(


x
ns

,

x
ns


)





]


,


r
l

=

[





R
l



(

x
,

x
1


)








R
l



(

x
,

x
2


)






M






R
l



(

x
,

x
ns


)





]


,



R
l



(


x

j





1


,

x

j





2



)


=




i
=
1

n








(


-

θ
li







x

j





1

i

-

x

j





2

i





)


p
li




,




θli and Pli represent unknown parameters corresponding to the lift coefficient, θli and Pli are obtained by calculating







max
(



-

ns
2




ln
(




(


Y
l

-

F







β
^

l



)

T




R
l

-
1




(


Y
l

-

F







β
^

l



)




n
s


)


-


1
2



ln


[



R
l



]




)

,



β
^

l

=



(


F
T



R
l

-
1



F

)


-
1




F
T



R
l

-
1




Y
l



,


R
d

=

[





R
d



(


x
1

,

x
1


)




L




R
d



(


x
1

,

x
ns


)






M


O


M






R
d



(


x
ns

,

x
1


)




L




R
d



(


x
ns

,

x
ns


)





]


,


r
d

=

[





R
d



(

x
,

x
1


)








R
d



(

x
,

x
2


)






M






R
d



(

x
,

x
ns


)





]


,



R
d



(


x

j





1


,

x

j





2



)


=




i
=
1

n








(


-

θ
di







x

j





1

i

-

x

j





2

i





)


p
di




,




θdi and Pdi represent unknown parameters corresponding to the drag coefficient, θdi and Pdi are obtained by calculating







max
(



-


n
s

2




ln
(




(


Y
d

-

F







β
^

d



)

T




R
d

-
1




(


Y
d

-

F







β
^

d



)




n
s


)


-


1
2



ln


[



R
d



]




)

,



β
^

d

=



(


F
T



R
d

-
1



F

)


-
1




F
T



R
d

-
1




Y
d



,


R
m

=

[





R
m



(


x
1

,

x
1


)




L




R
m



(


x
1

,

x
ns


)






M


O


M






R
m



(


x
ns

,

x
1


)




L




R
m



(


x
ns

,

x
ns


)





]


,


r
m

=

[





R
m



(

x
,

x
1


)








R
m



(

x
,

x
2


)






M






R
m



(

x
,

x
ns


)





]


,



R
m



(


x

j





1


,

x

j





2



)


=




i
=
1

n








(


-

θ
mi







x

j





1

i

-

x

j





2

i





)


p
mi




,




θmi and Pmi represent unknown parameters corresponding to the moment coefficient, θmi and Pmi are obtained by calculating







max
(



-


n
s

2




ln
(




(


Y
m

-

F







β
^

m



)

T




R
m

-
1




(


Y
d

-

F







β
^

m



)




n
s


)


-


1
2



ln


[



R
m



]




)

,



β
^

m

=



(


F
T



R
m

-
1



F

)


-
1




F
T



R
m

-
1




Y
m



,
and




42) Determine an unknown parameter in the mapping relationship by using the maximum likelihood estimation method, to obtain the trained mapping relationship.


43) Improve the trained mapping relationship by using the EI criterion, to obtain the trained mapping relationship.


Specifically, to improve prediction precision of the model, the EI criterion is used to find a new sample point with a maximum value of a mathematical expectation and a new sample point with an optimal objective function value in a process of training the Kriging model, and the two new sample points are used for training optimization of the model, till a convergence condition is satisfied. Ĉl(x) is used as an example, and Ĉd(x) and Ĉm(x) are similar to Ĉl(x). In the EI criterion, it is assumed that Ĉl(x) of any prediction point x is characterized by normal random distribution, in other words, satisfies Ĉl(x)−−N(μ,s2), where μ represents an average predicted value, and s represents a standard deviation. In calculation space, a mathematical expectation of any point can be calculated according to the following formula:







E


[

I


(
x
)


]


=

{






(



C
^



l
min


-


C
^



l


(
x
)




)



Φ
(




C
^



l
min


-


C
^



l


(
x
)




s

)


+


s


(
x
)




ϕ
(




C
^



l
min


-


C
^



l


(
x
)




s

)






s
>
0





0



s
=
0









where Ĉlmin represents a minimum value of the sample point, Φ represents a standard normal distribution function, and ϕ represents a density function of normal distribution. When a maximum value of the mathematical expectation E[I(x)] is found, a point with a relatively small predicted value or a point with relatively low prediction precision of the model can be found. In each interpolation process, an optimal improved point x can always be found. A specific process is as follows: (1) Substitute x into a CST formula to generate a new airfoil. (2) Generate a grid. (3) Calculate dynamic characteristics. (4) Add related data of the airfoil to the sample airfoil space, to re-establish the Kriging model. (5) Train the Kriging model by using the EI criterion, to obtain the optimal improved point x. Repeat the steps (1) to (5) till prediction precision of a surrogate model satisfies a requirement, so that the mapping relationship between xji and Cl, Cd, Cm is established.


Step S5: Determine an optimal sample point based on the trained mapping relationship by using NSGA-II. The optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic includes an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient.


The step S5 specifically includes the following substeps:


501) Randomly generate an initial population whose scale is N.


502) Calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in a current-generation population.


503) Perform fast non-dominated sorting on all individuals in the current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation population, and obtain a current-generation parent population by performing selection, crossover, and mutation based on a genetic algorithm, where a quantity of population generations is equal to a quantity of iterations of the genetic algorithm.


504) Calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the current-generation parent population.


505) Perform fast non-dominated sorting on all individuals in the current-generation parent population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation parent population, and obtain a current-generation child population by performing selection, crossover, and mutation based on the genetic algorithm.


506) Combine the current-generation parent population and the current-generation child population to obtain a combined current-generation population.


507) Calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the combined current-generation population.


508) Perform fast non-dominated sorting on all individuals in the combined current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the combined current-generation population, and calculate an individual crowdedness degree in each non-dominated layer.


509) Select, from the combined current-generation population based on a result of fast non-dominated sorting and the individual crowdedness degree, an individual that satisfies a preset condition, to constitute a next-generation parent population.


510) Determine whether a sequence number of a generation to which the next-generation parent population belongs is equal to a preset maximum quantity of iterations. If the sequence number of the generation to which the next-generation parent population belongs is not equal to the preset maximum quantity of iterations, the next-generation parent population is used as a current-generation parent population, and the step 504) is performed to calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the current-generation parent population. If the sequence number of the generation to which the next-generation parent population belongs is equal to the preset maximum quantity of iterations, an individual in the next-generation parent population is determined as the optimal sample point.


Refer to FIG. 3. Based on the foregoing ideas, a specific process of finding an optimal airfoil by using NSGA-II is as follows: (1) Randomly generate an initial population x1, x2, K, xN whose scale is N. (2) Calculate Cl, Cd, Cm of each individual in the population by using the mapping relationship established based on the Kriging model. (3) Perform non-dominated sorting on an individual in the population based on values of Cl, Cd, Cm, and obtain the first-generation parent population x11, x21, K, xN1 by performing selection, crossover, and mutation based on the genetic algorithm. (4) Repeat the steps (2) and (3) to obtain the second-generation child population x12, x22, K, xN2. (5) Combine the parent population and the child population. (6) Calculate Cl, Cd, Cm of each individual in a combined population by using the mapping relationship established based on the Kriging model. (7) Perform fast non-dominated sorting based on Cl, Cd, Cm, and calculate an individual crowdedness degree in each non-dominated layer. (8) Select a proper individual based on a non-dominated relationship and the individual crowdedness degree, to constitute a new parent population (x1it, x2it, K, xNit). (9) Repeat the steps (2) and (3) to obtain a new child population x1it+1, x2it+1, K, xNit+1. (10) Repeat the steps (5) to (9), till a quantity of iterations reaches a preset value and a final optimal sample point is obtained.


Step S6: Determining a rotor airfoil based on the optimal sample point.


In actual application, the optimized airfoil finally obtained by performing the steps S1 to S6 is characterized by a large leading edge radius on an upper airfoil surface, a flat transition part from 0.2 c to 0.4 c on the upper surface, a concave lower surface, and a slightly upturned trailing edge. The optimized airfoil is manufactured based on data provided in Table 1 A chord length of the airfoil is 1.









TABLE 1







Data of the optimized airfoil












Sequence


Sequence




number of a
Horizontal
Vertical
number of a
Horizontal
Vertical


coordinate
coordinate
coordinate
coordinate
coordinate
coordinate


point on the
(x/c) on the
(y/c) on the
point on the
(x/c) on the
(y/c) on the


lower airfoil
lower airfoil
lower airfoil
upper airfoil
upper airfoil
upper airfoil


surface
surface
surface
surface
surface
surface















1
1.00000
0.00000
1
1.00000
0.00000


2
0.99017
−0.00179
2
0.99008
0.00047


3
0.97206
−0.00458
3
0.97142
0.00007


4
0.94635
−0.00766
4
0.94499
0.00320


5
0.91376
−0.01057
5
0.91183
0.00961


6
0.87506
−0.01315
6
0.87275
0.01872


7
0.83104
−0.01537
7
0.82834
0.02946


8
0.78244
−0.01709
8
0.77916
0.04073


9
0.73002
−0.01807
9
0.72585
0.05173


10
0.67457
−0.01811
10
0.66915
0.06189


11
0.61685
−0.01723
11
0.60982
0.07079


12
0.55762
−0.01584
12
0.54867
0.07794


13
0.49762
−0.01455
13
0.48649
0.08296


14
0.43761
−0.01395
14
0.42418
0.08588


15
0.37835
−0.01434
15
0.36260
0.08733


16
0.32061
−0.01561
16
0.30259
0.08818


17
0.26515
−0.01736
17
0.24493
0.08891


18
0.21272
−0.01914
18
0.19043
0.08895


19
0.16407
−0.02060
19
0.13994
0.08651


20
0.11994
−0.02148
20
0.09471
0.07929


21
0.08109
−0.02144
21
0.05669
0.06610


22
0.04834
−0.01991
22
0.02789
0.04813


23
0.02268
−0.01599
23
0.00957
0.02849


24
0.00577
−0.00900
24
0.00154
0.01129


25
0.00011
−0.00100
25
0.00005
0.00099


26
0.00000
0.00000
26
0.00000
0.00000









A coordinate system and a profile of the airfoil are shown in FIG. 4. FIG. 5 compares an SC1095 airfoil and the optimized airfoil. It can be learned from FIG. 5 that the optimized airfoil has a larger camber, the leading edging radius on the upper airfoil surface of the optimized airfoil is obviously larger than that of the SC1095 airfoil, the flat transition part is available in the middle of the upper surface of the optimized airfoil, and the optimized airfoil is characterized by the concave lower surface and the slightly upturned trailing edge. FIG. 6 to FIG. 8 compare aerodynamic characteristics of the two airfoils in a coupled condition with a changing incoming flow and a changing angle of attack. It can be learned from FIG. 6 to FIG. 8 that, compared with the SC1095 airfoil, the optimized airfoil has been greatly improved in drag and moment characteristics, thereby alleviating drag divergence and moment divergence. In addition, an area of a lift coefficient hysteresis loop of the optimized airfoil is remarkably reduced, and this indicates that the optimized airfoil can effectively improve dynamic aerodynamic characteristics of the airfoil in a rotor environment and suppress dynamic stall. FIG. 9 to FIG. 10 compare flow field characteristics of the two airfoils under three angles of attack in the coupled condition with the changing incoming flow and the changing angle of attack. It can be learned from FIG. 9 to FIG. 10 that, when the angle of attack (α) is 17.56° in a pitching-up process, for both the SC1095 airfoil and the optimized airfoil, a vortex is generated only near a trailing edge, and a streamline attaches an airfoil surface at a leading edge and in the middle. As the angle of attack increases to 18.00°, a vortex is generated at the leading edge of the SC1095 airfoil. In addition, as the airfoil pitches down, the leading edge vortex is enhanced and gradually moves to the trailing edge until it sheds from the airfoil surface. For the optimized airfoil, an adverse pressure gradient of the leading edge is alleviated by changing the profile of the optimized airfoil, thereby suppressing the leading edge vortex. Therefore, the optimized airfoil as better dynamic characteristics.


Compared with the SC1095 airfoil, the optimized airfoil has a better dynamic stall characteristic. Reasons are as follows: (1) The optimized airfoil has a larger leading edge radius on the upper airfoil surface, and this can effectively alleviate the adverse pressure gradient of the leading edge of the airfoil (the optimized airfoil has a smaller peak negative pressure value at the leading edge, as shown in FIG. 10(a) and FIG. 10(b)), and weaken strength of the leading edge vortex and even suppress the leading edge vortex (streamline, as shown in FIG. 9), thereby alleviating and even avoiding a sharp drop in lift and divergence of drag and moments when the leading edge vortex moves or falls off. (2) The flat transition part from a position 0.2 c to a position 0.4 c is available on the upper surface of the optimized airfoil, and this can postpone drag divergence and moment divergence. (3) Positions 0.1 c to 0.7 c on the lower surface of the optimized airfoil are concave, increasing the camber of the airfoil and improving the lift of the airfoil. (4) The trailing edge of the optimized airfoil is slightly upturned, and this can postpone drag divergence and moment divergence, and can improve a reattaching characteristic of an airflow in the case of a large angle of attack.


According to the method for determining a helicopter rotor airfoil in this embodiment, in the condition with the changing incoming flow and the changing angle of attack, the leading edge vortex is suppressed in the entire motion process, and the vortex is generated only near the trailing edge in the case of the large angle of attack. In this way, the lift, the drag, and the moment of the airfoil are no longer affected by moving or shedding of the leading edge vortex in the dynamic stall process.


The present disclosure further provides a system for determining a helicopter rotor airfoil. FIG. 11 is a schematic structural diagram of a system for determining a helicopter rotor airfoil according to an embodiment of the present disclosure. Refer to FIG. 11. The system for determining a helicopter rotor airfoil in this embodiment includes a sample point generation module 1101, an airfoil surface equation determining module 1102, a flow field characteristic calculation module 1103, a mapping relationship calculation module 1104, a multi-objective optimization module 1105, and an airfoil determining module 1106.


The sample point generation module 1101 is configured to randomly generate a sample point by using an LHS method.


The airfoil surface equation determining module 1102 is configured to determine characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a CST method.


The flow field characteristic calculation module 1103 is configured to perform dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a CFD method, to obtain a flow field characteristic of the airfoil, where the flow field characteristic includes a lift coefficient, a drag coefficient, and a moment coefficient.


The mapping relationship calculation module 1104 is configured to establish a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and train the mapping relationship by using a maximum likelihood estimation method and an EI criterion, to obtain a trained mapping relationship.


The multi-objective optimization module 1105 is configured to determine optimal sample point based on the trained mapping relationship by using NSGA-II, where the optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic includes an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient.


The airfoil determining module 1106 is configured to determine a rotor airfoil based on the optimal sample point.


Each embodiment in this specification is described in a progressive manner, each embodiment focuses on the difference from other embodiments, and the same and similar parts between the embodiments may refer to each other. For the system disclosed in the embodiments, since the system corresponds to the method disclosed in the embodiments, the system is simply described, and reference can be made to the method description.


In this specification, several specific embodiments are used for illustration of the principles and implementations of the present disclosure. The description of the foregoing embodiments is used to help illustrate the method of the present disclosure and the core ideas thereof. In addition, those of ordinary skill in the art can make various modifications in terms of specific implementations and scope of application in accordance with the ideas of the present disclosure. In conclusion, the content of this specification shall not be construed as a limitation on the present disclosure.

Claims
  • 1. A method for determining a helicopter rotor airfoil, comprising: randomly generating a sample point by using a Latin hypercube sampling (LHS) method;determining characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a class shape transformation (CST) method;performing dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a computational fluid dynamics (CFD) method, to obtain a flow field characteristic of the airfoil, wherein the flow field characteristic comprises a lift coefficient, a drag coefficient, and a moment coefficient;establishing a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and training the mapping relationship by using a maximum likelihood estimation method and an expected improvement (EI) criterion, to obtain a trained mapping relationship;determining an optimal sample point based on the trained mapping relationship by using Non-dominated Sorting Genetic Algorithm II (NSGA-II), wherein the optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic comprises an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient; anddetermining a rotor airfoil based on the optimal sample point.
  • 2. The method for determining a helicopter rotor airfoil according to claim 1, wherein the characterization equations of the upper and lower airfoil surfaces are as follows:
  • 3. The method for determining a helicopter rotor airfoil according to claim 1, wherein the performing dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a CFD method, to obtain a flow field characteristic of the airfoil specifically comprises: determining coordinates of the sample point according to the characterization equations of the upper and lower airfoil surfaces;determining an initial grid, wherein the initial grid is generated based on the coordinates of the sample point;generating, based on the initial grid and by using a Poisson equation—based elliptic equation grid generation method, a C-shaped body-fitted grid around the airfoil, wherein the Poisson equation is as follows:
  • 4. The method for determining a helicopter rotor airfoil according to claim 1, wherein the establishing a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and training the mapping relationship by using a maximum likelihood estimation method and an EI criterion, to obtain a trained mapping relationship specifically comprises: establishing the mapping relationship between the sample point and the flow field characteristic by using the Kriging model, wherein the mapping relationship is as follows: Ĉl(x)=rlTRl−1Yl−(FTRl−1rl−f)T(FTRl−1F)−1FTRl−1Yl,Ĉd(x)=rdTRd−1Yd−(FTRd−1rd−f)T(FTRd−1F)−1FTRd−1Yd,Ĉm(x)=rmTRm−1Ym−(FTRm−1rm−f)T(FTRm−1F)−1FTRm−1Ym,wherein Ĉl represents a predicted value of the lift coefficient, Ĉd represents a predicted value of the drag coefficient, Ĉm represents a predicted value of the moment coefficient, Yl represents a vector constituted by lift coefficients Cl corresponding to sample points, Yd represents a vector constituted by drag coefficients Cd corresponding to the sample points, Ym represents a vector constituted by moment coefficients Cm corresponding to the sample points, Yl=[Cl1 Cl2 L Clns]T, Yd=[Cd1 Cd2 L Cdns]T, Ym=[Cm1 Cm2 L Cmns]T, ns represents a quantity of the sample points,
  • 5. The method for determining a helicopter rotor airfoil according to claim 1, wherein the determining an optimal sample point based on the trained mapping relationship by using NSGA-II specifically comprises: randomly generating an initial population whose scale is N;calculating, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in a current-generation population;performing fast non-dominated sorting on all individuals in the current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation population, and obtaining a current-generation parent population by performing selection, crossover, and mutation based on a genetic algorithm, wherein a quantity of population generations is equal to a quantity of iterations of the genetic algorithm;calculating, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the current-generation parent population;performing fast non-dominated sorting on all individuals in the current-generation parent population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation parent population, and obtaining a current-generation child population by performing selection, crossover, and mutation based on the genetic algorithm;combining the current-generation parent population and the current-generation child population to obtain a combined current-generation population;calculating, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the combined current-generation population;performing fast non-dominated sorting on all individuals in the combined current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the combined current-generation population, and calculating an individual crowdedness degree in each non-dominated layer;selecting, from the combined current-generation population based on a result of fast non-dominated sorting and the individual crowdedness degree, an individual that satisfies a preset condition, to constitute a next-generation parent population;determining whether a sequence number of a generation to which the next-generation parent population belongs is equal to a preset maximum quantity of iterations; andif the sequence number of the generation to which the next-generation parent population belongs is not equal to the preset maximum quantity of iterations, using the next-generation parent population as a current-generation parent population, and returning to the step of calculating, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the current-generation parent population; orif the sequence number of the generation to which the next-generation parent population belongs is equal to the preset maximum quantity of iterations, determining an individual in the next-generation parent population as the optimal sample point.
  • 6. A system for determining a helicopter rotor airfoil, comprising: a sample point generation module, configured to randomly generate a sample point by using an LHS method;an airfoil surface equation determining module, configured to determine characterization equations of upper and lower airfoil surfaces of an airfoil based on the sample point by using a CST method;a flow field characteristic calculation module, configured to perform dynamic characteristic simulation on the airfoil according to the characterization equations of the upper and lower airfoil surfaces by using a CFD method, to obtain a flow field characteristic of the airfoil, wherein the flow field characteristic comprises a lift coefficient, a drag coefficient, and a moment coefficient;a mapping relationship calculation module, configured to establish a mapping relationship between the sample point and the flow field characteristic by using a Kriging model, and train the mapping relationship by using a maximum likelihood estimation method and an EI criterion, to obtain a trained mapping relationship;a multi-objective optimization module, configured to determine an optimal sample point based on the trained mapping relationship by using NSGA-II, wherein the optimal sample point corresponds to an optimal flow field characteristic, and the optimal flow field characteristic comprises an optimal lift coefficient, an optimal drag coefficient, and an optimal moment coefficient; andan airfoil determining module, configured to determine a rotor airfoil based on the optimal sample point.
  • 7. The system for determining a helicopter rotor airfoil according to claim 6, wherein the characterization equations of the upper and lower airfoil surfaces are as follows:
  • 8. The system for determining a helicopter rotor airfoil according to claim 6, wherein the flow field characteristic calculation module specifically comprises: a coordinate determining unit, configured to determine coordinates of the sample point according to the characterization equations of the upper and lower airfoil surfaces;an initial grid determining unit, configured to determine an initial grid, wherein the initial grid is generated based on the coordinates of the sample point;a body-fitted grid generation unit, configured to generate, based on the initial grid and by using a Poisson equation—based elliptic equation grid generation method, a C-shaped body-fitted grid around the airfoil, wherein the Poisson equation is as follows:
  • 9. The system for determining a helicopter rotor airfoil according to claim 6, wherein the mapping relationship calculation module specifically comprises: a mapping relationship determining unit, configured to establish the mapping relationship between the sample point and the flow field characteristic by using the Kriging model, wherein the mapping relationship is as follows: Ĉl(x)=rlTRl−1Yl−(FTRl−1rl−f)T(FTRl−1F)−1FTRl−1Yl,Ĉd(x)=rdTRd−1Yd−(FTRd−1rd−f)T(FTRd−1F)−1FTRd−1Yd,Ĉm(x)=rmTRm−1Ym−(FTRm−1rm−f)T(FTRm−1F)−1FTRm−1Ym,wherein Ĉl represents a predicted value of the lift coefficient, Ĉd represents a predicted value of the drag coefficient, Ĉm represents a predicted value of the moment coefficient, Yl represents a vector constituted by lift coefficients Cl corresponding to sample points, Yd represents a vector constituted by drag coefficients Cd corresponding to the sample points, Ym represents a vector constituted by moment coefficients Cm corresponding to the sample points, Yl=[Cl1 Cl2 L Clns]T, Yd=[Cd1 Cd2 L Cdns]T, Ym=[Cm1 Cm2 L Cmns]T, ns represents a quantity of the sample points,
  • 10. The system for determining a helicopter rotor airfoil according to claim 6, wherein the multi-objective optimization module specifically comprises: an initial-population generation unit, configured to randomly generate an initial population whose scale is N;a first calculation unit, configured to calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in a current-generation population;a first sorting unit, configured to perform fast non-dominated sorting on all individuals in the current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation population, and obtain a current-generation parent population by performing selection, crossover, and mutation based on a genetic algorithm, wherein a quantity of population generations is equal to a quantity of iterations of the genetic algorithm;a second calculation unit, configured to calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the current-generation parent population;a second sorting unit, configured to perform fast non-dominated sorting on all individuals in the current-generation parent population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the current-generation parent population, and obtain a current-generation child population by performing selection, crossover, and mutation based on the genetic algorithm;a population combination unit, configured to combine the current-generation parent population and the current-generation child population to obtain a combined current-generation population;a third calculation unit, configured to calculate, by using the trained mapping relationship, a lift coefficient, a drag coefficient, and a moment coefficient that correspond to each individual in the combined current-generation population;a second sorting unit, configured to perform fast non-dominated sorting on all individuals in the combined current-generation population based on the lift coefficient, the drag coefficient, and the moment coefficient that correspond to each individual in the combined current-generation population, and calculate an individual crowdedness degree in each non-dominated layer;an individual selection unit, configured to select, from the combined current-generation population based on a result of fast non-dominated sorting and the individual crowdedness degree, an individual that satisfies a preset condition, to constitute a next-generation parent population; anda determining unit, configured to determine whether a sequence number of a generation to which the next-generation parent population belongs is equal to a preset maximum quantity of iterations; if the sequence number of the generation to which the next-generation parent population belongs is not equal to the preset maximum quantity of iterations, use the next-generation parent population as a current-generation parent population, and perform the operation performed by the second calculation unit; or if the sequence number of the generation to which the next-generation parent population belongs is equal to the preset maximum quantity of iterations, determine an individual in the next-generation parent population as the optimal sample point.
Priority Claims (1)
Number Date Country Kind
201911017404.0 Oct 2019 CN national
PCT Information
Filing Document Filing Date Country Kind
PCT/CN2020/107732 8/7/2020 WO 00