Dynamic aerodynamic load distribution method and device based on generalized chebyshev polynomial and electronic equipment

Abstract

The application discloses a dynamic aerodynamic load distribution method and device based on a generalized Chebyshev polynomial and electronic equipment, and the method comprises the following steps: constructing a fluid grid model and a finite element grid model of the same airfoil according to a computational fluid dynamics (CFD) outer shape grid and a finite element model surface grid of the airfoil structure; the two models of the airfoil are segmented respectively, and then a three-dimensional airfoil grid model is mapped to a two-dimensional plane by using an LSCM parameterization method to reduce the subsequent calculation amount; a time variable is introduced to the fluid model node load to obtain a dynamic aerodynamic load, and then a three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on the two-dimensional mapping plane of the fluid grid, so that a fitting coefficient matrix is obtained; and finally, the distribution load on the finite element grid is obtained by interpolation calculation according to the parameterized finite element grid node data. The application has mature theoretical basis and convenient calculation, and can effectively improve the work efficiency and distribution accuracy.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of aviation structure dynamics, and relates to a dynamic aerodynamic load distribution method and device based on a generalized Chebyshev polynomial and an electronic equipment. BACKGROUND

[0002] In the process of aircraft structure strength design, aerodynamic load is a key design input related to the whole. However, how to equivalently and efficiently distribute aerodynamic load to finite element nodes is still a pending issue, so it is mostly relied on manual operation of workers, which is a very important but very tedious work. At present, the most commonly used finite element calculation software in the aerospace industry, such as MSC.Patran / Nastran, Ansys, Abaqus, etc., cannot realize automatic application for the wing surface load with poor uniformity and large changes.

[0003] At present, the most commonly used method for this problem is the "multi-point row" method. This method distributes each aerodynamic load to several finite element nodes with the constraints of minimum strain energy and static force equivalence. However, when facing large wing surfaces and complex wing surfaces, the calculation process of this method will be very tedious. Since this kind of load is used for distributing the concentrated force of aerodynamic load directly integrated on the fluid grid, the modern large aircraft wing structure is complex, and the aerodynamic partition and structural partition are often inconsistent. Direct integration of concentrated load on the fluid grid may cross the area to transfer the load, resulting in an unrealistic force transmission route. Secondly, this method may artificially specify the distribution of some nodes in the face of some special cases, increasing the calculation error. Constructing the distribution function of aerodynamic load or using known functions to fit the distribution of aerodynamic load is also a common way, for example, Zhang Jianguang et al. use spline surface to fit the distribution of aerodynamic load, and get the pressure value on the finite element node from the surface, but there are still problems such as low calculation efficiency and difficult to control fitting accuracy. Finally, the current research is almost for static aerodynamic load, how to introduce the time variable and handle the dynamic aerodynamic load changing with time is also a research direction. SUMMARY

[0004] The purpose of the present application is to overcome the deficiencies of the prior art, provide a dynamic aerodynamic load distribution method, device and electronic equipment based on generalized Chebyshev polynomial. The method introduces a time variable, uses the generalized Chebyshev polynomial as the base function to fit the dynamic aerodynamic load distribution, uses the mature theory for calculation, ensures the calculation efficiency and accuracy; at the same time, in order to improve the calculation efficiency, the fluid grid wing surface and the finite element grid wing surface are first mapped into a two-dimensional plane by using the LSCM (Least Squares Conformal Maps) parameterization method, then the continuity distribution of the dynamic aerodynamic load is fitted on the two-dimensional mapping plane of the fluid grid according to the least square principle, and the load distribution is realized by interpolation calculation according to the parameterized finite element grid data.

[0005] To achieve the above purpose, the present application adopts the following technical solutions.

[0006] The dynamic aerodynamic load distribution method, device and electronic equipment based on generalized Chebyshev polynomial of the present application comprise the following steps:

[0007] Step 1, according to the selected computational fluid dynamics shape grid and finite element model surface grid of the wing surface structure, two virtual structure models of the same wing surface are constructed, i.e. the fluid grid wing surface model and the finite element grid wing surface model; at the same time, the coordinates of the two are unified, so that the Cartesian coordinates of the coincident nodes on the two models are the same. The process includes: first, according to the computational fluid dynamics shape of the wing surface structure, the fluid grid wing surface model of the wing surface is constructed by using modeling software, then the grid is divided on the model by using the Fluent software and the conditions are set for CFD calculation, and the original aerodynamic load is obtained; then, according to the finite element model surface grid of the wing surface structure, the finite element grid wing surface model of the wing surface is constructed by using the MSC.Patran software and is taken as the target carrier of load distribution.

[0008] Step 2, the fluid grid and the finite element grid of the wing surface are respectively mapped from three-dimensional space to two-dimensional space by using the LSCM parameterization method, so as to directly fit the load on the two-dimensional plane subsequently.

[0009] The LSCM (Least Squares Conformal Map) parameterization is the process of flattening the grid surface, which is through the least square approximation of Cauchy-Riemann equation, defines an energy function which makes the angle deformation minimum, correspondingly obtains a linear equation group, and the mapping relationship between the target surface and the plane grid can be constructed by solving the equation group. Because the LSCM has the angle-preserving property, the topological information of the original grid can be well preserved in the plane parameter domain.

[0010] The process of mapping the fluid grid model and the finite element grid model of the airfoil from three-dimensional space to two-dimensional space comprises:

[0011] Step 2.1. The airfoil fluid grid model and the airfoil finite element grid model are respectively segmented; there are two segmentation methods, one is to take three nodes on the airfoil fluid grid model, input the node number, and form a plane according to the coordinates of the three points, and the airfoil fluid grid model and the airfoil finite element grid model are divided into two parts by the plane; the other method is to rely on X-Y, X-Z and Y-Z planes, if X-Y plane is selected, the Z-axis coordinate of the plane is input to determine the position of the plane, so as to divide the airfoil grid into two parts.

[0012] Step 2.2. The four-edge element of the airfoil fluid grid model and the airfoil finite element grid model is uniformly converted into a triangular element, and the four node numbers of each four-edge element are divided into two groups, each group has three node numbers, which converts each four-edge element into two new triangular elements.

[0013] Step 2.3. The airfoil fluid grid model and the airfoil finite element grid model are mapped into a two-dimensional plane by using the LSCM parameterization method, and the grid distribution of the airfoil model fluid grid and the airfoil model finite element grid on the mapping plane is drawn respectively, so as to observe the grid quality after parameterization, and confirm whether the grid shape is consistent. The specific mapping process is as follows:

[0014] Suppose that each triangle has a local coordinate system (x, y), and the normal of the triangle is along the Z-axis direction in the coordinate system, and the coordinates of the three vertices are (x1, y1), (x2, y2), and (x2, y2). Assuming that U:(x, y)→(u, v) is a mapping from three-dimensional space to two-dimensional plane, then in the local coordinate system, the Cauchy-Riemann equation is discretized as:

[0015]

[0016] X is represented by complex number as:

[0017] X=x+iy (2)

[0018] For the mapped triangle, it can also be represented by a local coordinate system, and by analogy:

[0019] U=u+iv (3)

[0020] According to the inverse function derivative theorem, it is obtained that:

[0021]

[0022] The gradient in the triangle can be described by the following formula:

[0023]

[0024] where A is the area of the triangle. A can be calculated as follows:

[0025] A = (x1y2 - y1x2) + (x2y3 - y2x3) + (x2y3 - y2x3) (6)

[0026] Combining equation (4) with the gradient formula in the triangle, we can get:

[0027]

[0028] where the values of W1, W2, W3 have nothing to do with the selection of the local coordinate system of the triangle. They are vectors between the coordinates, and they are perpendicular to U j The expression of U

[0029]

[0030] Since most of the triangular meshes are not completely consistent before and after mapping in the parameterization process, i.e., they cannot completely satisfy equation (4). Therefore, this method also needs a conformal energy C to control the deformation degree of the triangle before and after mapping. When the conformal energy C is the smallest, the conformal mapping is satisfied as much as possible. By solving the unique minimum value of the conformal energy function, the two-dimensional parameterization of the three-dimensional surface is completed. The process of constructing the conformal energy is described as follows:

[0031] The conformal energy of each triangle is:

[0032]

[0033] where A T is the area of the triangle, n is the total number of triangles, and k can be regarded as the number of triangular elements.

[0034] For the convenience of calculating the minimum value of the conformal energy, it can be converted into a quadratic form, and the total conformal energy of the grid is written as:

[0035]

[0036] The matrix M in the above equation is (M kq ), where:

[0037]

[0038] Each triangular mesh node is numbered q = (1, 2, 3...l), and l is the total number of nodes. M and U in equation (10) are complex matrices. For triangle k, the three node numbers are q1, q2, q3, and the corresponding values of the M matrix in the kth row and q1, q2, q3 columns are The rest of the row is 0. Meanwhile, M and U can be written as complex numbers as follows:

[0039]

[0040] where A is the real part matrix of M, B is the imaginary part matrix of M, u is the real part vector of U, and v is the imaginary part vector of U. Then, substituting equation (12) into equation (10), the conformal energy can be written as:

[0041]

[0042] When the quadratic form takes the minimum value, there is a stationary point:

[0043]

[0044] At this point, the matrix L C The matrix U containing the mapped node coordinate information can be obtained according to equation (14) from the initial node and element node information.

[0045] Step 3, the original aerodynamic load is processed, a discrete time variable is introduced, and a dynamic aerodynamic load that changes nonlinearly with a given time node is obtained; a three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on the two-dimensional mapping plane of the fluid grid according to the parameterized fluid grid data, the coordinate fitting order and the time fitting order are set, and the least square optimization algorithm is used to calculate the polynomial fitting coefficient matrix; the polynomial fitting coefficient matrix and the parameterized finite element grid data are used for interpolation calculation, so as to obtain the dynamic distribution load of the finite element grid; the two parts of the finite element grid model load are combined and corresponded with the original finite element grid node, the resultant force and the pressure center position of the fluid grid model and the finite element grid model are calculated respectively, and the error is calculated. The maximum error of the resultant force and the pressure center position is set to 5%, if the calculation result error does not meet the requirement, the fitting order is adjusted and the above fitting steps are repeated.

[0046] The process specifically includes:

[0047] Step 3.1. Take the node load of the fluid grid model calculated by the Fluent software as P0, multiply it by the sine function f(t) = sint with time t as the independent variable to obtain the load that changes with time, denoted as P(t); take its time nodes in the interval (0, π / 2] at a certain interval to form a time node array Q = [t1, t2, t3,..., t s ]s is the number of time nodes, and take the load values corresponding to each time node to form a dynamic load array, denoted as

[0048] S = [P(t1), P(t2), P(t3),..., P(ts )] (15)

[0049] Step 3.2. Combine the parameterized fluid mesh model node new coordinates with its original corresponding node's nodal dynamic load, which is described as a dynamic load function q(x, y, t) with three variables; use three-dimensional first-type Chebyshev polynomials to fit the function, which based on the expression of the first m x n order Chebyshev terms (m is the coordinate fitting order, n is the time fitting order) can be written as:

[0050]

[0051] Where the i x j order expression is as follows:

[0052]

[0053] The general term formula of i order Chebyshev polynomial in the above formula is as follows:

[0054]

[0055] In the load function expression, T ikl (x, y, t) = T k (x)T i-k (y)T l (t) (k = 0, 1, 2... i, l = 0, 1,..., j), a ikl is the coefficient of each term in q(x, y, t), let:

[0056] A i = [a i00 a i01 ... a i0n a i10 a i11 ... a i1n ...... a ii0 a ii1 ... a iin ] T (19)

[0057] A = [A0 A1 A2... A m ] T (20)

[0058] T i = [T i00 T i01 ... T i0n T i10 T i11 ... T i1n ...... Tii0 T ii1 ... T iin ] (21)

[0059] T=[T0 T1 T2... T m ] (22)

[0060] then equation (16) can be written as:

[0061] q(x,y,t) = TA (23)

[0062] Step 3.3. Obtain the coefficient matrix of polynomial fitting by using least square method.

[0063] According to the theory of least square method, the function q(x,y,t) is not always through the known data points, but is obtained according to the difference between the value of fitting function at sampling and the actual value, i.e.:

[0064]

[0065] In the formula, r is the node number, p rs is the aerodynamic load value of node numbered r at time t s , N is the total number of nodes, M is the total number of time nodes, and E(q) is the calculation error square sum. In order to make E(q) of a time node t s take the minimum value, then:

[0066]

[0067] Expand equation (25) into linear transformation of linear equations and write it in matrix form, then:

[0068] BB T A = BP (26)

[0069] Where:

[0070]

[0071] P T = [p 1s p 2s ... p Ns ] (28)

[0072] Let E = (m+1)(m+2) / 2, F = n x N, B is an E x F matrix, A is an E-dimensional column vector, and P is an F-dimensional column vector. According to equation (26), each element in A can be obtained:

[0073] A = (BB T ) -1 BP (29)

[0074] Step 3.4. Calculate the load value of the finite element grid model node according to the obtained coefficient matrix and the finite element grid data. After obtaining the coefficient matrix A of the dynamic load function q(x, y, t), the load value of the finite element grid node is obtained by interpolation calculation using the finite element grid data, i.e.

[0075] R = (B') T A (30)

[0076] Wherein, R is the load value of each node on the finite element grid, B' is the B matrix calculated according to the node coordinates of the finite element grid on the mapping plane. At this point, the calculation of fitting and interpolation is completed.

[0077] Step 3.5. Since the obtained distribution load belongs to different partitions of the finite element grid model, it is first combined with the original nodes of the finite element grid model one by one to form the node load of the complete finite element grid model. According to the node load array and node and element information of the fluid grid model, the resultant force F c and the pressure center (X c , Y c ) of the aerodynamic load are calculated, and then the resultant force F' c and the pressure center (X' c , Y' c ) of the distribution load are calculated according to the node distribution load array and node and element information of the finite element grid model. The errors are calculated respectively:

[0078]

[0079]

[0080]

[0081] Wherein, W F is the error of the distribution load resultant force and the aerodynamic load resultant force, W X and W Y are the errors of the distribution load pressure center and the aerodynamic load pressure center in the chordwise and spanwise positions respectively, X c and Y c are the positions of the aerodynamic load pressure center in the chordwise and spanwise positions respectively, X' c and Y' c are the positions of the calculated distribution load pressure center in the chordwise and spanwise positions respectively, L a and L b are the maximum lengths of the chordwise and spanwise positions of the wing surface model respectively.

[0082] If the final calculated error fails to meet the predetermined requirement of 5%, the polynomial fitting order is adjusted to repeat steps 3.2 to 3.5 until the error finally meets the requirement.

[0083] The dynamic aerodynamic load distribution device based on the generalized Chebyshev polynomial comprises:

[0084] The wing surface model construction module is used to construct two virtual structure models of the same wing surface, i.e., a fluid grid wing surface model and a finite element grid wing surface model, according to the computational fluid dynamics shape grid and the finite element model surface grid of the selected wing surface structure; the Cartesian coordinates of the coincident nodes on the two models are made the same in the unified coordinates; wherein, the fluid grid wing surface model of the wing surface is constructed according to the computational fluid dynamics shape of the wing surface structure using modeling software, and then the Fluent software is used to divide the grid on the model and set the conditions for CFD calculation to obtain the original aerodynamic load; then, the finite element grid wing surface model of the wing surface is constructed using the MSC.Patran software according to the finite element model surface grid of the wing surface structure and is taken as the target carrier for load distribution;

[0085] The grid model mapping module is used to map the fluid grid model and the finite element grid model of the wing surface from three-dimensional space to two-dimensional space respectively by using the LSCM parameterization method, so as to perform load fitting directly on the two-dimensional plane subsequently;

[0086] The load distribution module is used to process the original aerodynamic load, introduce a discrete time variable, and obtain the dynamic aerodynamic load which changes nonlinearly with the given time nodes; the three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on the two-dimensional mapping plane of the fluid grid according to the parameterized fluid grid data, the coordinate fitting order and the time fitting order are set, and the least square optimization algorithm is used to calculate the polynomial fitting coefficient matrix; the polynomial fitting coefficient matrix and the parameterized finite element grid data are used for interpolation calculation, so as to obtain the dynamic distribution load of the finite element grid; the two parts of the finite element grid model load are combined and corresponded to the original finite element grid nodes, the resultant force and the pressure center position of the fluid grid model and the finite element grid model are calculated respectively, and the error is calculated; if the calculation result error does not meet the requirement, the fitting order is adjusted to repeat the above fitting steps.

[0087] The electronic device comprises a memory, a processor and a computer program stored on the memory and capable of running on the processor, and the processor implements the dynamic aerodynamic load distribution method based on the generalized Chebyshev polynomial when the computer program is executed.

[0088] A non-transitory computer-readable storage medium of the present invention stores a computer program thereon, wherein when the computer program is executed by a processor, the dynamic aerodynamic load distribution method based on generalized Chebyshev polynomials as described above is implemented.

[0089] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0090] 1. The present invention uses the LSCM parameterization method to first map the three-dimensional airfoil onto a two-dimensional plane, which greatly reduces the computational complexity of subsequent dynamic load fitting. Due to its inherent characteristics, it also offers the following advantages and benefits compared to other mapping methods:

[0091] (1) Reduces angular distortion and inconsistent scaling during the mapping process;

[0092] (2) There exists a unique conformal energy minimum, avoiding local optimal solutions;

[0093] (3) There is no need to determine the boundary mapping, so it can act on boundaries of arbitrary shapes.

[0094] 2. This invention uses generalized Chebyshev orthogonal polynomials as basis functions to fit dynamic aerodynamic load distribution. This method is simple and mature, achieving precise load distribution. Furthermore, due to the orthogonality of these polynomials, a large number of complex matrix operations can be eliminated, improving computational efficiency. Furthermore, the polynomial fitting order can be adjusted multiple times to achieve better results, fully ensuring the accuracy of the results.

[0095] 3. The present invention introduces a time variable into the original static aerodynamic load, attempting to fit and distribute the dynamic aerodynamic load that changes over time. This invention can be further developed to verify the feasibility of distributing dynamic aerodynamic loads and improve the dynamic aerodynamic load distribution technology to meet the demand for refined load distribution. BRIEF DESCRIPTION OF THE DRAWINGS

[0096] Figure 1 It is a method flow chart of an embodiment of the present invention.

[0097] Figure 2 1 is a schematic diagram of a fluid grid of an airfoil model according to an embodiment of the present invention, which includes 10,650 nodes and 10,500 grids.

[0098] Figure 3 3 is a schematic diagram of a finite element mesh of an airfoil model according to an embodiment of the present invention, which includes 3075 nodes and 3000 meshes.

[0099] Figure 4 This is a diagram of the original aerodynamic load distribution on the upper wing surface of an embodiment of the present invention.

[0100] Figure 5 is a lower surface original aerodynamic load distribution diagram of an embodiment of the present application.

[0101] Figure 6 is a load conversion error display diagram of different coordinate fitting orders of an embodiment of the present application, the force and the center of pressure position error decrease with the increase of the fitting order, the force error starts to stabilize after the 7th order, and the center of pressure position error always keeps a low level. DETAILED DESCRIPTION

[0102] In the load distribution between different grid models of the airfoil related to the structural dynamics technology, there are problems of complex calculation and low accuracy of the existing method and the distribution processing of the dynamic aerodynamic load changing with time. The present application firstly uses the LSCM parameterization method to map the three-dimensional airfoil grid to the two-dimensional plane to reduce the subsequent calculation amount, and introduces the time variable to the aerodynamic load of the existing fluid grid airfoil model to use the three-dimensional generalized Chebyshev polynomial to fit the dynamic node load function of the airfoil fluid grid model, and after obtaining the fitting coefficient matrix according to the least square principle, the distribution load on the finite element grid is calculated through interpolation according to the finite element grid model information.

[0103] The present application will be further described in detail below in combination with the drawings.

[0104] As shown in Figure 1 , the present application starts from the airfoil fluid grid model and the airfoil finite element grid model to be distributed, and is divided into two stages. The first stage is to segment and do LSCM parameterization processing on the two airfoil models, corresponding to step 2; the second stage is to use three-dimensional Chebyshev polynomial to fit the dynamic aerodynamic load on the plane after parameterization of the airfoil fluid grid, and then combine the airfoil finite element grid parameterization data to obtain the distribution load of each node on the airfoil finite element grid model, corresponding to step 3.

[0105] The method specifically comprises the following steps:

[0106] Step 1, according to the selected wing surface structure of computational fluid dynamics shape grid and finite element model surface grid, two virtual structure models of the same wing surface are constructed, namely the fluid grid wing surface model and the finite element grid wing surface model; at the same time, the coordinates of the two are unified, so that the Cartesian coordinates of the coincident nodes on the two models are the same. The process includes: first, according to the computational fluid dynamics shape of the wing surface structure, the fluid grid wing surface model of the wing surface is constructed using modeling software, then the Fluent software is used to divide the grid on the model and set the conditions for CFD calculation to obtain the original aerodynamic load; then, according to the finite element model surface grid of the wing surface structure, the finite element grid wing surface model of the wing surface is constructed using MSC.Patran software and used as the target carrier of load distribution.

[0107] Step 2, the fluid grid model and the finite element grid model of the wing surface are respectively mapped from three-dimensional space to two-dimensional space using the LSCM parameterization method, so that subsequent load fitting can be directly performed on the two-dimensional plane.

[0108] The process of mapping the fluid grid model and the finite element grid model of the wing surface from three-dimensional space to two-dimensional space includes:

[0109] Step 2.1. The wing surface fluid grid model and the wing surface finite element grid model are respectively segmented; the first segmentation method is selected in the method, three nodes on the wing surface fluid grid model are taken, a plane is formed according to the coordinates of the three points, and the wing surface fluid grid model and the wing surface finite element grid model are divided into two approximately symmetrical parts, called upper wing surface and lower wing surface; this step can simplify the subsequent calculation and speed up the calculation.

[0110] Step 2.2. The four-sided element of the wing surface fluid grid model and the wing surface finite element grid model is uniformly converted into a triangular element, and the four node numbers of each four-sided element are divided into two groups, each group having three node numbers, which converts each four-sided element into two new triangular elements.

[0111] Step 2.3. The wing surface fluid grid model and the wing surface finite element grid model are mapped to a two-dimensional plane using the LSCM parameterization method, and the grid distribution of the wing surface model fluid grid and the wing surface model finite element grid on the mapping plane is drawn respectively to observe the grid quality after parameterization and confirm whether the grid shape is consistent. The specific mapping process is as follows:

[0112] Let each triangle have a local coordinate system (x, y), in which the normal of the triangle is along the Z-axis direction, and the coordinates of the three vertices are (x1, y1), (x2, y2), and (x2, y2). Assuming that U:(x, y)→(u, v) is a mapping from three-dimensional space to two-dimensional plane, then in the local coordinate system, the Cauchy-Riemann equation is discretized as:

[0113]

[0114] X is expressed as a complex number:

[0115] X = x + iy (2)

[0116] For the mapped triangle, it can also be expressed in a local coordinate system, and similarly:

[0117] U = u + iv (3)

[0118] According to the inverse function derivative theorem, we have:

[0119]

[0120] The gradient inside the triangle can be described by the following formula:

[0121]

[0122] where A is the area of the triangle. It can be calculated according to the following formula:

[0123] A = (x1y2 - y1x2) + (x2y3 - y2x3) + (x2y3 - y2x3) (6)

[0124] Combining formula (4) with the gradient formula inside the triangle, we can derive:

[0125]

[0126] where W1, W2, W3 have nothing to do with the selection of the local coordinate system of the triangle. They are vectors between coordinates, and their relationship with U j is as follows:

[0127]

[0128] Since most of the triangular meshes are not completely consistent before and after mapping in the parameterization process, i.e., they do not completely satisfy formula (4). Therefore, this method also needs a conformal energy C to control the deformation degree of the triangle before and after mapping. When the conformal energy C is the smallest, it satisfies the conformal mapping as much as possible. By solving the unique minimum value of this conformal energy function, the two-dimensional parameterization of the three-dimensional surface is completed. The process of constructing the conformal energy is described as follows:

[0129] The conformal energy of each triangle is:

[0130]

[0131] where A T is the area of the triangle, n is the total number of triangles, and k can be regarded as the number of triangular elements.

[0132] For convenience of calculating the minimum value of conforming energy, it can be converted into a quadratic form, and the total conforming energy of the grid is written as:

[0133]

[0134] The matrix M in the above formula is (M kq ), wherein:

[0135]

[0136] Each triangular grid node is numbered q=(1, 2, 3...l), and l is the total number of nodes. M and U in formula (10) are complex matrices. For triangular grid k, the three node numbers are q1, q2, q3, and the corresponding values of the M matrix in the kth row and q1, q2, q3 columns are The rest of the positions in this row are 0. At the same time, M and U can be written in the following complex form:

[0137]

[0138] Wherein A is the real matrix of matrix M, B is the imaginary matrix of matrix M, u is the real part vector of U, and v is the imaginary part vector of U. Then formula (12) is substituted into formula (10), and the conforming energy can be written as:

[0139]

[0140] When the quadratic form takes the minimum value, there is a stationary point:

[0141]

[0142] At this point, the matrix L C Can be obtained from the initial node and element node information, and according to formula (14), the column matrix U containing the mapped node coordinate information can be obtained.

[0143] Step 3, the original aerodynamic load is processed, and a discrete time variable is introduced to obtain a dynamic aerodynamic load that varies nonlinearly with a given time node; a three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on the two-dimensional mapping plane of the fluid grid according to the parameterized fluid grid data, the coordinate fitting order and the time fitting order are set, and the least square optimization algorithm is used to calculate the polynomial fitting coefficient matrix; the polynomial fitting coefficient matrix and the parameterized finite element grid data are used for interpolation calculation, so as to obtain the dynamic distribution load of the finite element grid; the two parts of the finite element grid model load are combined and corresponded with the original finite element grid node, the resultant force and the pressure center position of the fluid grid model and the finite element grid model are calculated respectively, and the error is calculated. The maximum error of the resultant force and the pressure center position is set to 5%, if the calculation result error does not meet the requirement, the fitting order is adjusted and the above fitting steps are repeated. The process specifically includes:

[0144] Step 3.1. Take the node load of the fluid grid model calculated by the Fluent software as P0, multiply it by the sine function f(t)=sint with time t as the independent variable to obtain the load varying with time, denoted as P(t); take the time nodes in the interval (0, π / 2] at a certain interval to form a time node array Q=[t1,t2,t3,...,t s ], s is the number of time nodes, and take the load values corresponding to each time node to form a dynamic load array, denoted as

[0145] S=[P(t1),P(t2),P(t3),...,P(t s )] (15)

[0146] Step 3.2. Combine the parameterized fluid grid model node new coordinates with the node dynamic load of the original corresponding node, and describe it as a dynamic load function q(x,y,t) containing three variables; a three-dimensional first type Chebyshev polynomial is used to fit the function, which is based on the expression of the first m×n order Chebyshev polynomial (m is the coordinate fitting order, and n is the time fitting order), and can be written as:

[0147]

[0148] Wherein, the i×j order expression is as follows:

[0149]

[0150] The general term formula of the i order Chebyshev polynomial in the above formula is as follows:

[0151]

[0152] In the expression of the load function, T ikl (x,y,t) = T k (x)T i-k (y)T l (t) (k = 0, 1, 2...i, l = 0, 1,...,j), a ikl is the coefficient of each term in q(x,y,t), let:

[0153] A i = [a i00 a i01 ... a i0n a i10 a i11 ... a i1n ...... a ii0 a ii1 ... a iin ] T (19)

[0154] A = [A0 A1 A2... A m ] T (20)

[0155] T i = [T i00 T i01 ... T i0n T i10 T i11 ... T i1n ...... T ii0 T ii1 ... T iin ] (21)

[0156] T = [T0 T1 T2... T m ] (22)

[0157] Then, equation (16) can be written as:

[0158] q(x,y,t) = TA (23)

[0159] Step 3.3. Obtain the coefficient matrix of the polynomial fitting by using the least square method.

[0160] According to the theory of the least square method, the function q(x,y,t) is not generally obtained by all the known data points, but is obtained by making the square sum of the difference between the value of the fitting function at the sampling point and the actual value minimum, i.e.:

[0161]

[0162] In the equation, r is the node number, prs The aerodynamic load value of the node numbered r at time t s N is the total number of nodes, M is the total number of time nodes, and E(q) is the calculation error square sum. To make E(q) of a time node t s take the minimum value, then:

[0163]

[0164] The linear transformation of formula (25) into a linear equation group and writing it in a matrix form, then:

[0165] BB T A = BP (26)

[0166] Wherein:

[0167]

[0168] P T = [p 1s p 2s ... p Ns ] (28)

[0169] Let E = (m+1)(m+2) / 2, F = n x N, B is an E x F matrix, A is an E-dimensional column vector, and P is an F-dimensional column vector. According to formula (26), each element in A can be solved:

[0170] A = (BB T ) -1 BP (29)

[0171] Step 3.4. According to the obtained coefficient matrix and finite element grid data, the load value of the finite element grid model node is obtained. After obtaining the coefficient matrix A of the dynamic load function q(x, y, t), the load value of the finite element grid node is obtained by interpolation calculation using the finite element grid data, that is:

[0172] R = (B′) T A (30)

[0173] Wherein, R is the load value of each node on the finite element grid, and B' is the B matrix calculated according to the node coordinates of the finite element grid on the mapping plane. Thus, the calculation of fitting and interpolation is completed.

[0174] Step 3.5. Since the obtained distribution load belongs to different partitions of the finite element grid model, it is first combined with the original nodes of the finite element grid model one by one to form the node load of the complete finite element grid model. According to the node load array and node and element information of the fluid grid model, the resultant force F c and the pressure center (Xc ,Y c ), and then calculate the resultant force F' of the distributed load according to the node distribution load array and node and unit information of the finite element mesh model c and pressure center (X' c ,Y' c ). Find their respective errors:

[0175]

[0176]

[0177]

[0178] Among them, W F is the error between the resultant force of the distributed load and the resultant force of the aerodynamic load, W X 、W Y are the errors of the distribution load pressure center and the aerodynamic load pressure center in the chord direction and span direction of the wing, respectively, X c 、Y c are the positions of the pressure center of the aerodynamic load in the chord direction and span direction of the airfoil, X' c 、Y' c are the calculated positions of the distributed load pressure center in the chord direction and span direction of the airfoil, L a 、L b are the maximum lengths of the airfoil model in the chord direction and span direction, respectively.

[0179] If the final calculated error fails to meet the predetermined requirement of 5%, the polynomial fitting order is adjusted and steps 3.2 to 3.5 are repeated until the error finally meets the requirement.

[0180] A verification example based on the method of the present invention:

[0181] The ONERA M6 standard trapezoidal wing (M6 wing for short) was selected, and its specific parameters are: the model half-span (end) length is 1196.3mm, the root chord length is 810.49mm, the blade tip chord length (end) is 455.91mm, the leading edge curvature is 330°, the trailing edge curvature is 15.69°, and the fillet tip span is 1218.535mm.

[0182] Figure 2 This is a schematic diagram of the fluid mesh model of the wing surface of this verification example, which contains 10,650 aerodynamic load nodes and 10,500 meshes. Figure 3 Schematic diagram of the structural finite element mesh model of the wing, which contains 3075 finite element nodes and 3000 units.

[0183] Figure 4 、 Figure 5The original aerodynamic load distribution diagram of the present verification example is shown in Figure 1. The load along the chord direction generally decreases first and then increases, and the phenomenon is more obvious in the upper surface part.

[0184] Figure 6 The load conversion error under different coordinate fitting orders is shown. Since the influence of time order on fitting error in the present method is far less than that of coordinate order, only the influence of different coordinate orders on load conversion error is shown. It can be seen that the force and pressure center position error decrease with the increase of fitting order. The force error starts to stabilize after the 7th order, while the pressure center position error remains at a low level.

[0185] The aerodynamic load is distributed by the load distribution method based on generalized Chebyshev polynomial at each time node. When the fitting order is 7x10, the fitting error has started to stabilize, and when the fitting order is 14x10, the fitting accuracy is already good, and the specific distribution of the load is also roughly the same. The specific results are shown in Tables 1-3 below.

[0186] Table 1. Force and error of load distribution at different time points

[0187]

[0188] Table 2. Chordwise pressure center position and error of load distribution at different time points

[0189]

[0190] Table 3. Spanwise pressure center and error of load distribution at different time points

[0191]

[0192] From the above results, it can be seen that the present method maintains high accuracy when fitting and distributing aerodynamic loads at different time nodes, and thus has the ability to distribute dynamic aerodynamic loads that change with time. At the same time, the segmentation and parameterization of the wing surface model can effectively solve the problem of time-consuming load distribution calculation caused by complex model structure and large number of nodes.

Claims

1. A dynamic aerodynamic load distribution method based on generalized Chebyshev polynomials, characterized by, The method comprises the following steps: Step 1, according to the selected airfoil structure of the computational fluid dynamics outer shape grid and the finite element model surface grid, two virtual structure models of the same airfoil surface are constructed, namely a fluid grid airfoil model and a finite element grid airfoil model; the coordinates are unified so that the Cartesian coordinates of the coincident nodes on the two models are the same; Step 2, the fluid grid and the finite element grid of the airfoil model are respectively mapped from three-dimensional space to two-dimensional space by using the LSCM parameterization method, so as to directly perform load fitting on the two-dimensional plane subsequently; Step 3, the original aerodynamic load is processed, a discrete time variable is introduced, and a dynamic aerodynamic load changing nonlinearly with a given time node is obtained; a three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on the two-dimensional mapping plane of the fluid grid according to the parameterized fluid grid data, the coordinate fitting order and the time fitting order are set, and the least square method is used to calculate the polynomial fitting coefficient matrix; the polynomial fitting coefficient matrix and the parameterized finite element grid data are used for interpolation calculation, so as to obtain the dynamic distribution load of the finite element grid; the two parts of the finite element grid model are combined and corresponded with the original finite element grid nodes, the resultant force and the pressure center position of the fluid grid model and the finite element grid model are calculated respectively, and the error is calculated; the maximum error of the resultant force and the pressure center position is set to 5%, if the calculation result error does not meet the requirement, the fitting order is adjusted and the above fitting steps are repeated; In step 2, the process of mapping the fluid grid model and the finite element grid model of the airfoil from three-dimensional space to two-dimensional space comprises: Step 2.1, the airfoil fluid grid model and the airfoil finite element grid model are respectively divided; there are two division methods, one is to take three nodes on the airfoil fluid grid model, input the node number, and form a plane according to the coordinates of the three points, so as to divide the airfoil fluid grid model and the airfoil finite element grid model into two parts; the other method is to rely on X-Y, X-Z and Y-Z planes, if the X-Y plane is selected, the Z-axis coordinate of the plane is input to determine the position of the plane, so as to divide the airfoil grid into two parts; Step 2.2, the four-edge element of the airfoil fluid grid model and the airfoil finite element grid model is uniformly converted into a triangular element, and the four node numbers of each four-edge element are divided into two groups, each group has three node numbers, so as to convert each four-edge element into two new triangular elements; Step 2.3, the LSCM parameterization method is used to map the airfoil fluid grid model and the airfoil finite element grid model to a two-dimensional plane, and the grid distribution of the airfoil model fluid grid and the airfoil model finite element grid on the mapping plane is drawn respectively, so as to observe the grid quality after parameterization and confirm whether the grid shapes are consistent; the specific mapping process is as follows: Let each triangle has a local coordinate system (x, y), in which the normal of the triangle is along the Z-axis direction, and the coordinates of the three vertices are (x1, y1), (x2, y2), (x2, y2) respectively; suppose U: (x, y)→(u, v) is the mapping from three-dimensional space to two-dimensional plane, then in the local coordinate system, the Cauchy-Riemann equation is discretized as: X is expressed by complex number as: X = x + iy (2) For the mapped triangle, also use the local coordinate system to represent, and similarly have: U = u + iv (3) According to the inverse function derivative theorem, it is obtained that: The gradient in the triangle can be described by the following formula: In the formula, A is the area of the triangle; A can be calculated according to the following formula: A = (x1y2 - y1x2) + (x2y3 - y2x3) + (x2y3 - y2x3) (6) The formula (4) is combined with the gradient formula in the triangle to derive: Wherein, the value of W1, W2, W3 has nothing to do with the selection of the local coordinate system of the triangle, they are vectors between coordinates, and U j The expression of is as follows: U j = u j + iv j (j = 1, 2, 3) Since most of the triangular meshes are not completely consistent before and after the mapping in the parameterization process, that is, formula (4) cannot be completely satisfied; therefore, a conformal energy C is needed to control the deformation degree of the triangle before and after the mapping, and when the conformal energy C is the minimum, the conformal mapping is as possible as to meet the condition; by solving the unique minimum value of the conformal energy function, the two-dimensional parameterization of the three-dimensional surface is completed; the process of constructing the conformal energy is described as follows: The conformal energy of each triangle is: where A T is the area of the triangle, n is the total number of triangles, and k is considered as the number of the triangle element; In order to facilitate the calculation of the minimum value of the conformal energy, it is converted into a quadratic form, and the total conformal energy of the grid is written as: The matrix M = (M kq ) in the above equation, where: Each triangular mesh node is numbered q = (1, 2, 3...l), l is the total number of nodes; M and U in formula (10) are complex matrices, for a triangular mesh k, three node numbers are q1, q2, q3, then the kth row of the M matrix has corresponding values in the q1, q2, q3 columns The rest of the positions in the row are 0; meanwhile, M and U are written in the following complex form: Where A is the real matrix of matrix M, B is the imaginary matrix of matrix M, u is the real part vector of U, and v is the imaginary part vector of U; then formula (12) is substituted into formula (10), and the conformal energy is written as: When the quadratic form takes the minimum value, there is a stationary point: ▽C = L C U = 0 (14) Thus far, the matrix L C From the initial node and cell node information, the matrix U containing the mapped node coordinate information is obtained according to equation (14).

2. The dynamic airload distribution method based on generalized Chebyshev polynomials as claimed in claim 1, wherein, In step 1, the process of constructing the fluid grid wing surface model and the finite element grid wing surface model of the same wing surface includes: First, according to the computational fluid dynamics shape of the wing surface structure, a fluid grid wing surface model of the wing surface is constructed using modeling software, then the Fluent software is used to divide the grid on the model and set the conditions for CFD calculation to obtain the original aerodynamic load; then, according to the surface grid of the finite element model of the wing surface structure, the MSC.Patran software is used to construct the finite element grid wing surface model of the wing surface and take it as the target carrier of load distribution.

3. The dynamic airload distribution method based on generalized Chebyshev polynomials as claimed in claim 1, wherein, In step 3, the process of constructing and distributing the dynamic aerodynamic load includes: Step 3.

1. Take the node load of the fluid grid model calculated by Fluent software as P0, multiply it by the sine function f(t) = sint with time t as the independent variable to obtain the load varying with time, denoted as P(t); take its time nodes in the interval (0, π / 2] at a certain interval to form a time node array Q = [t1, t2, t3,..., t s ]s is the number of nodes, and take the load values corresponding to each time node to form a dynamic load array, denoted as S = [P(t1), P(t2), P(t3),..., P(t s )] (15) Step 3.

2. Combine the node coordinates of the parameterized fluid grid model with the node dynamic aerodynamic load of the original corresponding node, and describe it as a dynamic load function q(x, y, t) containing three variables; use the three-dimensional first type Chebyshev polynomial to fit the function, which is based on the expression of the first m×n order Chebyshev polynomial (m is the coordinate fitting order, and n is the time fitting order) and is written as: Where the i×j order expression is as follows: The general term formula of the i order Chebyshev polynomial in the above formula is as follows: In the expression of the load function, T ikl (x,y,t) = T k (x)T i-k (y)T l (t)(k = 0, 1, 2...i, l = 0, 1,...,j), a ikl is the coefficient of each term in q(x,y,t), let: A i = [a i00 a i01 ...a i0n a i10 a i11 ...a i1n ......a ii0 a ii1 ...a iin ] T (19) A = [A0 A1 A2...A m ] T (20) T i = [T i00 T i01 ...T i0n T i10 T i11 ...T i1n ......T ii0 T ii1 ...T iin ] (21) T = [TO T1 T2...TN] (21) m ] (22) Then formula (16) can be written as: q(x, y, t) = TA (23) Step 3.

3. The coefficient matrix of the polynomial fitting is solved by the least square method; According to the least square method theory, the function q(x, y, t) is not obtained by the known data points, but is obtained by making the square sum of the difference between the value of the fitting function at the sampling point and the actual value minimum, that is: where r is the node number, p rs is the aerodynamic load value of the node numbered r at time t s , N is the total number of nodes, M is the total number of time nodes, E(q) is the calculation error square sum; to make E(q) of a time node t s take the minimum value, then: The linear transformation of formula (25) is expanded into a linear equation group, and is written in a matrix form, and has: BB T A = BP (26) Wherein: P T = [p 1s p 2s ...p Ns ] (28) Let E=(m+1)(m+2) / 2, F=n×N, B is an E×F matrix, A is an E-dimensional column vector, and P is an F-dimensional column vector, and each element in A is obtained according to formula (26): A = (BB T ) -1 BP (29) Step 3.

4. According to the obtained coefficient matrix and the finite element grid node data, the load distribution value of the finite element grid model node is obtained; after obtaining the coefficient matrix A of the dynamic load function q(x, y, t), the load distribution value of the finite element grid node is obtained by using the finite element grid data for interpolation calculation, that is: R = (B') T A (30) Wherein, R is the load value of each node on the finite element grid, and B' is the B matrix calculated according to the node coordinates of the finite element grid on the plane; at this time, the calculation of fitting and interpolation is completed; Step 3.

5. Since the obtained distribution load belongs to different partition of finite element grid model, it is combined with the original node of finite element grid model one by one to form the node load of complete finite element grid model; the resultant force F of aerodynamic load is calculated according to the node load array and node and unit information of fluid grid model c and pressure center (X c , Y c ), and the resultant force F of distribution load is calculated according to the node distribution load array and node and unit information of finite element grid model c ' and pressure center (X c ', Y c ); the respective error is calculated: where W F is the error of the distribution load resultant force and the aerodynamic load resultant force, W X , W Y are the errors of the distribution load pressure center and the aerodynamic load pressure center in the chordwise and spanwise positions of the wing, respectively, X c , Y c are the positions of the pressure center of the aerodynamic load in the chordwise and spanwise of the airfoil, respectively, X c ', Y c ' are the positions of the calculated distribution load pressure center in the chordwise and spanwise of the airfoil, respectively, L a , L b are the maximum lengths of the airfoil model in the chordwise and spanwise, respectively. If the error of the final calculation cannot meet the predetermined requirement of 5%, the order of the polynomial fitting is adjusted, and steps 3.2 to 3.5 are repeated until the error meets the requirement.

4. A dynamic aerodynamic load distribution device based on generalized Chebyshev polynomials, characterized by Comprise: The wing surface model construction module is used for constructing two virtual structure models of the same wing surface, i.e. the fluid grid wing surface model and the finite element grid wing surface model, according to the selected computational fluid dynamics shape grid and the finite element model surface grid of the wing surface structure; the Cartesian coordinates of the coincident nodes on the two models are the same; wherein, the fluid grid wing surface model of the wing surface is constructed according to the computational fluid dynamics shape of the wing surface structure, the modeling software is used to construct the fluid grid wing surface model of the wing surface, the Fluent software is used to divide the grid on the model and set the conditions for CFD calculation, and the original aerodynamic load is obtained; then, the finite element grid wing surface model of the wing surface is constructed according to the finite element model surface grid of the wing surface structure, and the model is used as a target carrier for load distribution; The grid model mapping module is used for mapping the fluid grid model and the finite element grid model of the wing surface from three-dimensional space to two-dimensional space by using the LSCM parameterization method, so as to directly perform load fitting on the two-dimensional plane subsequently; The grid model mapping module is used for mapping the fluid grid model and the finite element grid model of the wing surface from three-dimensional space to two-dimensional space by using the LSCM parameterization method, so as to directly perform load fitting on the two-dimensional plane subsequently; The load distribution module is used for processing the original aerodynamic load, introducing a discrete time variable, and obtaining a dynamic aerodynamic load which varies nonlinearly with a given time node; a three-dimensional generalized Chebyshev polynomial is used to fit the continuity distribution of the dynamic aerodynamic load on a two-dimensional mapping plane of the fluid grid data according to parameterized fluid grid data, coordinate fitting order and time fitting order are set, and a least square method optimization algorithm is used to calculate a polynomial fitting coefficient matrix; interpolation calculation is performed by using the polynomial fitting coefficient matrix and the parameterized finite element grid data, so that the dynamic distribution load of the finite element grid is obtained; the two divided finite element grid model loads are combined and corresponded with the original finite element grid nodes, the resultant force and the pressure center position of the fluid grid model and the finite element grid model are calculated respectively, and the error is calculated; if the error of the calculation result does not meet the requirement, the fitting order is adjusted and the above fitting steps are repeated; The fluid grid model and the finite element grid model of the airfoil are mapped from three-dimensional space to two-dimensional space, and the process includes: Step 2.

1. The airfoil fluid grid model and the airfoil finite element grid model are divided respectively; there are two ways of division, one is to take three nodes on the airfoil fluid grid model, input the node number, and form a plane according to the coordinates of the three points, and divide the airfoil fluid grid model and the airfoil finite element grid model into two parts by the plane; the other way is to rely on X-Y, X-Z and Y-Z planes, if X-Y plane is selected, the Z-axis coordinate of the plane is input to determine the position of the plane, so as to divide the airfoil grid into two parts; Step 2.

2. The four-edge element of the airfoil fluid grid model and the airfoil finite element grid model is uniformly converted into a triangular element, and the four node numbers of each four-edge element are divided into two groups, each group has three node numbers, and each four-edge element is converted into two new triangular elements; Step 2.

3. The LSCM parameterization method is used to map the airfoil fluid grid model and the airfoil finite element grid model to a two-dimensional plane, and the grid distribution of the airfoil model fluid grid and the airfoil model finite element grid on the mapping plane is drawn respectively, so as to observe the grid quality after parameterization and confirm whether the grid shape is consistent; the specific mapping process is as follows: Assume that each triangle has a local coordinate system (x, y), and the normal of the triangle is along the Z-axis direction in the coordinate system, and the coordinates of the three vertices are (x1, y1), (x2, y2) and (x2, y2); suppose that U:(x, y)→(u, v) is a mapping from three-dimensional space to two-dimensional plane, then in the local coordinate system, the Cauchy-Riemann equation is discretized as: X is represented by a complex number as: X=x+iy (2) For the mapped triangle, also use the local coordinate system, and similarly have: U=u+iv (3) According to the inverse function derivative theorem, it is obtained that: The gradient in the triangle can be described by the following formula: In the formula, A is the area of the triangle; A can be calculated according to the following formula: A=(x1y2-y1x2)+(x2y3-y2x3)+(x2y3-y2x3) (6) Combining the gradient formula in the triangle, we get: Wherein, the value of W1, W2, W3 has nothing to do with the selection of the local coordinate system of the triangle, they are vectors between coordinates, and U j The expression of U is as follows: U j = u j + iv j (j = 1, 2, 3) Because most of the triangular mesh in the parameterization process is not completely consistent before and after the mapping, that is, it cannot completely satisfy formula (4); therefore, a conformal energy C is needed to control the deformation degree of the triangle before and after the mapping, and when the conformal energy C is the smallest, the conformal mapping is as possible to satisfy. The two-dimensional parameterization of the three-dimensional surface is completed by solving the unique minimum value of the conformal energy function; the process of constructing the conformal energy is described as follows: The conformal energy of each triangle is: where A T is the area of the triangle, n is the total number of triangles, and k is considered as the number of the triangle element; In order to facilitate the calculation of the minimum value of the conformal energy, it is converted into a quadratic form, and the total conformal energy of the grid is written as: The matrix M = (M kq ) in the above equation, where: Each triangular mesh node is numbered q = (1, 2, 3...l), l is the total number of nodes; M and U in formula (10) are complex matrices, for a triangular mesh k, three node numbers are q1, q2, q3, then the kth row of the M matrix has corresponding values in the q1, q2, q3 columns The rest of the positions in the row are 0; meanwhile, M and U are written in the following complex form: Where A is the real part matrix of the matrix M, B is the imaginary part matrix of the matrix M, u is the real part vector of U, and v is the imaginary part vector of U; then formula (12) is substituted into formula (10), and the conformal energy is written as: When the quadratic form takes the minimum value, there is a stationary point: ▽C = L C U = 0 (14) Thus far, the matrix L C From the initial node and cell node information, the matrix U containing the mapped node coordinate information is obtained according to equation (14).

5. An electronic device comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that, The processor implements the dynamic aerodynamic load distribution method based on the generalized Chebyshev polynomial according to any one of claims 1 to 3 when executing the computer program. 6.A non-transitory computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to implement the dynamic aerodynamic load distribution method based on the generalized Chebyshev polynomial according to any one of claims 1 to 3.

Images