[go: up one dir, main page]

CN111199105B - An optimization method for flapping motion parameters - Google Patents

An optimization method for flapping motion parameters 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
motion
aerodynamic
precision
optimization
flapping
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

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及计算机仿真技术领域，尤其涉及一种扑翼运动参数优化方法。该方法包括：使用实验设计在参数空间内预设采样点，使用数值模拟方法计算采样点对应的气动力系数，使用高精度与低精度样本构建代理模型，结合贝叶斯优化对所述扑翼进行运动参数优化，优化收敛判断。当最优点的真实气动力系数不满足误差要求时，将最优点计算结果作为高精度样本，使用主动学习策略选取候选采样点，通过低精度数值模拟计算对应气动力系数作为低精度样本，更新修正代理模型。本发明有效利用了不同精度的数据样本。整合优化流程，根据目标气动力系数，科学高效地对扑翼的运动参数进行了优化，以实现特定的气动性能。

Figure 202010005675

The invention relates to the technical field of computer simulation, in particular to a method for optimizing flapping motion parameters. The method includes: using an experimental design to preset sampling points in a parameter space, using a numerical simulation method to calculate the aerodynamic coefficients corresponding to the sampling points, using high-precision and low-precision samples to construct a surrogate model, and combining Bayesian optimization for the flapping wing Perform motion parameter optimization and optimize convergence judgment. When the true aerodynamic coefficient of the optimal point does not meet the error requirements, the calculation result of the optimal point is used as a high-precision sample, and the active learning strategy is used to select candidate sampling points, and the corresponding aerodynamic coefficient is calculated by low-precision numerical simulation as a low-precision sample, and the update correction proxy model. The present invention effectively utilizes data samples of different precisions. Integrating the optimization process, according to the target aerodynamic coefficient, the motion parameters of the flapping wing are scientifically and efficiently optimized to achieve specific aerodynamic performance.

Figure 202010005675

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.一种扑翼运动参数优化方法，其特征在于，该方法包括：1. a flapping motion parameter optimization method, is characterized in that, the method comprises: (1)采样点预设：在扑翼运动参数构成的参数空间中预设采样点，组成采样点集X1；(1) Sampling point preset: preset sampling points in the parameter space formed by flapping motion parameters to form a sampling point set X 1 ; (2)训练数据获取：(2) Training data acquisition: 选取多套网格密度不同的计算网格作为不同精度网格，采用数值模拟方法计算步骤(1)预设的采样点x对应的气动力系数
Figure FDA0003300758960000016
获得不同精度的气动力系数组成训练样本集
Figure FDA0003300758960000017
Figure FDA0003300758960000018
t表示精度等级，i表示不同气动力系数，x∈X1；Select multiple sets of calculation grids with different grid densities as grids of different precision, and use the numerical simulation method to calculate the aerodynamic coefficient corresponding to the preset sampling point x in step (1).
Figure FDA0003300758960000016
Obtain aerodynamic coefficients of different precisions to form a training sample set
Figure FDA0003300758960000017
Figure FDA0003300758960000018
t represents the accuracy level, i represents different aerodynamic coefficients, x∈X 1 ; (3)代理模型构建：(3) Proxy model construction: 根据获取的训练样本集Si通过多精度高斯过程回归构建代理模型；其中，代理模型的输入为运动参数，输出为气动力系数，所述多精度高斯过程回归包括两个假设：所有气动力系数服从高斯分布、不同精度样本的分布满足如下关系：According to the obtained training sample set S i , a surrogate model is constructed through multi-precision Gaussian process regression; the input of the surrogate model is motion parameters, and the output is aerodynamic coefficients. The multi-precision Gaussian process regression includes two assumptions: all aerodynamic coefficients The distribution of samples with different precisions obeying a Gaussian distribution satisfies the following relationship: Zt(x)＝ρt-1(x,γ)Zt-1(x)+δt(x)Z t (x)=ρ t-1 (x,γ)Z t-1 (x)+δ t (x) δt(x)⊥Zt-1(x)δ t (x)⊥Z t-1 (x) Zt(x)为t精度等级里所有x对应的其中一种气动力系数Ct(x)服从的高斯分布函数，ρ(x,γ)为参数函数，γ为超参数，δ(x)为独立的高斯分布；最低精度等级样本的高斯分布Z1(x)与所述独立高斯分布公式如下：Z t (x) is the Gaussian distribution function of one of the aerodynamic coefficients C t (x) corresponding to all x in the t precision level, ρ(x, γ) is the parameter function, γ is the hyperparameter, δ(x) is an independent Gaussian distribution; the Gaussian distribution Z 1 (x) of the lowest precision level sample and the independent Gaussian distribution formula are as follows:
Figure FDA0003300758960000011
Figure FDA0003300758960000011
Figure FDA0003300758960000012
Figure FDA0003300758960000012
 fT(x)为回归函数，η与σ为超参数，r(x,x′)为协方差函数，协方差函数为高斯函数或径向基函数；f T (x) is a regression function, η and σ are hyperparameters, r(x, x′) is a covariance function, and the covariance function is a Gaussian function or a radial basis function; 结合不同精度样本的高斯分布可构造联合高斯分布，写作
Figure FDA0003300758960000014
Figure FDA0003300758960000015
为联合高斯分布，μ(t)为联合高斯分布的期望，V(t)为联合高斯分布的协方差矩阵；具体公式如下：Combining the Gaussian distributions of different precision samples can construct a joint Gaussian distribution, which is written as
Figure FDA0003300758960000014
Figure FDA0003300758960000015
is the joint Gaussian distribution, μ (t) is the expectation of the joint Gaussian distribution, and V (t) is the covariance matrix of the joint Gaussian distribution; the specific formula is as follows:
Figure FDA0003300758960000013
Figure FDA0003300758960000013
Figure FDA0003300758960000028
为精度等级为t的其中一种气动力系数的先验期望，Cov代表对应协方差矩阵；
Figure FDA0003300758960000028
is the prior expectation of one of the aerodynamic coefficients with an accuracy level of t, and Cov represents the corresponding covariance matrix; 对于高斯过程回归，所述的超参数可以通过最小化负对数边缘似然函数计算：For Gaussian process regression, the hyperparameters described can be computed by minimizing the negative log marginal likelihood function:
Figure FDA0003300758960000021
Figure FDA0003300758960000021
Figure FDA0003300758960000022
Figure FDA0003300758960000022
 Θ为超参数集，包括γ，η与σ；NLML为负对数边缘似然函数；Θ is the hyperparameter set, including γ, η and σ; NLML is the negative logarithmic edge likelihood function; 已知多精度样本的先验联合高斯分布
Figure FDA0003300758960000027
对于任意的运动参数输入x,可通过贝叶斯定理计算概率后验分布，得到代理模型，获得对应气动力系数的期望μpost与方差
Figure FDA0003300758960000023
Prior joint Gaussian distribution for known multi-precision samples
Figure FDA0003300758960000027
For any motion parameter input x, the probability posterior distribution can be calculated by Bayes' theorem, the surrogate model can be obtained, and the expected μ post and variance of the corresponding aerodynamic coefficient can be obtained
Figure FDA0003300758960000023
Figure FDA0003300758960000024
Figure FDA0003300758960000024
 R为运动参数空间；根据不同气动力系数
Figure FDA0003300758960000029
建立对应的代理模型：R is the motion parameter space; according to different aerodynamic coefficients
Figure FDA0003300758960000029
Create the corresponding proxy model:
Figure FDA0003300758960000025
Figure FDA0003300758960000025
 得到相应的期望μpost(i)与方差
Figure FDA0003300758960000026
期望μpost(i)即为该运动参数x的代理模型输出气动力系数Ci'(x)；get the corresponding expected μpost(i) and variance
Figure FDA0003300758960000026
The expectation μ post(i) is the output aerodynamic coefficient C i '(x) of the proxy model of the motion parameter x; (4)扑翼运动参数优化：以实现特定的气动力系数为优化目标，对扑翼运动参数进行优化，在优化过程中，计算最优点真实气动力系数与优化目标的相对误差，判断是否满足误差要求，包括如下步骤：(4) Optimization of flapping motion parameters: To achieve a specific aerodynamic coefficient as the optimization goal, optimize the flapping motion parameters. During the optimization process, calculate the relative error between the real aerodynamic coefficient of the optimal point and the optimization target, and judge whether it is satisfied Error requirements, including the following steps: (4.1)寻找最优点和候选采样点：(4.1) Find the optimal point and candidate sampling points: 在运动参数空间内的点集X2中寻找最优点；所述点集X2由三个部分组成：基于网格搜索理论的点集、基于随机搜索理论的点集、基于扰动理论的点集；Find the optimal point in the point set X 2 in the motion parameter space; the point set X 2 consists of three parts: the point set based on grid search theory, the point set based on random search theory, and the point set based on perturbation theory ; 基于代理模型，结合贝叶斯优化理论，优化运动参数实现特定气动力系数，所述的贝叶斯优化使用Expected improvement策略为采集函数：Based on the surrogate model, combined with the Bayesian optimization theory, the motion parameters are optimized to achieve specific aerodynamic coefficients. The Bayesian optimization uses the Expected improvement strategy as the acquisition function:
Figure FDA0003300758960000031
Figure FDA0003300758960000031
 Fb(x)＝min(abs(Ci'(x)-<Ci>O(x)))，为当前最小值；F b (x)=min(abs(C i '(x)-<C i > O (x))), which is the current minimum value;
Figure FDA0003300758960000032
Figure FDA0003300758960000032
 abs代表绝对值，<Ci>O(x)代表气动力系数目标值，i表示不同气动力系数,m值为气动力系数的个数，Si为训练样本集，μi(x)为该运动参数x的气动力系数与目标值之差的绝对值，Φ(M)为概率分布函数，σpost(i)(x)为方差，
Figure FDA0003300758960000035
为概率密度函数，使用Cμ与Cσ两个参数分别控制期望与方差的权重，Cμ+Cσ＝1，Cμ≥0.5；abs represents the absolute value, <C i > O (x) represents the target value of aerodynamic coefficients, i represents different aerodynamic coefficients, m is the number of aerodynamic coefficients, S i is the training sample set, and μ i (x) is The absolute value of the difference between the aerodynamic coefficient of the motion parameter x and the target value, Φ(M) is the probability distribution function, σ post(i) (x) is the variance,
Figure FDA0003300758960000035
is a probability density function, using C μ and C σ to control the weight of expectation and variance respectively, C μ +C σ =1, C μ ≥0.5; 选取EI值最大的点作为最优点，并基于主动学习策略，将采样点按方差从大到小排序，选取前任意个采样点作为候选采样点；The point with the largest EI value is selected as the optimal point, and based on the active learning strategy, the sampling points are sorted according to the variance from large to small, and any previous sampling point is selected as the candidate sampling point; (4.2)优化收敛判断；(4.2) Optimization convergence judgment; 计算最优点真实气动力系数与优化目标的相对误差，判断是否满足误差要求；若不满足误差要求，使用步骤二中最高精度数值模拟计算最优点的气动力系数，剩余不同精度等级数值模拟计算候选采样点气动力系数，最优点和对应的气动力系数加入最高精度等级的样本中，候选采样点及气动力系数加入对应的精度等级的样本中组成新的训练样本集，重复步骤(3)-(4)，更新修正代理模型和优化运动参数，直至满足误差要求，运动参数优化完成，最优点的运动参数即为优化后的运动参数。Calculate the relative error between the real aerodynamic coefficient of the optimal point and the optimization target, and judge whether the error requirements are met; if the error requirements are not met, use the highest precision numerical simulation in step 2 to calculate the aerodynamic coefficient of the optimal point, and the remaining candidates for numerical simulation calculation of different precision levels The aerodynamic coefficient of the sampling point, the optimal point and the corresponding aerodynamic coefficient are added to the samples of the highest accuracy level, and the candidate sampling points and aerodynamic coefficients are added to the samples of the corresponding accuracy level to form a new training sample set, and repeat step (3)- (4), update and correct the surrogate model and optimize the motion parameters until the error requirements are met, the motion parameter optimization is completed, and the motion parameters at the optimal point are the optimized motion parameters. 2.根据权利要求1所述的运动参数优化方法，其特征在于，所述扑翼的运动形式如下：2. motion parameter optimization method according to claim 1, is characterized in that, the motion form of described flapping wing is as follows: y(t)/c＝hcos(2πft)y(t)/c=hcos(2πft)
Figure FDA0003300758960000033
Figure FDA0003300758960000033
Figure FDA0003300758960000034
Figure FDA0003300758960000034
 其中，在机身坐标系下，所述运动形式由三自由度运动组成；y(t)/c为y方向上的俯仰振荡运动，c为机翼弦长，h为俯仰运动的幅值，0.05<h<1.25；f为扑动频率，0.06<f<0.3；x(t)为x方向上的振荡运动，其幅值由俯仰运动幅值h与冲程角β控制，50°<β<135°；α(t)为绕机翼弦上一旋转轴o的振荡转动,α(t)在前半个周期做类正弦转动，α0为转动幅值，0<α0<50°；在后半个周期值恒为零；在地面坐标系中，扑翼跟随机身以恒定速度向左运动；α(t)为机翼弦线方向与运动方向的夹角，θ(t)为机翼弦线与水平线的夹角；t为瞬时涡量场对应的时刻，T为扑翼的扑动周期，值为扑动频率的倒数；其运动参数包括：扑动频率f、俯仰运动幅值h、冲程角β、转动幅值α0；运动参数空间由满足上述公式及约束条件的扑翼运动参数构成。Among them, in the fuselage coordinate system, the motion form is composed of three degrees of freedom motion; y(t)/c is the pitch oscillation motion in the y direction, c is the wing chord length, h is the amplitude of the pitch motion, 0.05<h<1.25; f is the flutter frequency, 0.06<f<0.3; x(t) is the oscillating motion in the x direction, and its amplitude is controlled by the pitch motion amplitude h and the stroke angle β, 50°<β<135°; α(t) is the oscillating rotation around the rotation axis o of the wing chord, α(t) does a quasi-sinusoidal rotation in the first half cycle, α 0 is the rotation amplitude, 0<α 0 <50°; The value of the second half period is always zero; in the ground coordinate system, the flapping wing follows the fuselage and moves to the left at a constant speed; α(t) is the angle between the direction of the wing chord and the direction of movement, and θ(t) is the The angle between the chord line and the horizontal line; t is the moment corresponding to the instantaneous vorticity field, T is the flapping period of the flapping wing, and the value is the reciprocal of the flapping frequency; its motion parameters include: flapping frequency f, pitching motion amplitude h, stroke angle β, rotation amplitude α 0 ; the motion parameter space consists of flapping motion parameters that satisfy the above formulas and constraints. 3.根据权利要求1所述的运动参数优化方法，其特征在于，所述步骤(1)中，使用拉丁超立方采样方法在扑翼运动参数构成的参数空间中预设采样点。3. motion parameter optimization method according to claim 1, is characterized in that, in described step (1), use Latin hypercube sampling method to preset sampling point in the parameter space that flapping motion parameter forms. 4.根据权利要求1所述的运动参数优化方法，其特征在于，所述步骤(2)中，数值模拟采用基于径向基函数插值的动网格方法或者浸入边界法。4 . The motion parameter optimization method according to claim 1 , wherein, in the step (2), the numerical simulation adopts a dynamic grid method based on radial basis function interpolation or an immersion boundary method. 5 . 5.根据权利要求1所述的运动参数优化方法，其特征在于，优化目标包括时均升力系数与时均推力系数；升力系数计算公式
Figure FDA0003300758960000041
推力系数计算公式
Figure FDA0003300758960000042
L与T为升力与推力，ρ为流体密度，S为特征面积，U∞为自由来流速度。5. The motion parameter optimization method according to claim 1, wherein the optimization target comprises a time-averaged lift coefficient and a time-averaged thrust coefficient; a lift coefficient calculation formula
Figure FDA0003300758960000041
Thrust coefficient calculation formula
Figure FDA0003300758960000042
L and T are the lift and thrust, ρ is the fluid density, S is the characteristic area, and U ∞ is the free flow velocity.
CN202010005675.0A 2020-01-03 2020-01-03 An optimization method for flapping motion parameters Active CN111199105B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010005675.0A CN111199105B (en) 2020-01-03 2020-01-03 An optimization method for flapping motion parameters

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010005675.0A CN111199105B (en) 2020-01-03 2020-01-03 An optimization method for flapping motion parameters

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 An optimization method for flapping motion parameters

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 Design Method for Suspending Lower Wings of a Hummingbird-like Micro Air Vehicle
CN109783978A (en) * 2019-02-11 2019-05-21 南京航空航天大学 An aerodynamic numerical simulation method of micro flapping-wing aircraft based on ANSYS Workbench
CN110008639A (en) * 2019-04-24 2019-07-12 东莞理工学院 An intelligent parametric design method for the wings of a micro flapping aircraft
CN110334401A (en) * 2019-05-31 2019-10-15 南京航空航天大学 An optimization method for propulsion efficiency of double flapping wings based on tandem arrangement

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 Design Method for Suspending Lower Wings of a Hummingbird-like Micro Air Vehicle
CN109783978A (en) * 2019-02-11 2019-05-21 南京航空航天大学 An aerodynamic numerical simulation method of micro flapping-wing aircraft based on ANSYS Workbench
CN110008639A (en) * 2019-04-24 2019-07-12 东莞理工学院 An intelligent parametric design method for the wings of a micro flapping aircraft
CN110334401A (en) * 2019-05-31 2019-10-15 南京航空航天大学 An optimization method for propulsion efficiency of double flapping wings based on tandem arrangement

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) An optimization method for flapping motion parameters
Jeong et al. Data mining for aerodynamic design space
CN104881510B (en) A kind of lifting airscrew/tail-rotor aerodynamic interference numerical value emulation method
CN105426955A (en) Disturbance-based elite reverse learning particle swarm optimization implementation method
CN110610050A (en) Airfoil aerodynamic drag reduction method based on improved radial basis function deformation algorithm
CN110516318B (en) Airfoil design method based on radial basis function neural network proxy model
CN111581784B (en) A data-driven adaptive quasi-steady-state model for the optimization of flapping wing motion parameters
CN113238666A (en) Ship motion attitude prediction method based on sparrow search algorithm optimization GRU
Peng et al. Multi-fidelity nonlinear unsteady aerodynamic modeling and uncertainty estimation based on Hierarchical Kriging
CN107633105B (en) Improved hybrid frog-leaping algorithm-based quad-rotor unmanned aerial vehicle parameter identification method
CN115809594A (en) Power optimization method and system for floating wind farm based on proxy model
CN113627075A (en) Projectile aerodynamic coefficient identification method based on adaptive particle swarm optimization extreme learning
CN107563107B (en) Static aeroelasticity design method of aircraft structure based on sequential optimization thought
CN110765706B (en) Aerofoil unsteady stall aerodynamic coefficient modeling method based on OHNGBM (1, 1)
Wei et al. Diffairfoil: An efficient novel airfoil sampler based on latent space diffusion model for aerodynamic shape optimization
Sunny et al. An artificial neural network residual kriging based surrogate model for shape and size optimization of a stiffened panel
Trad et al. Airfoils generation using neural networks, CST curves and aerodynamic coefficients
Xiang et al. A manifold-based airfoil geometric-feature extraction and discrepant data fusion learning method
Mukesh et al. Influence of search algorithms on aerodynamic design optimisation of aircraft wings
Thinakaran et al. Predicting the 2-dimensional airfoil by using machine learning methods
CN117556727A (en) Flapping time-varying aerodynamic force prediction method and device based on deep learning
CN114397818A (en) MAPIO-based spacecraft cluster orbit reconstruction path planning method
Jia et al. Optimization of flexible flapping wings for thrust and efficiency
Wei et al. Automatic parameterization for aerodynamic shape optimization via deep geometric learning
CN114186477A (en) Elman neural network-based orbit prediction algorithm

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