Intelligent unmanned aerial vehicle pneumatic structure reduced order analysis model

Abstract

An intelligent unmanned aerial vehicle pneumatic structure reduced-order analysis model comprises a pneumatic reduced-order model and a structure reduced-order model. The aerodynamic structure order reduction analysis process comprises the steps that an aerodynamic order reduction model and a structure order reduction model are initialized; under the given calculation working condition, the aerodynamic load acting on the intelligent unmanned aerial vehicle structure is obtained through calculation of the aerodynamic reduced-order model; transmitting the load on the pneumatic grid to the structural grid through a force interpolation algorithm; solving deformation response by the structure order reduction model under the action of aerodynamic load; the structural deformation is fed back to the pneumatic grid through displacement interpolation, and the pneumatic configuration is updated; and performing iterative calculation until a convergence criterion is met, finally obtaining a pneumatic structure analysis result in an equilibrium state, and outputting a pneumatic coefficient and structural deformation and stress. According to the method, the CFD / CSD model is replaced by the equivalent model, so that the calculation efficiency is remarkably improved on the premise of ensuring the calculation precision.

Description

Technical Field

[0001] This invention relates to the field of aerodynamic structure analysis technology for aircraft, specifically a reduced-order analysis model for the aerodynamic structure of an intelligent unmanned aerial vehicle. Background Technology

[0002] The CFD / CSD coupled method offers advantages such as high accuracy and high confidence in calculating aeroelasticity. With the rapid development of computer technology, its application areas are becoming increasingly widespread. However, the pace of computer development cannot keep up with the demands of fluid-structure interaction researchers, and the computational cost of CFD methods is increasing. Therefore, from the perspectives of aerodynamics and structures, establishing physically meaningful reduced-order aerodynamic structural analysis models to achieve rapid calculations of lift, drag, and structural response, thereby significantly reducing computational complexity, is of great significance. Summary of the Invention

[0003] In view of the increasing computational overhead of current CFD / CSD coupled methods, which makes it difficult to meet the requirements of rapid calculation of lift, drag and structural response of aircraft, this invention proposes a reduced-order analysis model of aerodynamic structure for intelligent unmanned aerial vehicles.

[0004] The technical solution of this invention is as follows:

[0005] A reduced-order aerodynamic structure analysis model for intelligent unmanned aerial vehicles, including a reduced-order aerodynamic model and a reduced-order structural model;

[0006] The aerodynamic structure order reduction analysis process is as follows:

[0007] Initialize the aerodynamic reduced-order model and the structural reduced-order model;

[0008] Under given computational conditions, the aerodynamic reduced-order model calculates the aerodynamic loads acting on the structure of the intelligent unmanned aerial vehicle; the loads on the aerodynamic mesh are transferred to the structural mesh through a force interpolation algorithm.

[0009] Then, the reduced-order structural model is used to solve the deformation response under aerodynamic loads; and the structural deformation is fed back to the aerodynamic mesh through displacement interpolation to update the aerodynamic configuration.

[0010] By iteratively calculating until the convergence criterion is met, the aerodynamic structure analysis results under equilibrium state are finally obtained, and the aerodynamic coefficients, structural deformation, and stress are output.

[0011] A further preferred embodiment is that the method for constructing the aerodynamic order reduction model is as follows:

[0012] First, reasonable sampling data is obtained in the design space using the maximum area adaptive sampling method of triangulation.

[0013] Secondly, the POD basis containing the features of the original data is obtained through the intrinsic orthogonal decomposition method;

[0014] Ultimately, the modeling and training of the aerodynamic reduced-order model were achieved.

[0015] A further preferred embodiment of the triangulation maximum area adaptive sampling method includes the following process:

[0016] ① Provide the initial sampling points using any sampling method;

[0017] ② Perform Delaunay triangulation on the existing sampling points;

[0018] ③ Calculate the area of ​​all triangles, mark the triangle with the largest area and add a sampling point at its centroid; and add the triangle side length constraint to the judgment of the sampling point location. When the sum of the side lengths of two sides is less than the side length of the third side, the location of the new sampling point is defined as the midpoint of the third side.

[0019] ④ Repeat steps ② and ③ until the required total number of sampling points is reached.

[0020] A further preferred embodiment, wherein obtaining the POD basis containing the features of the original data through the intrinsic orthogonal decomposition method includes the following process:

[0021] Considering the discrete compressible Navier-Stokes equations, for steady-state problems, removing the time term yields the residual governing equations of the following form:

[0022]

[0023] In the formula, These are the conserved variables of the computational domain boundary and the grid cells. w = [ r , r u , r v , r w , r E , ⋯ ] T The dimension of variable x is determined by the number of design parameters. x = [ AOA , Yes , H , ⋯ ] T ;

[0024] The flow field variables are represented using the POD basis as follows:

[0025]

[0026] In the formula It is composed of the first n POD bases, and the subspace spanned by it is called the trial space. The result is obtained using the POD base coefficients. Low-dimensional control equations with variables as variables:

[0027]

[0028] This equation has N nonlinear equations and n unknown POD basis coefficients; a set of reduced-order basis coefficients is selected. The vector field spanned is called the test space. Projecting the control equations onto the test space yields the reduced-order model equations:

[0029]

[0030] Solving using Newton's iteration method

[0031]

[0032] In the formula, For each Newton iteration step, K is determined by the convergence criterion. This is the initial guess value. These represent the increment of the base coefficients and the Jacobian matrix at the k-th iteration step, respectively. Is The step size in the search direction can be a unit length. Initial value It can be determined by a given initial flow field ,Right now .

[0033] After iterative convergence, the flow field solution obtained by applying the POD reduction method to the governing equations is:

[0034] .

[0035] A further optimized approach modifies the process of obtaining the POD basis containing the features of the original data through the intrinsic orthogonal decomposition method:

[0036] Will Perform a Taylor expansion at step k, neglecting higher-order error terms, to obtain...

[0037]

[0038] To minimize the above expression, we define a sum of squares function.

[0039]

[0040] It also introduces constraint terms and constraint residuals for boundary conditions.

[0041]

[0042] In the formula, These represent the increment of the flow field variables and the deviation between the predicted and accurate values ​​in the boundary region at the k-th iteration, respectively; the sum of squares functions are weighted and combined into a minimization problem.

[0043]

[0044] Will This can be written in the form of a weighted matrix to obtain a modified sum of squares function.

[0045]

[0046] in, For constraint terms, the weight matrix The contribution of each grid cell to the constraint term is determined, and its form is as follows:

[0047]

[0048] In the formula, This refers to the number of boundary grid cells, including the incoming flow boundary and the object surface boundary. To calculate the maximum distance between all grid cells and the boundary within the domain, This represents the total number of grid cells; for the weight matrix, it is corrected using a Gaussian weighting function.

[0049]

[0050] in

[0051]

[0052] In the formula, This is a scaling factor related to the grid cell position y. All are greater than 0;

[0053] Obtain the constrained iterative formula

[0054] .

[0055] In a further preferred embodiment, the structural order reduction model uses a polynomial to condense the displacements of the finite element nodes on the cross section of the full-order model to the centroid of the cross section, thereby establishing a reduced basis matrix for model order reduction.

[0056] A further preferred embodiment is that the displacement of any finite element node i on section j is expressed as:

[0057] u i j = { u ix j u iy j u iz j } = [ F t F t F t ] { u ¯ cx t j u ¯ cy t j u ¯ Czech t j } = R i j u ¯ c j , t = 1 , 2 , ⋯ , W

[0058] in, Let i represent the displacement vector of finite element node i on section j. This represents the displacement transformation matrix of node i in the finite element method. Let represent the generalized displacement vector of the centroid of the section containing node i in the finite element method. Representing an Nth-order polynomial:

[0059]

[0060] in Let be the number of terms in the polynomial, and x and y represent the differences in the x-direction and y-direction coordinates between the finite element node i and the centroid of its cross section, respectively.

[0061] A further optimized approach involves dividing the structure into several beam segments along the axial direction. Using the centroids of the cross-sections at both ends of each beam segment as principal nodes, a displacement transformation relationship between the finite element nodes and the principal nodes of each beam segment is established using polynomials and beam element shape functions, thereby reducing the model order.

[0062] The expression for the displacement transformation relationship between finite element node i and the centroid j of its cross section within a beam segment is as follows:

[0063] u i j = { u ix j u iy j u iz j } = [ F t F t F t ] { u ¯ cx t j u ¯ cy t j u ¯ Czech t j } = R i j u ¯ c j , t = 1 , 2 , ⋯ , W

[0064] The displacement of centroid j can be approximated by the displacement of the main nodes of the beam segment, where the axial coordinates of main nodes k and k+1 are respectively... and The axial coordinate of centroid j is The displacement of the centroid is then written as:

[0065]

[0066] in, and It is the shape function of the beam element.

[0067] [ S k S k + 1 ] = [ z k + 1 − z j z k + 1 − z k z j − z k z k + 1 − z k ]

[0068] The displacement transformation relationship between finite element node i and the master node within the beam segment is as follows:

[0069]

[0070] Assuming there are m finite element nodes in the k-th beam segment, the displacement transformation relationship of the m nodes is as follows:

[0071]

[0072] Assuming the entire structure is divided into q beam segments along the axial direction, the displacement transformation equation of the structure is:

[0073] U = { U 1 ⋮ U k ⋮ U q } = [ T ˜ 1 ⋱ T ˜ k ⋱ T ˜ q ] { u ˜ 1 u ˜ 2 ⋮ u ˜ k u ˜ k + 1 ⋮ u ˜ q u ˜ q + 1 } = T ˜ N q U ˜

[0074] in, This represents the displacement vector of all nodes in the finite element method. Represents q beam segments and N The reduced basis matrix is ​​established by a polynomial of order 1.

[0075] A further preferred approach involves modifying the reduced basis matrix as follows:

[0076]

[0077] Where K and M are the stiffness matrix and mass matrix of the full-order model, respectively.

[0078] Beneficial effects

[0079] This invention proposes a reduced-order aerodynamic structure analysis model for intelligent unmanned aerial vehicles (UAVs), employing an equivalent model to replace the CFD / CSD model, significantly improving computational efficiency while maintaining computational accuracy. Examples demonstrate that the time-domain propulsion fluid-structure interaction calculation method used in this invention possesses high reliability.

[0080] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0081] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the description of the embodiments taken in conjunction with the following drawings, in which:

[0082] Figure 1 Flowchart of the aerodynamic POD reduced-order model;

[0083] Figure 2 Initial sampling point distribution map;

[0084] Figure 3 Delaunay triangulation diagram;

[0085] Figure 4 Distribution map after marking the maximum area;

[0086] Figure 5 Delaunay triangulation diagram after adding new sampling points;

[0087] Figure 6 : Delaunay triangulation diagram of sampling points in the DMAS method after adding boundary corner points. Solid blue dots represent original sampling points, solid red dots represent newly added sampling points, and the sequence number indicates the order in which sampling points are added (corresponding iteration step); (when adding a sampling point, the cyan-blue shading represents the boundary triangle).

[0088] Figure 7 : Delaunay triangulation diagram of sampling points in DMAS method after adding boundary corner points. Solid blue dots represent original sampling points, and solid red dots represent newly added sampling points. The index indicates the order in which the sampling points are added (corresponding to the iteration step); (10 sampling points added);

[0089] Figure 8 : Schematic diagram of Galerkin-POD method error;

[0090] Figure 9: A schematic diagram of polynomial-based model order reduction;

[0091] Figure 10 : A schematic diagram of polynomial-based model order reduction;

[0092] Figure 11 Diagram showing the relationship between main nodes and centroid within a beam segment;

[0093] Figure 12 The number of degrees of freedom in the reduced-order model;

[0094] Figure 13 MAC similarity assessment;

[0095] Figure 14 Aerodynamic structure reduced-order analysis model;

[0096] Figure 15 Geometric model parameters;

[0097] Figure 16 Finite element computational mesh;

[0098] Figure 17 Schematic diagram of the first-order mode shape of the AGARD 445.6 wing, frequency f=9.2808Hz;

[0099] Figure 18 Schematic diagram of the second-order mode shape of the AGARD 445.6 wing, frequency f=38.014Hz;

[0100] Figure 19 : Schematic diagram of the third-order mode shape of the AGARD 445.6 wing, frequency f=47.853Hz;

[0101] Figure 20 Schematic diagram of the fourth mode shape of the AGARD 445.6 wing, frequency f=91.753Hz;

[0102] Figure 21 : AGARD 445.6 wing Ma=0.901 structural displacement response history. Detailed Implementation

[0103] The embodiments of the present invention are described in detail below. These embodiments are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0104] This embodiment proposes a reduced-order aerodynamic structure analysis model for intelligent unmanned aerial vehicles, based on aerodynamics and structural science, specifically including:

[0105] 1. Aerodynamic calculation model and order reduction method

[0106] This invention uses CFD methods to calculate sample points, and then establishes a POD degradation model to calculate the aerodynamic data required for the loitering configuration of intelligent unmanned aerial vehicles.

[0107] 1.1 CFD Computational Model

[0108] The Navier-Stokes equations are used to solve for the physical quantities of fluid particles at different locations and times within a flow field. They include three main governing equations: mass, energy, and Newton's second law (momentum). Each governing equation can be written in three-dimensional Cartesian coordinates as follows:

[0109]

[0110] in For time; , , These are the three direction axes of a three-dimensional rectangular coordinate system; As a conserved variable, , , For inviscid flux; , , This is the viscous flux.

[0111] Conservation variable matrix Inside, in order, are: mass density, momentum density in three directions, and energy density.

[0112] Q = [ r , r u , r v , r w , r e ] T

[0113]

[0114] Inviscid flux matrix , , The contents, in order, are: mass flow rate, momentum and pressure from the inflowing mass in three directions, energy from the inflowing mass and work done by pressure; viscous flux matrix. , , Inside, in order, are: the impulses in three directions provided by viscous force, the heat input due to thermal conduction, and the work done by viscous force. Their specific forms are as follows:

[0115] E = [ r u r u 2 + p r UV r oh ( r e + p ) u ] , E v = [ 0 t xx t xyz t xz u t xx + v t xyz + w t xz − q ˙ x ]

[0116] F = [ r v r wow r v 2 + p r VW ( r e + p ) v ] , F v = [ 0 t x t yyy t yz u t x + v t yyy + w t yz − q ˙ y ]

[0117] G = [ r w r wu r wv r w 2 + p ( r e + p ) w ] , G v = [ 0 t zx t zy t zzz u t zx + v t zy + w t zzz − q ˙ z ]

[0118] The expression for the shear stress term in the formula is as follows:

[0119]

[0120] in

[0121]

[0122]

[0123]

[0124]

[0125] In the formula, The laminar viscosity coefficient can be obtained from the Sutherland formula; The turbulent viscosity coefficient is given by the turbulence model; the values ​​of the constants for the gas are as follows: For specific heat ratio, ; For free-flowing Mach number; Let be the Prandtl constant for laminar flow, and take . ; Let be the Prandtl constant for turbulence, and take... .

[0126] Turbulence equations are adopted The two-equation model is as follows:

[0127] ∂ ( r k ) ∂ t + u i ∂ ( r k ) ∂ x i = P k − β k r k oh + ∂ ∂ x i [( m l + m t s k ) ∂ k ∂ x i ] ∂ ( I am ) ∂ t + u i ∂ ( I am ) ∂ x i = C oh P oh − β oh I am 2 + ∂ ∂ x i [( m l + m t s oh ) ∂ oh ∂ x i ] + 2 r ( 1 − E ) 1 s oh 2 1 oh ∂ k ∂ x i ∂ oh ∂ x i

[0128] Among them, the generating term and the dissipating term are respectively

[0129] P k = m l [( ∂ u i ∂ x j )( ∂ u i ∂ x j + ∂ u j ∂ x i ) − 2 3 ( ∂ u k ∂ x k ) 2 ] P oh = gr [( ∂ u i ∂ x j )( ∂ u i ∂ x j + ∂ u j ∂ x i ) − 2 3 ( ∂ u k ∂ x k ) 2 ]

[0130] Turbulent viscosity coefficient is defined as

[0131]

[0132] The SSTk-ω turbulence model performs well across various flow types, including low Reynolds number flows, flows with adverse pressure gradients, boundary layer flows, and separated flows. It exhibits high accuracy in boundary layer simulation and relatively good computational efficiency, thus finding wide application in numerous fields such as aerospace, mechanical engineering, and automotive engineering. This invention targets an intelligent unmanned aerial vehicle employing a coupled-wing design, with a Reynolds number typically on the order of Re = 1.07 × 10⁻⁶. 6 Since it is a turbulent flow, this invention uses the SSTk-ω turbulent model to solve the RANS equations, where the time-progression scheme adopts implicit time-progression.

[0133] 1.2 Aerodynamic order reduction model

[0134] Supported by CFD numerical simulation technology, data analysis-based order reduction models are an effective way to achieve rapid flow field prediction. This invention analyzes the experimental design method for flow field data sampling, combines intrinsic orthogonal decomposition (POD) dimensionality reduction technology, establishes a rapid numerical prediction method for flow fields, develops a parallel computational order reduction model framework suitable for complex flow fields, and conducts research on rapid numerical prediction of large-scale flow fields that meets engineering application requirements.

[0135] Data-driven approximation models are an effective approach for solving time-critical engineering problems. These methods aim to generate a low-dimensional, low-complexity computational model that captures as closely as possible the relationship between system inputs (such as design variables and shape parameters) and the outputs of high-fidelity model calculations (such as lift / drag, flow field information, and performance indicators). To reduce model dimensionality, these methods first perform an "offline" phase, acquiring flow data within a specific design space through high-fidelity CFD numerical simulation to calculate a low-dimensional subspace that captures the model's main physical behaviors. Then, in the "online" prediction phase, the optimal solution is searched within this low-dimensional subspace. For example... Figure 1 As shown.

[0136] This invention addresses the following aspects of approximate models based on data analysis: First, it obtains reasonable sampling data in the design space through an adaptive sampling method using triangulation, i.e., it uses high-fidelity numerical simulation to obtain a certain amount of sampling data to capture as much information as possible from the design space; second, it obtains a POD basis containing the features of the original data through the Proper Orthogonal Decomposition (POD) method; and finally, it realizes the modeling and training of the approximate model.

[0137] 1) Adaptive sampling method for maximum area of ​​triangulation

[0138] Delaunay triangulation with large-area adaptive sampling is an extremely important preprocessing technique for both numerical analysis and computer graphics. It provides a special definition for triangular meshes, possessing the excellent properties of empty circles and maximizing the minimum angle, thus avoiding the generation of elongated triangles.

[0139] The minimum angle maximization property of Delaunay triangulation means that among the interior angles that can form a triangular network from the same set of scattered points, the minimum interior angle of the triangular network formed by Delaunay triangulation is the largest. The empty circle property means that the circumcircle of the triangle formed by Delaunay triangulation does not contain a fourth data point.

[0140] The DMAS steps are as follows:

[0141] ① Provide the initial sampling points using any sampling method (such as LHS). For example... Figure 2 Ten sampling points were generated in the [0,1] design space using the LHS method, and are represented by red dots.

[0142] ② Perform Delaunay triangulation on the existing sampling points. For example... Figure 3 .

[0143] ③ Calculate the area of ​​all triangles, and mark the triangle with the largest area, adding a sampling point at its centroid. For example... Figure 4 Newly added sampling points are indicated by red pentagrams. Figure 5 The new Delaunay triangulation after adding new sampling points is given.

[0144] ④ Repeat steps ② and ③ until the required total number of sampling points is reached.

[0145] As can be seen from the above, the DMAS method can only sample within the initial sampling points, therefore it is necessary to manually add reasonable boundary sampling points, which makes the DMAS method less concise. If 2 are added initially... D If there are several boundary corner points, the final sampling result will always appear as follows: Figure 6 and Figure 7 As shown, the triangles obtained by the Delaunay boundary partitioning, with the boundary design variable range as their sides, will result in elongated triangles, such as... Figure 6 The triangle is marked with a medium blue shading. This invention incorporates the triangle side length constraint into the determination of the location of the additional sampling point. That is, when the sum of the lengths of two sides is less than a certain threshold compared to the length of the third side, the location of the new sampling point is defined as the midpoint of the third side, which can effectively solve this problem.

[0146] 2) POD order reduction method for flow control equations

[0147] Projecting the governing equations and searching for the optimal flow field solution that satisfies them in the subspace spanned by the POD basis allows for consideration of the actual flow physics, resulting in a more accurate representation of the actual flow conditions. Through this process, the originally complex nonlinear partial differential equations can be transformed into a set of nonlinear ordinary differential equations for solving the basis coefficients, making them easier to solve. In unsteady conditions, time is typically used as a design variable to generate sampled solutions; in steady conditions, time is disregarded, and factors such as angle of attack, Mach number, rudder deflection, and geometry are used as design variables to generate sampled solutions. Considering the discrete compressible Navier-Stokes equations as follows, removing the time term for steady problems yields the residual governing equations of the following form.

[0148]

[0149] where is the boundary of the computational domain, and the conserved variables of the grid cell w = [ r , r u , r v , r w , r E , ⋯ ] T , the dimension of the variable x is determined by the number of design parameters, x = [ AOA , Yes , H , ⋯ ] T

[0150] The flow field variables are represented in terms of the POD basis as:

[0151]

[0152] where is composed of the first n POD bases, and the spanned subspace is called the trial space, obtaining the low-dimensional variable control equation with the POD basis coefficients as variables:

[0153]

[0154] This equation has N non-linear equations and n unknowns (POD basis coefficients). Generally, the number of bases used is much smaller than the degrees of freedom of the initial problem (n << N), that is, the above equation is an overdetermined system of equations, and no POD basis coefficients can make it hold exactly, but the residual of the reduced-order control equation can be made as small as possible.

[0155] Select a set of reduced-order bases , and the spanned vector field is called the test space , project the control equation onto this test space to obtain the reduced-order model equation:

[0156]

[0157] Since the residual R is non-linear, it can be solved by the Newton iteration method

[0158]

[0159] where is the Newton iteration step, K is determined by the convergence criterion, is the initial guess value, are the increments of the basis coefficients and the Jacobian matrix at the k-th iteration respectively, is at the step size in the search direction, which can take the unit length. The initial value can be obtained from the given initial flow field , that is .

[0160] After the iteration converges, the flow field solution obtained by applying the POD reduced-order method of the control equation is

[0161]

[0162] The above-mentioned method for reducing the order of the governing equations using POD is the Galerkin-POD method.

[0163] To select the optimal test space, a qualitative analysis of the error of the Galerkin-POD method is performed. The error of Galerkin-POD is defined as the accurate value of the flow field. The value obtained from the reduced subspace difference

[0164]

[0165] It can be decomposed into two mutually orthogonal errors.

[0166]

[0167] In the formula, This indicates that we are currently in a testing phase. The error component, Indicates parallel to The error components, and their error relationships are as follows: Figure 8 As shown:

[0168] If it were merely a data compression problem, the flow field solution would reside in a low-dimensional trial space. Approximately, the vertical error When the residual control equations are projected onto the test space... When this happens, additional errors are generated. Since the perpendicular error is orthogonal to the test space, it is only related to the test space; while the parallel error is also related to the test space.

[0169] To minimize or even eliminate the parallelism error, i.e., to obtain the optimal test space, when the Jacobian matrix is ​​a symmetric positive definite matrix, Galerkin projection can be applied to determine the test basis. This is widely used in structural mechanics; for the Jacobian matrix in the usual case, such as in fluid mechanics, generally speaking... Since it is not a symmetric positive definite matrix, applying the least squares Petrov-Galerkin projection, we take... ,Right now

[0170]

[0171] Mathematically, the above equation is equivalent to solving it using the Gauss-Newton method. It should be noted that the boundary conditions must be met during the iterative process; therefore, the Galerkin-POD method does not handle the boundary conditions separately during iterative solving.

[0172] 1.3 Constraint Modification of the Galerkin-POD Method

[0173] To address the convergence issue encountered in the Galerkin-POD method, constraint terms are added. To add appropriate constraints, we consider the mathematical principles of the Gauss-Newton iteration method. Perform a Taylor expansion at step k, neglecting higher-order error terms, to obtain...

[0174]

[0175] To minimize the above expression, we define a sum of squares function.

[0176]

[0177] By taking the partial derivative of the above equation and setting the result to 0, it can be seen that treating boundary elements and mesh elements within the computational domain of the flow field equally during the iteration process fails to distinguish between boundary and internal elements. This may lead to the calculated flow field deviating from the constraints of the boundary conditions, thus converging to an incorrect solution. Therefore, this invention introduces a constraint term, namely the constraint residual of the boundary conditions.

[0178]

[0179] In the formula, Let represent the increment of the flow field variables at the k-th iteration and the deviation between the predicted and accurate values ​​in the boundary region, respectively. These can be approximated, and the two minimization problems can be weighted and then merged into a single minimization problem.

[0180]

[0181] Will This can be written in the form of a weighted matrix to obtain a modified sum of squares function.

[0182]

[0183] in, For constraint terms, the weight matrix This determines the contribution of each grid cell to the constraint term. Its form is as follows:

[0184]

[0185] In the formula, This refers to the number of boundary grid cells, including the incoming flow boundary and the object surface boundary. To calculate the maximum distance between all grid cells and the boundary within the domain, This represents the total number of grid cells. This invention provides an improved Gaussian weighting function, defined as...

[0186]

[0187] In the formula, This is a scaling factor related to the grid cell position y. All values ​​are greater than 0, affecting the shape of the weight function. Therefore, the weight matrix is ​​corrected as follows:

[0188]

[0189] This yields the constrained iterative formula.

[0190] .

[0191] 2. Structural calculation model and order reduction method

[0192] 2.1 CSD Calculation Model

[0193] The wing undergoes deformation and motion under the influence of dynamic aerodynamic forces, and its governing equations are as follows:

[0194]

[0195] in, , , These are the overall mass matrix, overall damping matrix, and overall stiffness matrix, which describe the mass distribution, drag characteristics, and elastic deformation of the wing, respectively; F is the external load borne by the structure. , , These represent nodal acceleration, velocity, and displacement, respectively.

[0196] By employing the Newmark implicit time integration method to solve the structural dynamics control equations, the velocity and acceleration within one time step are first expressed as:

[0197]

[0198] Where Δt is the time step, These are the integration parameters of Newmark, which determine the integration accuracy and stability. In this embodiment, the parameters selected for calculation are... =0.167, =0.5, at which point the Newmark integral has second-order accuracy.

[0199] Further, by organizing, we can obtain information about The system of equations is as follows:

[0200]

[0201] By solving The solution can be obtained from the above formula. , This means that the solution to the structural dynamics equations is completed within one time step, and this solution process is repeated until the final time step.

[0202] 2.2 Structural order reduction model

[0203] By using polynomials to condense the displacements of finite element nodes on the cross section of the full-order model to the centroid of that cross section, a reduced basis matrix for model order reduction is quickly established. For structures with arbitrary meshes or complex geometries, beam element shape functions are introduced to further improve the polynomial-based model order reduction method, solving the problems that the polynomial-based model order reduction method requires the full-order model to have a regular mesh and the reduced model to be too large.

[0204] This invention establishes the displacement transformation relationship between finite element nodes and the centroid of any cross section using polynomials. Figure 9 A schematic diagram of the order reduction model is given, where yellow dots represent the centroids of the cross sections, and red dots represent the finite element nodes on the cross sections. The basic idea of ​​this order reduction method is to approximate the displacements of the finite element nodes on each cross section as a linear combination of a set of polynomials, and the combination coefficients are the generalized degrees of freedom of the reduced model. The displacement of any finite element node i on cross section j can be expressed as:

[0205] u i j = { u ix j u iy j u iz j } = [ F t F t F t ] { u ¯ cx t j u ¯ cy t j u ¯ Czech t j } = R i j u ¯ c j , t = 1 , 2 , ⋯ , W

[0206] in, Let i represent the displacement vector of finite element node i on section j. This represents the displacement transformation matrix of node i in the finite element method. Let represent the generalized displacement vector of the centroid of the section containing node i in the finite element method. Representing an Nth-order polynomial, the polynomial functions are shown in Table 1 below, where... Let be the number of terms in the polynomial, and x and y represent the differences in the x-direction and y-direction coordinates between the finite element node i and the centroid of its cross section, respectively.

[0207] Table 1 Polynomial Functions

[0208]

[0209] The number of degrees of freedom of the reduced-order model established by the above polynomials is related to the number of cross sections p of the full-order model. When the number of cross sections along the axial direction is large, the number of degrees of freedom of the reduced-order model will still be large. At the same time, the model reduction method requires the structure to have a regular mesh and also needs to be modified for more complex configurations. Furthermore, we divide the structure into several beam segments along the axial direction, with the centroids of the cross sections at both ends of the beam segment as the main nodes. We use polynomials and beam element shape functions to establish the displacement transformation relationship between the finite element nodes in each beam segment and the main nodes of the beam segment, thereby realizing the model reduction. Figure 10 A schematic diagram of model order reduction based on polynomials and beam element shape functions is given, in which yellow dots represent the centroid of the cross section, red dots represent the finite element nodes on the cross section, and green dots represent the main nodes of the beam segment.

[0210] The expression for the displacement transformation relationship between finite element node i and the centroid j of its cross section within a beam segment is as follows:

[0211] u i j = { u ix j u iy j u iz j } = [ F t F t F t ] { u ¯ cx t j u ¯ cy t j u ¯ Czech t j } = R i j u ¯ c j , t = 1 , 2 , ⋯ , W

[0212] The displacement of centroid j can be approximated by the displacement of the main nodes of the beam segment, where the axial coordinates of main nodes k and k+1 are respectively... and The axial coordinate of centroid j is ,like Figure 11 As shown.

[0213] The displacement of the centroid can then be written as:

[0214]

[0215] in, and It is the shape function of the beam element.

[0216] [ S k S k + 1 ] = [ z k + 1 − z j z k + 1 − z k z j − z k z k + 1 − z k ]

[0217] The displacement transformation relationship between finite element node i and the master node within the beam segment is as follows:

[0218]

[0219] Assuming there are m finite element nodes in the k-th beam segment, the displacement transformation relationship of the m nodes is as follows:

[0220]

[0221] Assuming the entire structure is divided into q beam segments along the axial direction, the displacement transformation equation of the structure is:

[0222] U = { U 1 ⋮ U k ⋮ U q } = [ T ˜ 1 ⋱ T ˜ k ⋱ T ˜ q ] { u ˜ 1 u ˜ 2 ⋮ u ˜ k u ˜ k + 1 ⋮ u ˜ q u ˜ q + 1 } = T ˜ N q U ˜

[0223] in, This represents the displacement vector of all nodes in the finite element method. Represents q beam segments and N The reduced basis matrix is ​​established using a polynomial of order 1. The degrees of freedom of the reduced-order model are as follows: Figure 12 As shown.

[0224] Depend on Figure 12 It can be seen that the number of degrees of freedom in the reduced-order model is related to the number of beam segments q and the order of the polynomial N. Since the number of beam segments q is much smaller than the number of sections p in the finite element model, the improved reduced-order model has fewer degrees of freedom and higher computational efficiency. To ensure the computational accuracy of the reduced-order model, the reduced basis matrix is ​​further modified as follows:

[0225]

[0226] Where K and M are the stiffness matrix and mass matrix of the full-order model, respectively.

[0227] 2.3 Similarity Assessment of Structural Reduction Models

[0228] In the process of model reduction, it is necessary not only to reduce the analysis scale of the structural dynamic model, but also to ensure that the reduced model has good similarity with the original full-order model. Commonly used indicators for evaluating the prediction accuracy of the reduced model include: frequency error between the reduced model and the full-order model, modal similarity between the reduced model and the full-order model, etc.

[0229] The non-uniqueness of the MAC values ​​for commonly used repetition modes makes traditional modal confidence criteria unsuitable for assessing the similarity of structural repetition modes. MAC similarity assessment, as follows... Figure 13 As shown.

[0230] If there exists a combination coefficient such that the MAC value is not less than 0.9, then the similarity between the two modes is relatively good. Therefore, this invention proposes an improved modal confidence criterion to calculate the maximum value of the MAC value, and then uses the iMAC value to evaluate the similarity between the repetitive frequency modes. The expression is as follows:

[0231]

[0232] when At that time, the maximum value of MAC was

[0233]

[0234] when At that time, the maximum value of MAC was

[0235]

[0236] Based on the above formula, the linear combination coefficients of the maximum MAC value are obtained as follows:

[0237]

[0238]

[0239] The maximum value is:

[0240]

[0241] Considering both of the above situations, It can be represented as:

[0242]

[0243] when When the reduced-order model mode and the full-order model mode have good similarity, when At that time, the modal similarity between the reduced-order model and the full-order model is poor.

[0244] 3. Fluid-structure interaction analysis

[0245] In the multidisciplinary optimization design of intelligent unmanned aerial vehicles (UAVs), fluid-structure interaction exists. The analysis of these two coupled physical fields is called fluid-structure interaction (FSI). FSI is divided into unidirectional and bidirectional FSI. Unidirectional FSI does not consider the changes in the solid structure caused by the fluid acting on the solid, thus affecting the flow field. It usually does not require excessive hardware resources when calculating unidirectional FSI. However, bidirectional FSI considers the changes in the flow field caused by changes in the solid structure at each calculation step, often requiring huge computational resources.

[0246] Intelligent unmanned aerial vehicles (UAVs) generate distributed aerodynamic loads in the flow field, which cause structural deformation when acting on the structure. The main deformable component of intelligent UAVs, especially projectile-launched vehicles, is the wing, which typically has high stiffness and low wing loading, so the impact of structural deformation on the flow field can be disregarded. Therefore, this invention employs a unidirectional fluid-structure interaction research method that balances computational efficiency and accuracy, realistically simulating the aerodynamic forces and structural deformation of intelligent UAVs.

[0247] The one-way coupling simulation approach of this invention is as follows:

[0248] Step 1: In the flow field calculation, the aerodynamic load of the intelligent unmanned aerial vehicle is predicted by the aerodynamic reduced-order model.

[0249] Step 2: Then, the predicted aerodynamic loads are transferred to the intelligent unmanned aerial vehicle structure using the fluid-structure interaction surface.

[0250] Step 3: After the intelligent unmanned aerial vehicle is subjected to aerodynamic loads, the response of the solid structure is calculated using a structural reduced-order model to obtain the deformation and stress results of the intelligent unmanned aerial vehicle.

[0251] Ensuring the accuracy of data mapping is one of the main issues that needs to be considered. Commonly used data interpolation methods have four main requirements: 1) ensuring smooth interpolation; 2) accurately converting information between two sets of grids; 3) being able to complete interpolation even when the grid is sparse; and 4) being applicable to complex shapes.

[0252] There are many methods for transferring information at coupling interfaces, but they can be broadly categorized into two types: local interpolation and global interpolation. Local interpolation refers to the method where, during mesh data exchange, different fluid and solid meshes can be interpolated through partial elements or known physical quantities existing at mesh points. However, local interpolation has significant drawbacks, especially when dealing with complex curved surface boundaries, as it excessively increases computational complexity and makes the process cumbersome. Common local interpolation methods include weighted residual methods, constant volume conversion methods, and mapping interpolation. Global interpolation, on the other hand, requires knowledge of all known physical quantities of all meshes or elements to calculate the mapping values ​​between different meshes during mesh data transfer. Its advantage is that it eliminates the need for searching for local interpolation points, saving computational effort. Its disadvantage is that it is prone to significant fluctuations when there are near points in the interpolation points, especially with traditional spline and spline surface fitting methods. Global interpolation methods include various spline function methods, radial basis function methods, and spline surface interpolation methods.

[0253] This paper applies morphological interpolation theory to large deformation structures / aerodynamic coupling, transferring motion information from the structural interface to the aerodynamic interface. Given the coordinates of n structural nodes and the deformation vector... It is necessary to obtain the deformation vector corresponding to m aerodynamic mesh nodes. Based on the coordinates of the structural nodes, the equation for the interpolation function coefficients is as follows:

[0254]

[0255] In the formula This is a coefficient matrix containing the coordinates of the structural nodes. The interpolation coefficient matrix is ​​undetermined. It is a coefficient matrix containing the coordinates of the aerodynamic grid nodes.

[0256] The displacement vector at the aerodynamic grid node can be expressed as:

[0257]

[0258] In the formula This is a coefficient matrix containing the coordinates of the aerodynamic mesh nodes. This is the displacement interpolation matrix between structural nodes and aerodynamic nodes.

[0259] Using the principle of virtual work, the force interpolation relationship between the aerodynamic model and the structural model can be obtained as follows:

[0260]

[0261] In the formula The load is at the aerodynamic interface. The load is on the structural interface.

[0262] 4. Aerodynamic structure reduced-order analysis model

[0263] Combining the aforementioned aerodynamic / structural reduced-order model and aerodynamic / structural coupling method, an efficient aeroelastic analysis workflow is constructed: First, the aerodynamic and structural models are initialized. Under a given computational condition, the aerodynamic reduced-order model calculates the aerodynamic loads acting on the structure. The loads on the aerodynamic mesh are transferred to the structural mesh using a force interpolation algorithm. The structural reduced-order model solves for the deformation response under the aerodynamic loads. Subsequently, the structural deformation is fed back to the aerodynamic model through displacement interpolation to update the aerodynamic configuration. The system iterates until the convergence criterion is met, ultimately obtaining the aeroelastic analysis results under equilibrium conditions, outputting lift / drag coefficients, structural deformation, and stress. This method uses an equivalent model to replace the CFD / CSD model, significantly improving computational efficiency while ensuring computational accuracy. The specific algorithm flow is as follows: Figure 14 As shown.

[0264] This embodiment uses the AGARD 445.6 standard model to verify the time-domain two-way fluid-structure interaction analysis method. Referring to the transonic wind tunnel flutter experiment at NASA Langley Research Center, the wing material properties are shown in Table 2, and the geometric parameters are as follows: Figure 15 As shown.

[0265] Table 2 Material Properties of AGARD 445.6

[0266]

[0267] First, the modal frequencies of the AGARD 445.6 wing configuration established in this paper were verified to ensure that the structural characteristics of its experimental model remained consistent. The computational mesh was generated as follows: Figure 16 As shown in Table 3, the comparison between the calculated frequencies of the first four modalities and the experimental values ​​reveals that the error range for the first four modalities is only 0.227% to 3.317%, and their mode shapes are within... Figure 17~Figure 20 The pattern of "one bend, two twists, three bends, and four twists" is consistent with the experimental results.

[0268] Table 3 Comparison of the first four modal frequencies of AGARD 445.6

[0269]

[0270] Furthermore, the flutter boundary of the established model was obtained using a fluid-structure interaction calculation method. To maintain consistency with the experimental conditions, the dimensionless velocity calculation method is shown in the following equation, where... For dimensionless velocity, For the incoming flow velocity, For reference length, For first-order torsional frequency, For mass ratio, For wing mass, For incoming flow density, This refers to volume parameters.

[0271]

[0272] Then, based on the experimental results, dimensionless velocities around the flutter critical velocity can be selected, and the time history of the structural displacement response can be calculated using the time-domain propagation method. Figure 21 The structural displacement response history of the AGARD 445.6 wing at Ma=0.901 shows that the CFD / CSD fluid-structure interaction method used in this invention calculates the structural displacement response with high accuracy, closely matching experimental data. Table 4 compares the calculated flutter critical velocity with experimental values, showing that the relative error is within 2.6% in the range of Ma=0.499~0.901. In conclusion, the time-domain propulsion fluid-structure interaction calculation method used in this invention possesses high reliability.

[0273] Table 4. Comparison of experimental results and calculation results of flutter boundary in AGARD 445.6 standard mode with those in this paper.

[0274]

[0275] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention without departing from the principles and spirit of the present invention.

Claims

1. A reduced-order aerodynamic structure analysis model for an intelligent unmanned aerial vehicle, characterized in that: This includes aerodynamic order reduction models and structural order reduction models; The aerodynamic structure order reduction analysis process is as follows: Initialize the aerodynamic reduced-order model and the structural reduced-order model; Under given computational conditions, the aerodynamic reduced-order model calculates the aerodynamic loads acting on the structure of the intelligent unmanned aerial vehicle; the loads on the aerodynamic mesh are transferred to the structural mesh through a force interpolation algorithm. Then, the reduced-order structural model is used to solve the deformation response under aerodynamic loads; and the structural deformation is fed back to the aerodynamic mesh through displacement interpolation to update the aerodynamic configuration. By iteratively calculating until the convergence criterion is met, the aerodynamic structure analysis results under equilibrium state are finally obtained, and the aerodynamic coefficients, structural deformation, and stress are output.

2. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 1, characterized in that: The method for constructing the aerodynamic order reduction model is as follows: First, reasonable sampling data is obtained in the design space using the maximum area adaptive sampling method of triangulation. Secondly, the POD basis containing the features of the original data is obtained through the intrinsic orthogonal decomposition method; Ultimately, the modeling and training of the aerodynamic reduced-order model were achieved.

3. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 2, characterized in that: The triangulation maximum area adaptive sampling method The process includes the following: ① Provide the initial sampling points using any sampling method; ② Perform Delaunay triangulation on the existing sampling points; ③ Calculate the area of ​​all triangles, mark the triangle with the largest area and add a sampling point at its centroid; and add the triangle side length constraint to the judgment of the sampling point location. When the sum of the side lengths of two sides is less than the side length of the third side, the location of the new sampling point is defined as the midpoint of the third side. ④ Repeat steps ② and ③ until the required total number of sampling points is reached.

4. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 2, characterized in that: The process of obtaining the POD basis containing the features of the original data through the intrinsic orthogonal decomposition method includes the following steps: Considering the discrete compressible Navier-Stokes equations, for steady-state problems, removing the time term yields the residual governing equations of the following form: In the formula, These are the conserved variables of the computational domain boundary and the grid cells. The dimension of variable x is determined by the number of design parameters. ; The flow field variables are represented using the POD basis as follows: In the formula It is composed of the first n POD bases, and the subspace spanned by it is called the trial space. The result is obtained using the POD base coefficients. Low-dimensional control equations with variables as variables: This equation has N nonlinear equations and n unknown POD basis coefficients; a set of reduced-order basis coefficients is selected. The vector field spanned is called the test space. Projecting the control equations onto the test space yields the reduced-order model equations: Solving using Newton's iteration method In the formula, For each Newton iteration step, K is determined by the convergence criterion. This is the initial guess value. These represent the increment of the base coefficients and the Jacobian matrix at the k-th iteration step, respectively. Is Step size in the search direction; After iterative convergence, the flow field solution obtained by applying the POD reduction method to the governing equations is: 。 5. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 4, characterized in that: The process of obtaining the POD basis containing the features of the original data through the intrinsic orthogonal decomposition method is modified as follows: Will Perform a Taylor expansion at step k, neglecting higher-order error terms, to obtain... To minimize the above expression, we define a sum of squares function. It also introduces constraint terms and constraint residuals for boundary conditions. In the formula, These represent the increment of the flow field variables and the deviation between the predicted and accurate values ​​in the boundary region at the k-th iteration, respectively; the sum of squares functions are weighted and combined into a minimization problem. Will This can be written in the form of a weighted matrix to obtain a modified sum of squares function. in, For constraint terms, the weight matrix The contribution of each grid cell to the constraint term is determined, and its form is as follows: In the formula, This refers to the number of boundary grid cells, including the incoming flow boundary and the object surface boundary. To calculate the maximum distance between all grid cells and the boundary within the domain, This represents the total number of grid cells; for the weight matrix, it is corrected using a Gaussian weighting function. in In the formula, This is a scaling factor related to the grid cell position y. All are greater than 0; Obtain the constrained iterative formula 。 6. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 1, characterized in that: The structural order reduction model uses polynomials to condense the displacements of finite element nodes on the cross section of the full-order model to the centroid of the cross section, thereby establishing a reduced basis matrix for model order reduction.

7. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 6, characterized in that: The displacement of any finite element node i on section j is expressed as: in, Let i represent the displacement vector of finite element node i on section j. This represents the displacement transformation matrix of node i in the finite element method. Let represent the generalized displacement vector of the centroid of the section containing node i in the finite element method. Representing an Nth-order polynomial: in Let be the number of terms in the polynomial, and x and y represent the differences in the x-coordinate and y-coordinate of the finite element node i and the centroid of its cross section, respectively.

8. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 7, characterized in that: The structure is divided into several beam segments along the axial direction. Using the centroids of the two ends of each beam segment as principal nodes, a displacement transformation relationship between the finite element nodes and the principal nodes of each beam segment is established using polynomials and beam element shape functions, thus reducing the model order. The expression for the displacement transformation relationship between finite element node i and the centroid j of its cross section within a beam segment is as follows: The displacement of centroid j can be approximated by the displacement of the main nodes of the beam segment, where the axial coordinates of main nodes k and k+1 are respectively... and The axial coordinate of centroid j is The displacement of the centroid is then written as: in, and It is the shape function of the beam element. The displacement transformation relationship between finite element node i and the master node within the beam segment is as follows: Assuming there are m finite element nodes in the k-th beam segment, the displacement transformation relationship of the m nodes is as follows: Assuming the entire structure is divided into q beam segments along the axial direction, the displacement transformation equation of the structure is: in, This represents the displacement vector of all nodes in the finite element method. Represents q beam segments and N The reduced basis matrix is ​​established by a polynomial of order 1.

9. The reduced-order analysis model of the aerodynamic structure of an intelligent unmanned aerial vehicle according to claim 8, characterized in that: The reduced basis matrix is ​​modified as follows: Where K and M are the stiffness matrix and mass matrix of the full-order model, respectively.

Images