CN111199105B - Flapping wing motion parameter optimization method - Google Patents

Flapping wing motion parameter optimization method Download PDF

Info

Publication number
CN111199105B
CN111199105B CN202010005675.0A CN202010005675A CN111199105B CN 111199105 B CN111199105 B CN 111199105B CN 202010005675 A CN202010005675 A CN 202010005675A CN 111199105 B CN111199105 B CN 111199105B
Authority
CN
China
Prior art keywords
aerodynamic
motion
precision
flapping
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010005675.0A
Other languages
Chinese (zh)
Other versions
CN111199105A (en
Inventor
郑鸿宇
谢芳芳
季廷炜
朱灶旭
郑耀
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN202010005675.0A priority Critical patent/CN111199105B/en
Publication of CN111199105A publication Critical patent/CN111199105A/en
Application granted granted Critical
Publication of CN111199105B publication Critical patent/CN111199105B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Abstract

The invention relates to the technical field of computer simulation, in particular to a flapping wing motion parameter optimization method. The method comprises the following steps: the method comprises the steps of presetting sampling points in a parameter space by using experimental design, calculating aerodynamic coefficients corresponding to the sampling points by using a numerical simulation method, constructing a proxy model by using high-precision and low-precision samples, optimizing motion parameters of the flapping wings by combining Bayesian optimization, and optimizing convergence judgment. And when the real aerodynamic coefficient of the optimal point does not meet the error requirement, taking the calculation result of the optimal point as a high-precision sample, selecting a candidate sampling point by using an active learning strategy, calculating the corresponding aerodynamic coefficient as a low-precision sample through low-precision numerical simulation, and updating the correction agent model. The invention effectively utilizes data samples with different precisions. And an optimization process is integrated, and the motion parameters of the flapping wings are scientifically and efficiently optimized according to the target aerodynamic coefficient so as to realize specific aerodynamic performance.

Description

Flapping wing motion parameter optimization method
Technical Field
The invention relates to the technical field of computer simulation, in particular to a flapping wing motion parameter optimization method.
Background
With the development of bionics, the aerodynamic mechanisms behind flying organisms in nature are attracting attention due to their inspiring effects on the design and control of micro-aircraft. Unlike conventional fixed wing aircraft, the unsteady effects of micro aircraft are very complex and significant at low reynolds numbers due to their small size and low velocity. The motion parameters of the flapping wings can effectively influence the generation of the motion track and vortex of the flapping wings, and greatly influence the unsteady effect.
The traditional flapping wing motion parameter optimization method combines an optimization algorithm with different aerodynamic prediction models to optimize the motion parameters, and has the following defects:
a quasi-steady-state model or a constant vortex lattice method is used for the prediction model, and the precision of the model is low;
the prediction model is subjected to numerical simulation or experiment means, so that the cost is high;
disclosure of Invention
In order to overcome the defects of the prior art, the invention provides an efficient flapping wing motion parameter optimization method based on a proxy model, and the flapping wing motion parameter optimization under different target aerodynamic coefficients is effectively realized.
In order to achieve the above purpose, the solution of the invention is:
a method of optimizing flapping wing motion parameters, the method comprising:
(1) presetting sampling points: presetting sampling points in a parameter space formed by flapping wing motion parameters to form a sampling point set X1。
(2) Acquiring training data:
selecting a plurality of sets of calculation grids with different grid densities as grids with different accuracies, and calculating aerodynamic coefficients corresponding to the sampling points x preset in the step (1) by adopting a numerical simulation method respectively
Figure BDA0002355198490000012
Obtaining aerodynamic coefficients with different precisions to form a training sample set
Figure BDA0002355198490000011
t denotes the accuracy class, i denotes the different aerodynamic coefficients, X ∈ X1。
(3) Constructing a proxy model:
according to the obtained training sample set SiAnd constructing a proxy model through multi-precision Gaussian process regression. The input of the proxy model is a motion parameter, the output of the proxy model is an aerodynamic coefficient, and the multi-precision Gaussian process regression comprises two assumptions: all the aerodynamic coefficients follow Gaussian distribution, and the distribution of samples with different accuracies meets the following relation:
Zt(x)＝ρt-1(x,γ)Zt-1(x)+δt(x)
δt(x)⊥Zt-1(x)
Zt(x) Is composed ofOne of the aerodynamic coefficients C corresponding to all x in the t precision gradet(x) And (3) obeying a Gaussian distribution function, wherein rho (x, gamma) is a parameter function, gamma is a hyperparameter, and delta (x) is an independent Gaussian distribution. Gaussian distribution Z of lowest precision class samples1(x) The independent gaussian distribution formula is as follows:
Figure BDA0002355198490000021
Figure BDA0002355198490000022
fT(x) For regression functions, η and σ are hyper-parameters, r (x, x') is a covariance function, and commonly used covariance functions are gaussian functions or radial basis functions, etc.
Combining the Gaussian distributions of the samples with different precisions to construct a combined Gaussian distribution for writing
Figure BDA0002355198490000023
Figure BDA0002355198490000024
For a joint Gaussian distribution, mu(t)For the expectation of the joint Gaussian distribution, V(t)Is a covariance matrix of the joint gaussian distribution. The specific formula is as follows:
Figure BDA0002355198490000025
Figure BDA0002355198490000026
cov represents the corresponding covariance matrix for the a priori expectation of one of the aerodynamic coefficients with an accuracy rating t.
For gaussian process regression, the hyperparameters may be calculated by minimizing the negative log-edge likelihood function:
Figure BDA0002355198490000027
Figure BDA0002355198490000028
Θ is a hyper-parameter set, comprising γ, η and σ. NLML is a negative log-edge likelihood function.
A priori joint gaussian distribution of known multi-precision samples
Figure BDA0002355198490000029
For any motion parameter input x, probability posterior distribution can be calculated through Bayes theorem to obtain a proxy model, and then expected mu corresponding to the aerodynamic coefficient is obtainedpostAnd variance
Figure BDA00023551984900000210
Figure BDA00023551984900000211
Figure BDA00023551984900000212
Figure BDA00023551984900000213
And R is a motion parameter space. According to different aerodynamic coefficients
Figure BDA00023551984900000214
Establishing a corresponding agent model:
Figure BDA0002355198490000031
Figure BDA0002355198490000032
Figure BDA0002355198490000033
to obtain the corresponding expected mupost(i)And variance
Figure BDA0002355198490000034
Expectation of mupost(i)Namely the output aerodynamic coefficient C of the proxy model of the motion parameter xi'(x)。
Fourthly, optimizing flapping wing motion parameters: the method comprises the following steps of optimizing flapping wing motion parameters by taking the realization of a specific aerodynamic coefficient as an optimization target, calculating the relative error between the true aerodynamic coefficient of an optimal point and the optimization target in the optimization process, and judging whether the error requirement is met, wherein the method comprises the following steps:
(4.1) searching an optimal point and a candidate sampling point:
set of points X in the space of motion parameters2Finding the optimal point. The point set X2Consists of three parts: the method comprises the steps of a point set based on a grid search theory, a point set based on a random search theory and a point set based on a disturbance theory.
Based on an agent model, combining a Bayesian optimization theory, optimizing motion parameters to realize a specific aerodynamic coefficient, wherein the Bayesian optimization uses an expecteimprovement strategy as an acquisition function:
Figure BDA0002355198490000035
Fb(x)＝min(abs(Ci'(x)-<Ci>O(x) ) is the current minimum value.
Figure BDA0002355198490000036
abs represents the absolute value of the signal and,<Ci>O(x) Representing the target value of aerodynamic coefficient, i represents different aerodynamic coefficients, m is the number of aerodynamic coefficients, SiFor the training sample set, μi(x) Phi (M) is the absolute value of the difference between the aerodynamic coefficient of the motion parameter x and the target value, and sigma is the probability distribution functionpost(i)(x) Is the variance of the received signal and the received signal,
Figure BDA0002355198490000037
for the probability density function, use CμAnd CσTwo parameters control the weight of the expectation and variance, respectively, Cμ+Cσ＝1，Cμ≥0.5。
And selecting a point with the largest EI value as an optimal point, sequencing the sampling points according to the variance from large to small based on an active learning strategy, and selecting any previous sampling point as a candidate sampling point.
(4.2) optimizing convergence judgment;
and calculating the relative error between the true aerodynamic coefficient of the optimal point and the optimization target, and judging whether the error requirement is met. If the error requirement is not met, performing numerical simulation calculation on the aerodynamic coefficient of the optimal point by using the highest precision numerical value in the step two, performing numerical simulation calculation on the aerodynamic coefficients of the candidate sampling points by using the rest numerical values with different precision levels, adding the optimal point and the corresponding aerodynamic coefficient into the sample with the highest precision level, adding the selected sampling point and the aerodynamic coefficient into the sample with the corresponding precision level to form a new training sample set, repeating the steps (3) to (4), updating the modified proxy model and optimizing the motion parameter until the error requirement is met, completing the optimization of the motion parameter, and obtaining the motion parameter of the optimal point as the optimized motion parameter.
Further, the flapping wings move in the following mode:
y(t)/c＝hcos(2πft)
Figure BDA0002355198490000041
Figure BDA0002355198490000042
and under the coordinate system of the machine body, the motion form is composed of three-degree-of-freedom motion. y (t)/c is the pitching oscillating movement in the y direction, c is the wing chord length, h is the amplitude of the pitching movement, 0.05<h<1.25; f is the flapping frequency, 0.06<f<0.3. x (t) is oscillatory motion in the x direction, the amplitude of which is controlled by the amplitude of the pitching motion h and the stroke angle beta, 50 °<β<135 deg. Alpha (t) is an oscillatory rotation about a rotation axis o on the wing chord, alpha (t) being a sinusoidal-like rotation in the first half-cycle, alpha0To the amplitude of rotation, 0<α0<50 degrees; the value is constant at zero in the second half cycle. In the ground coordinate system, the flapping wing follows the fuselage at a constant velocity to the left. Alpha (t) is the included angle between the direction of the wing chord line and the motion direction, and theta (t) is the included angle between the wing chord line and the horizontal line. T is the moment corresponding to the instantaneous vorticity field, T is the flapping cycle of the flapping wing, and the value is the reciprocal of the flapping frequency. The motion parameters comprise: flapping frequency f, pitching motion amplitude h, stroke angle beta and rotation amplitude alpha0(ii) a The motion parameter space is composed of flapping wing motion parameters which meet the formula and the constraint condition.
Further, in the step (1), sampling points are preset in a parameter space formed by flapping wing motion parameters by using a Latin hypercube sampling method.
Further, in the step (2), the numerical simulation adopts a dynamic grid method or an immersion boundary method based on radial basis function interpolation.
Further, the optimization objective comprises a time-averaged lift coefficient and a time-averaged thrust coefficient; lift coefficient calculation formula
Figure BDA0002355198490000043
Thrust coefficient calculation formula
Figure BDA0002355198490000044
L and T are lift and thrust, ρ is fluid density, S is characteristic area, U∞Is the free incoming flow velocity.
The invention has the following beneficial effects: the invention effectively utilizes samples with different precisions to improve the precision of the proxy model; candidate sampling points are selected from the parameter space efficiently to improve the accuracy of the proxy model, so that the motion parameters of the flapping wings can be optimized by the proxy model efficiently; the method can obtain the optimal motion parameter under the specific aerodynamic coefficient, and the aerodynamic coefficient of the flapping wing aircraft is closer to the target aerodynamic coefficient under the motion parameter. The obtained optimal motion parameters have an auxiliary effect on the control aspect of the flapping wing aircraft, for example, the corresponding optimal motion parameters are obtained by taking the corresponding aerodynamic coefficients required by the hovering and hovering tasks of the aircraft as target aerodynamic coefficients according to the optimization method of the invention, and the aircraft is set to the corresponding optimal motion parameters, so that the aircraft can well execute the tasks.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a diagram of the motion pattern adopted by the flapping wings, wherein a is a coordinate system of the fuselage and b is a coordinate system of the ground;
FIG. 3 is a diagram of the optimized instantaneous flow field vorticity provided by the present invention;
fig. 4 is a diagram of the optimized flapping wing motion trail and stress condition provided by the invention.
Detailed Description
The present invention will be described in detail below with reference to the accompanying drawings.
FIG. 1 is a flow chart of the optimization of flapping wing motion parameters provided by the present application.
The motion parameter optimization method of the embodiment comprises the following steps:
s1: presetting sampling points in a parameter space of the flapping wing by using experimental design;
fig. 2 shows the motion form of the flapping wing of the present embodiment, but is not limited thereto. The mathematical expression of the motion form is as follows:
y(t)/c＝hcos(2πft)
Figure BDA0002355198490000051
Figure BDA0002355198490000052
and under the coordinate system of the machine body, the motion form is composed of three-degree-of-freedom motion. y (t)/c is the pitching oscillating movement in the y direction, c is the wing chord length, h is the amplitude of the pitching movement, 0.05<h<1.25; f is the flapping frequency, 0.06<f<0.3. x (t) is oscillatory motion in the x direction, the amplitude of which is controlled by the amplitude of the pitching motion h and the stroke angle beta, 50 °<β<135 deg. Alpha (t) is an oscillatory rotation about a rotation axis o on the wing chord, alpha (t) being a sinusoidal-like rotation in the first half-cycle, alpha0To the amplitude of rotation, 0<α0<50 degrees; the value is constant at zero in the second half cycle. In the ground coordinate system, the flapping wing follows the fuselage at a constant velocity to the left. Alpha (t) is the included angle between the direction of the wing chord line and the motion direction, and theta (t) is the included angle between the wing chord line and the horizontal line. T is the moment corresponding to the instantaneous vorticity field, T is the flapping cycle of the flapping wing, and the value is the reciprocal of the flapping frequency. The motion parameters comprise: flapping frequency f, pitching motion amplitude h, stroke angle beta and rotation amplitude alpha0(ii) a The motion parameter space is composed of flapping wing motion parameters which meet the formula and the constraint condition.
As an optional implementation mode of the invention, the sampling point experiment is designed to be Latin hypercube sampling. T groups of sampling points are obtained through t times of sampling to form a sampling point set X1And the number of sampling points in each group is sequentially reduced, and t is the highest precision level of the multi-precision proxy model. In this embodiment, t is 2.
S2: and calculating the aerodynamic coefficient corresponding to the sampling point by using a numerical simulation method.
The flapping wing flow field can be numerically simulated by using a dynamic grid method or an immersion boundary method based on radial basis function interpolation, in the embodiment, the dynamic grid method based on open source software OpenFOAM is adopted for simulation, grid independence verification is carried out, appropriate t grids with different accuracies are selected from the dynamic grid method and are respectively used as grids with high accuracy and low accuracy, and aerodynamic coefficients of sampling points are calculated according to lift force, thrust force and the like obtained through numerical simulation and are used as training samples Si(x,
Figure BDA0002355198490000066
) I denotes the different aerodynamic coefficients, X ∈ X1. As an alternative implementation manner of the present invention, the aerodynamic coefficient includes a lift coefficient and a thrust coefficient. Lift coefficient calculation formula
Figure BDA0002355198490000061
Thrust coefficient calculation formula
Figure BDA0002355198490000062
L and T are lift and thrust, ρ is fluid density, S is characteristic area, U∞Is the free incoming flow velocity.
It should be noted that, the higher the accuracy level of the numerical simulation, the lower the number of sampling points corresponding to the sampling point group.
S3: a proxy model is constructed using the multi-precision samples.
The method for constructing the proxy model is multi-precision Gaussian process regression. The input of the proxy model is a motion parameter, and the output is an aerodynamic coefficient. The multiple precision gaussian process regression includes two important assumptions. All aerodynamic coefficients follow a gaussian distribution, and the distribution of samples of different accuracies satisfies the following relationship:
Zt(x)＝ρt-1(x,γ)Zt-1(x)+δt(x)
δt(x)⊥Zt-1(x)
Zt(x) One of high-precision samples
Figure BDA0002355198490000067
Of a Gaussian distribution function of, Zt-1(x) Is one of low-precision samples
Figure BDA0002355198490000063
The gaussian distribution function of (1) is marked with different precision levels, ρ (x, γ) is a parameter function, and the embodiment adopts a constant function, and takes γ as a hyperparameter, and δ (x) is an independent gaussian distribution. Gaussian distribution Z of lowest precision class samples1(x) With said independent Gaussian fractionThe cloth formula is as follows:
Figure BDA0002355198490000064
Figure BDA0002355198490000065
fT(x) In this embodiment, the values are set to zero, η and σ are hyper-parameters, r (x, x') is a covariance function, and the common covariance function is a gaussian function or a radial basis function.
Combining the Gaussian distributions of the samples with different precisions to construct a combined Gaussian distribution for writing
Figure BDA0002355198490000071
Figure BDA0002355198490000072
With combined Gaussian distribution, mu, of high and low precision(t)For the expectation of the joint Gaussian distribution, V(t)Is a covariance matrix of the joint gaussian distribution.
The specific formula is as follows:
Figure BDA0002355198490000073
Figure BDA0002355198490000074
cov represents the corresponding covariance matrix for the a priori expectation of one of the aerodynamic coefficients with an accuracy rating t.
For gaussian process regression, the hyperparameters may be calculated by minimizing the negative log-edge likelihood function.
Figure BDA0002355198490000075
Figure BDA0002355198490000076
Θ is a hyper-parameter set, comprising γ, η and σ. NLML is a negative log-edge likelihood function.
A priori joint gaussian distribution of known multi-precision samples
Figure BDA0002355198490000077
For any motion parameter input x, probability posterior distribution can be calculated through Bayes theorem to obtain a proxy model, and then expected mu corresponding to the aerodynamic coefficient is obtainedpostAnd variance
Figure BDA0002355198490000078
Figure BDA0002355198490000079
Figure BDA00023551984900000710
Figure BDA00023551984900000711
And R is a motion parameter space. Establishing corresponding proxy models according to different aerodynamic coefficients:
Figure BDA00023551984900000712
Figure BDA00023551984900000713
Figure BDA00023551984900000714
to obtain the corresponding expected mupost(i)And variance
Figure BDA00023551984900000715
Expectation of mupost(i)Namely the output aerodynamic coefficient C of the proxy model of the motion parameter xi'(x)。
S4: and performing multi-objective aerodynamic coefficient optimization on the flapping wings, wherein the optimization problem can be expressed as:
Figure BDA0002355198490000081
abs represents the absolute value of the signal and,<Ci>O(x) Representing the target value of the aerodynamic coefficient, according to the stress requirement of the aircraft, and the value of m is 2.
In the optimization process, calculating the relative error between the true aerodynamic coefficient of the optimal point and the optimization target, and judging whether the error requirement is met, wherein the method comprises the following steps:
(4.1) searching an optimal point and a candidate sampling point:
set of points X in parameter space2Finding the optimal point. The point set X2Is composed of three parts
A point set based on a grid search theory, wherein points in the point set are uniformly distributed in a parameter space;
a point set based on a random search theory, wherein points in the point set are randomly distributed in a parameter space;
and based on the point set of the perturbation theory, the points in the point set are randomly distributed near the optimal point in the parameter space.
Based on a proxy model, combining a Bayes optimization theory, optimizing motion parameters to realize specific aerodynamic coefficients, wherein the optimized aerodynamic coefficients comprise time-averaged lift coefficients and time-averaged thrust coefficients. The Bayesian optimization uses an expectedprovement strategy as an acquisition function:
Figure BDA0002355198490000082
Fb(x)＝min(abs(Ci'(x)-<Ci>O(x) ) is the current minimum value.
Figure BDA0002355198490000083
SiFor the training sample set, μi(x) Phi (M) is the absolute value of the difference between the aerodynamic coefficient of the motion parameter x and the target value, and sigma is the probability distribution functionpost(i)(x) Is the variance of the received signal and the received signal,
Figure BDA0002355198490000084
for the probability density function, use CμAnd CσTwo parameters control the weight of expectation and variance respectively, in this embodiment, Cμ＝0.5，Cσ＝0.5。
And selecting a point with the largest EI value as an optimal point, sequencing the sampling points according to the variance from large to small based on an active learning strategy, and selecting any previous sampling point as a candidate sampling point.
(4.2) optimizing convergence judgment; and (3) calculating the real aerodynamic coefficient of the optimal point by using high-precision numerical simulation, and judging whether the relative error between the real value and the target aerodynamic coefficient meets the error requirement, wherein the error requirement is that the relative error is less than 10%.
And (3) performing numerical simulation calculation on the aerodynamic coefficient of the optimal point by using the highest precision numerical value in the step two, performing numerical simulation calculation on candidate sampling point aerodynamic coefficients by using the rest numerical values with different precision grades, adding the optimal point and the corresponding aerodynamic coefficient into the sample with the highest precision grade, adding the selected sampling point and the aerodynamic coefficient into the sample with the corresponding precision grade to form a new training sample set, repeating the steps (3) to (4), updating and correcting the proxy model and optimizing the motion parameters until the error requirements are met, and completing the optimization of the motion parameters.
The optimization results of this embodiment, in which the time-averaged lift coefficient is 3 and the time-averaged thrust coefficient is 0 as the optimization target, are shown in the following table:
table 1: optimization results table
Number of iterations 1
Optimum flapping frequency (Hz) 0.1660
Optimal amplitude of pitching motion 1.1166
Optimum stroke angle (rad) 0.9972
Optimum rotation amplitude (rad) 0.6271
Time-average lift coefficient 2.8694
Time-averaged thrust coefficient -0.0126
It can be seen that after one iteration, the aerodynamic coefficient at the optimal point already meets the error requirement. Fig. 3-4 show the optimization results of the simulation of this embodiment, from which the instantaneous flow field vortex diagram of the flapping wing after optimization, the trajectory diagram of the flapping wing in one cycle, and the aerodynamic vector arrows can be seen. Under the motion parameter, the aerodynamic coefficient of the flapping wing aircraft is closer to the target aerodynamic coefficient. The motion parameters of the flapping wing air vehicle are set to the numerical values in the table, so that the air vehicle can well execute the hovering task. The method effectively utilizes the samples with different precisions to improve the precision of the proxy model, and simultaneously selects the candidate sampling points in the parameter space efficiently to improve the precision of the proxy model, thereby optimizing the flapping wing motion parameters by efficiently utilizing the proxy model. The obtained optimal motion parameters have an auxiliary effect on the control of the flapping wing aircraft.

Claims (5)

1. A method for optimizing flapping wing motion parameters, the method comprising:
(1) presetting sampling points: presetting sampling points in a parameter space formed by flapping wing motion parameters to form a sampling point set X1；
(2) Acquiring training data:
selecting a plurality of sets of calculation grids with different grid densities as grids with different precisions, and calculating aerodynamic coefficients corresponding to sampling points x preset in the step (1) by adopting a numerical simulation method
Figure FDA0003300758960000016
Obtaining aerodynamic coefficients with different precisions to form a training sample set
Figure FDA0003300758960000017
Figure FDA0003300758960000018
t denotes the accuracy class, i denotes the different aerodynamic coefficients, X ∈ X1；
(3) Constructing a proxy model:
according to the obtained training sample set SiConstructing a proxy model through multi-precision Gaussian process regression; the input of the proxy model is a motion parameter, the output of the proxy model is an aerodynamic coefficient, and the multi-precision Gaussian process regression comprises two assumptions: all the aerodynamic coefficients follow Gaussian distribution, and the distribution of samples with different accuracies meets the following relation:
Zt(x)＝ρt-1(x,γ)Zt-1(x)+δt(x)
δt(x)⊥Zt-1(x)
Zt(x) One of the aerodynamic coefficients C corresponding to all x in the t precision gradet(x) The obeyed Gaussian distribution function, rho (x, gamma) is a parameter function, gamma is a hyper-parameter, and delta (x) is independent Gaussian distribution; gaussian distribution Z of lowest precision class samples1(x) The independent gaussian distribution formula is as follows:
Figure FDA0003300758960000011
Figure FDA0003300758960000012
fT(x) Is a regression function, eta and sigma are hyper-parameters, r (x, x') is a covariance function, and the covariance function is a Gaussian function or a radial basis function;
combining the Gaussian distributions of the samples with different precisions to construct a combined Gaussian distribution for writing
Figure FDA0003300758960000014
Figure FDA0003300758960000015
For a joint Gaussian distribution, mu(t)For the expectation of the joint Gaussian distribution, V(t)A covariance matrix which is a joint gaussian distribution; the specific formula is as follows:
Figure FDA0003300758960000013
Figure FDA0003300758960000028
cov represents a corresponding covariance matrix, which is a priori expectation of one of the aerodynamic coefficients with the accuracy grade t;
for gaussian process regression, the hyperparameters may be calculated by minimizing the negative log-edge likelihood function:
Figure FDA0003300758960000021
Figure FDA0003300758960000022
Θ is a hyper-parameter set, including γ, η, and σ; NLML is a negative log-edge likelihood function;
a priori joint gaussian distribution of known multi-precision samples
Figure FDA0003300758960000027
For any motion parameter input x, probability posterior distribution can be calculated through Bayes theorem to obtain a proxy model, and expected mu corresponding to the aerodynamic coefficient is obtainedpostAnd variance
Figure FDA0003300758960000023
Figure FDA0003300758960000024
R is a motion parameter space; according to different aerodynamic coefficients
Figure FDA0003300758960000029
Establishing a corresponding agent model:
Figure FDA0003300758960000025
to obtain the corresponding expected mupost(i)And variance
Figure FDA0003300758960000026
Expectation of mupost(i)Is the generation of the motion parameter xOutput aerodynamic coefficient C of physical modeli'(x)；
(4) Optimizing the motion parameters of the flapping wings: the method comprises the following steps of optimizing flapping wing motion parameters by taking the realization of a specific aerodynamic coefficient as an optimization target, calculating the relative error between the true aerodynamic coefficient of an optimal point and the optimization target in the optimization process, and judging whether the error requirement is met, wherein the method comprises the following steps:
(4.1) searching an optimal point and a candidate sampling point:
set of points X in the space of motion parameters2Searching for an optimal point; the point set X2Consists of three parts: a point set based on a grid search theory, a point set based on a random search theory and a point set based on a disturbance theory;
based on an agent model, combining a Bayesian optimization theory, optimizing motion parameters to realize a specific aerodynamic coefficient, wherein the Bayesian optimization uses an expecteimprovement strategy as an acquisition function:
Figure FDA0003300758960000031
Fb(x)＝min(abs(Ci'(x)-<Ci>O(x) ) is the current minimum value;
Figure FDA0003300758960000032
abs represents the absolute value of the signal and,<Ci>O(x) Representing the target value of aerodynamic coefficient, i represents different aerodynamic coefficients, m is the number of aerodynamic coefficients, SiFor the training sample set, μi(x) Phi (M) is the absolute value of the difference between the aerodynamic coefficient of the motion parameter x and the target value, and sigma is the probability distribution functionpost(i)(x) Is the variance of the received signal and the received signal,
Figure FDA0003300758960000035
for the probability density function, use CμAnd CσTwo parameters control the expectation and the square respectivelyWeight of difference, Cμ+Cσ＝1，Cμ≥0.5；
Selecting a point with the largest EI value as an optimal point, sequencing sampling points according to variance from large to small based on an active learning strategy, and selecting any previous sampling point as a candidate sampling point;
(4.2) optimizing convergence judgment;
calculating the relative error between the true aerodynamic coefficient of the optimal point and the optimization target, and judging whether the error requirement is met; if the error requirement is not met, performing numerical simulation calculation on the aerodynamic coefficient of the optimal point by using the highest precision numerical value in the step two, performing numerical simulation calculation on the aerodynamic coefficients of the candidate sampling points by using the rest numerical values with different precision levels, adding the optimal point and the corresponding aerodynamic coefficient into the sample with the highest precision level, adding the candidate sampling point and the aerodynamic coefficient into the sample with the corresponding precision level to form a new training sample set, repeating the steps (3) to (4), updating and correcting the proxy model and optimizing the motion parameter until the error requirement is met, completing the optimization of the motion parameter, and taking the motion parameter of the optimal point as the optimized motion parameter.
2. The method of claim 1, wherein the flapping wings move in the form of:
y(t)/c＝hcos(2πft)
Figure FDA0003300758960000033
Figure FDA0003300758960000034
wherein, under the coordinate system of the fuselage, the movement form is composed of three-degree-of-freedom movement; y (t)/c is the pitching oscillating movement in the y direction, c is the wing chord length, h is the amplitude of the pitching movement, 0.05<h<1.25; f is the flapping frequency, 0.06<f<0.3; x (t) is oscillatory motion in the x direction, the amplitude of which is controlled by the amplitude of the pitching motion h and the stroke angle beta, 50 °<β<135°；α(t) is an oscillatory rotation about an axis of rotation o on the wing chord, alpha (t) being a sine-like rotation in the first half-cycle, alpha0To the amplitude of rotation, 0<α0<50 degrees; the value is constantly zero in the latter half period; in a ground coordinate system, the flapping wing moves leftwards at a constant speed along with the fuselage; alpha (t) is an included angle between the direction of the wing chord line and the motion direction, and theta (t) is an included angle between the wing chord line and a horizontal line; t is the moment corresponding to the instantaneous vorticity field, T is the flapping cycle of the flapping wing, and the value is the reciprocal of the flapping frequency; the motion parameters comprise: flapping frequency f, pitching motion amplitude h, stroke angle beta and rotation amplitude alpha0(ii) a The motion parameter space is composed of flapping wing motion parameters which meet the formula and the constraint condition.
3. The method for optimizing kinetic parameters of claim 1 wherein in step (1), the sampling points are preset in the parameter space formed by the flapping-wing kinetic parameters using the Latin hypercube sampling method.
4. The method of claim 1, wherein in the step (2), the numerical simulation employs a dynamic grid method based on radial basis function interpolation or an immersion boundary method.
5. The method of claim 1, wherein the optimization objectives include time-averaged lift coefficients and time-averaged thrust coefficients; lift coefficient calculation formula
Figure FDA0003300758960000041
Thrust coefficient calculation formula
Figure FDA0003300758960000042
L and T are lift and thrust, ρ is fluid density, S is characteristic area, U∞Is the free incoming flow velocity.
CN202010005675.0A 2020-01-03 2020-01-03 Flapping wing motion parameter optimization method Active CN111199105B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010005675.0A CN111199105B (en) 2020-01-03 2020-01-03 Flapping wing motion parameter optimization method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010005675.0A CN111199105B (en) 2020-01-03 2020-01-03 Flapping wing motion parameter optimization method

Publications (2)

Publication Number Publication Date
CN111199105A CN111199105A (en) 2020-05-26
CN111199105B true CN111199105B (en) 2022-03-22

Family

ID=70746788

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010005675.0A Active CN111199105B (en) 2020-01-03 2020-01-03 Flapping wing motion parameter optimization method

Country Status (1)

Country Link
CN (1) CN111199105B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112407139B (en) * 2020-11-14 2023-05-23 西北工业大学 Active drag reduction method for flapping-wing wake flow control of underwater vehicle
CN116502566A (en) * 2023-06-27 2023-07-28 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) Multi-objective optimization method for performance of combustion chamber of gas turbine based on Bayesian optimization
CN116522068B (en) * 2023-07-03 2023-09-15 西安羚控电子科技有限公司 Test parameter generation method and system

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103699757A (en) * 2014-01-06 2014-04-02 西北工业大学 Micro flapping-wing analysis system and method involved with pneumatic and structural coupling properties
CN104568373A (en) * 2014-12-20 2015-04-29 浙江大学 Testing device and testing method for mass force of minitype ornithopter
CN107729639A (en) * 2017-10-10 2018-02-23 东莞理工学院 A kind of outstanding design method for flying lower wing of imitative hummingbird flapping-wing modal
CN109783978A (en) * 2019-02-11 2019-05-21 南京航空航天大学 A kind of micro flapping wing air vehicle aerodynamic numerical emulation mode based on ANSYS Workbench
CN110008639A (en) * 2019-04-24 2019-07-12 东莞理工学院 A kind of micro flapping wing air vehicle wing intelligent parameter design method
CN110334401A (en) * 2019-05-31 2019-10-15 南京航空航天大学 A kind of double flapping wings based on tandem arrangement promote the optimization method of efficiency

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2001015971A2 (en) * 1999-08-30 2001-03-08 Smith Michael J C Wing-drive mechanism and vehicle employing same
JP4115349B2 (en) * 2002-07-12 2008-07-09 シャープ株式会社 Flapping levitation moving device
FR2944896B1 (en) * 2009-04-24 2011-07-08 Airbus France METHOD FOR PREDICTING AERODYNAMIC BEHAVIOR OF A STRUCTURE OF AN AIRCRAFT
CN109446688A (en) * 2018-11-07 2019-03-08 上海海事大学 One kind is based on two-dimensional surface hovering flapping wing Aerodynamic characteristics method
CN110254742A (en) * 2019-07-02 2019-09-20 深圳灵动牛科技有限责任公司 A kind of Wing design method of flapping-wing type aircraft

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103699757A (en) * 2014-01-06 2014-04-02 西北工业大学 Micro flapping-wing analysis system and method involved with pneumatic and structural coupling properties
CN104568373A (en) * 2014-12-20 2015-04-29 浙江大学 Testing device and testing method for mass force of minitype ornithopter
CN107729639A (en) * 2017-10-10 2018-02-23 东莞理工学院 A kind of outstanding design method for flying lower wing of imitative hummingbird flapping-wing modal
CN109783978A (en) * 2019-02-11 2019-05-21 南京航空航天大学 A kind of micro flapping wing air vehicle aerodynamic numerical emulation mode based on ANSYS Workbench
CN110008639A (en) * 2019-04-24 2019-07-12 东莞理工学院 A kind of micro flapping wing air vehicle wing intelligent parameter design method
CN110334401A (en) * 2019-05-31 2019-10-15 南京航空航天大学 A kind of double flapping wings based on tandem arrangement promote the optimization method of efficiency

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于双BP神经网络的扑翼飞行器气动参数辨识;韩建福等;《计算机应用》;20191230;229-302 *

Also Published As

Publication number Publication date
CN111199105A (en) 2020-05-26

Similar Documents

Publication Publication Date Title
CN111199105B (en) Flapping wing motion parameter optimization method
Jeong et al. Data mining for aerodynamic design space
Han et al. Online policy iteration ADP-based attitude-tracking control for hypersonic vehicles
Qu et al. An improved TLBO based memetic algorithm for aerodynamic shape optimization
Wickramasinghe et al. Designing airfoils using a reference point based evolutionary many-objective particle swarm optimization algorithm
CN108090621B (en) Short-term wind speed prediction method and system based on staged overall optimization
CN108171315B (en) Multi-unmanned aerial vehicle task allocation method based on SMC particle swarm algorithm
CN110516318B (en) Airfoil design method based on radial basis function neural network proxy model
CN110705029B (en) Flow field prediction method of oscillating flapping wing energy acquisition system based on transfer learning
CN110543727B (en) Improved particle swarm algorithm-based omnidirectional mobile intelligent wheelchair robot parameter identification method
CN114355777B (en) Dynamic gliding method and system based on distributed pressure sensor and sectional attitude control
CN111581784B (en) Flapping wing motion parameter optimization method based on data-driven self-adaptive quasi-steady-state model
Chen et al. Optimization of flatback airfoils for wind-turbine blades using a genetic algorithm
CN113777931A (en) Icing wing type pneumatic model construction method, device, equipment and medium
CN108427428A (en) Based on the adaptive sliding moding structure Spacecraft Attitude Control method for improving iterative algorithm
Liu et al. Optimization of nano-rotor blade airfoil using controlled elitist NSGA-II
CN114818477A (en) Aircraft aerodynamic characteristic prediction method, device and medium
CN110765706A (en) Airfoil unsteady stall aerodynamic coefficient modeling method based on OHNGBM (1,1)
Bouhelal et al. Blade Element Momentum Theory Coupled with Machine Learning to Predict Wind Turbine Aerodynamic Performances
CN112016162A (en) Four-rotor unmanned aerial vehicle PID controller parameter optimization method
Mukesh et al. Influence of search algorithms on aerodynamic design optimisation of aircraft wings
CN114462149B (en) Aircraft pneumatic parameter identification method based on pre-training and incremental learning
Trad et al. Airfoils generation using neural networks, CST curves and aerodynamic coefficients
CN115963855A (en) Unpowered reentry aircraft landing area prediction method based on deep learning
Xiang et al. A Manifold-Based Airfoil Geometric-Feature Extraction and Discrepant Data Fusion Learning Method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant