A two-way fluid-structure interaction flexible flapping wing solving method based on LBM and FEM

Abstract

This invention relates to flapping-wing aerodynamics, specifically to a two-way fluid-structure interaction (FSI) solution method for flexible flapping wings based on LBM and FEM. The two-way FSI method for solving flexible flapping wings includes establishing a dynamic model, simulating the dynamic model using the LBM method (fluidic solver), solving the highly nonlinear transient problems in the dynamic model using the FEM method (structural solids solver), establishing data exchange between the LBM and FEM methods, and performing the FSI solution. The two-way FSI method is used to numerically calculate flapping wings with different Young's moduli, comparing and analyzing the influence of flow field on the flexible flapping wing under different Young's moduli and the influence of flexible flapping wing deformation on the flow field. This invention overcomes the limitation of existing traditional techniques that can only simulate and calculate the flow field under rigid flapping wing motion, thus revealing the influence of flexible flapping wing deformation on the flow field.

Description

Technical Field

[0001] This invention relates to flapping wing aerodynamics, and more specifically to a solution method for flexible flapping wings based on two-way fluid-structure interaction (LBM) and FEM. Background Technology

[0002] In recent years, significant progress has been made both domestically and internationally in the study of the unsteady aerodynamic characteristics of flapping-wing aircraft. This research primarily employs experimental methods or numerical simulations using the Immersed Boundary-Lattice Boltzmann Method (IB-LBM), revealing the dynamic evolution of leading-edge vortices (LEV), wingtip vortices (WTV), and wake vortices (WV) in free-flight states, and their coupling mechanism with lift generation. Current research mainly focuses on the numerical analysis of rigid flapping-wing models. This lack of high-precision aerodynamic calculation methods addresses the impact of deformation on complex flexible structures. Flexible deformation affects the aerodynamic forces of flapping wings, and for flexible flapping wings, severe fluid-structure interaction effects are present. Therefore, studying only rigid flapping wings will result in problems such as low aerodynamic efficiency, preventing them from flying as efficiently as naturally occurring flying organisms.

[0003] The LBM (Liquid Probe) method, a fluid solver, combines the advantages of both microscopic and macroscopic methods. Its simple boundary condition settings allow it to handle complex and irregular structural boundary problems, and it exhibits good parallelism. However, it neglects the large deformation of the flapping fin, limiting its numerical simulation to rigid flapping fins. Furthermore, the large deformation of the flapping fin significantly affects the surrounding flow field distribution, resulting in inaccurate accuracy and lack of universality for flexible flapping fins.

[0004] While the IB-LBM method can effectively handle complex and irregular boundary problems of butterfly-shaped flapping wings by introducing a virtual boundary layer into the mesh and combining the surface geometry of the object with the solution of hydrodynamic equations, it cannot solve flapping wings under large flexible deformations.

[0005] Therefore, it is necessary to design a numerical simulation method for studying the influence of the flow field on the flexible flapping wing under flexible deformation and the influence of the deformation of the flexible flapping wing on the flow field. Summary of the Invention

[0006] The purpose of this invention is to at least solve one of the technical problems existing in the prior art, and to provide a solution method for bidirectional fluid-structure interaction flexible flapping wings based on LBM and FEM.

[0007] To achieve the above objectives, the technical solution adopted by the present invention is as follows: Step 1, Establish a dynamic model: By establishing a quantitative relationship between force and motion, the dynamic behavior of the system under different conditions can be predicted; Step 2: Use the fluid solver LBM method to simulate the dynamic model: Based on the Boltzmann equation, by discretizing the space and velocity space, the macroscopic motion of the continuous medium is described as the motion of discrete particles, thereby simulating the evolution of the fluid in space and time. Step 3: Solve the highly nonlinear transient problem in the dynamic model using the structural solid solver FEM method: Discretize the time into a series of small increment steps and solve using the central difference method; Step 4: Establish a bidirectional data exchange mechanism between the fluid solver LBM and the structural solid solver FEM and perform fluid-structure interaction (FSI) solution: Through the coupling interface and bidirectional iteration, the fluid load is applied to the structure and the structural deformation is fed back to the fluid domain; the fluid solver LBM is set as the master solver, and the synchronization time step of the coupled calculation is negotiated and determined.

[0008] Furthermore, in step 1, the process of establishing the model specifically includes: Step 101, establish a flexible flapping wing model: select a butterfly-shaped flexible flapping wing. The butterfly-shaped flexible flapping wing adopts an isotropic flapping wing configuration. The body and the wing root connection are both set as rigid bodies. The structure of the butterfly-shaped flexible flapping wing consists of a pair of wings and a body. Step 102, mesh generation and boundary conditions setting: periodic boundary conditions are selected, and wall boundary conditions are selected for the remaining lateral and upper and lower end face boundaries. The butterfly-shaped flexible flapping wing adopts a hexahedral structure mesh. Step 103, setting the motion of the butterfly-shaped flexible flapping wings: using Euler angles to describe... coordinate system relative to Rotational motion of the coordinate system, The three unit vectors along the X, Y, and Z axes in a coordinate system are defined as follows: , and ,exist The unit vector in a coordinate system is defined as , and vector group , , A vector group can be obtained through continuous orthogonal transformations. , , as follows:

[0009] ; In the above formula, , , The definition is as follows:

[0010] ;

[0011] ;

[0012] ; In the above formula, This represents the angle of attack of the flapping wing at time t. The flapping amplitude of the wing at time t is defined as follows:

[0013] ;

[0014] ; In the above formula, The maximum angle of attack for flapping wings, The amplitude of wing flapping. The flapping cycle, in which Its magnitude is equal to half of the full-stroke amplitude.

[0015] Furthermore, in step 2, the process of simplifying the model specifically includes: Step 201: Based on the Boltzmann equation, construct an evolution model of the particle distribution function; The evolution formula of the particle distribution function is expressed as:

[0016] ; In the above formula, It is a collision term, velocity distribution function. It is a function of the velocity, space, and time of the object being considered. Represents the acceleration term. Represents spatial coordinates The partial derivatives, Represents the velocity component The partial derivatives, Indicates time The partial derivatives; Step 202, perform a linear approximation on the Boltzmann equation, i.e., the Boltzmann-BGK approximation equation:

[0017] ; In the above formula, It is a relaxation factor It is the equilibrium distribution function; Step 203, to enable the model to correctly recover macroscopic fluid dynamics, uses the following equilibrium distribution function equation:

[0018] ; In the above formula, Represented as weighting coefficients in different discrete velocity directions. Represented as lattice speed of sound, macroscopic density velocity distribution , Indicates different discrete velocity directions. Indicates the first A discrete velocity in directional components, Indicates the first A discrete velocity in A vector of direction; The Kronecker symbol is defined as follows: ={1, α=β0, α≠β}; Step 204: Define an additional matrix containing macroscopic quantities in the fluid system, and assign a relaxation time to each macroscopic quantity. The multi-relaxation-time model formula is as follows:

[0019] ; In the above formula, It is a distribution function. It is a relaxation matrix that controls the relaxation rate of various macroscopic quantities; It is a collision item; Representation of the collision matrix in moment space:

[0020] ; In the above equation, the middle rectangle A rectangle is formed after the collision. , This is the quantity that the velocity distribution function is transformed and mapped to the space matrix. Representing the equilibrium moment function in moment space, as shown in the above equation, the particle distribution function... By transforming the matrix Mapping to a matrix Obtain the transformed matrix , is represented as:

[0021] ; In the above formula, The definition is as follows:

[0022] ; In the above formula, This is called the diagonal collision relaxation matrix. It is a relaxation factor; ; Transformation matrix The matrix has a dimension of 19×19, where each row represents a different physical quantity. The rows of the matrix mainly correspond to the following physical quantities: The first row corresponds to the density. The second row corresponds to the kinetic energy term; the third row corresponds to the higher-order energy term; the fourth to sixth rows respectively represent... The seventh to ninth rows represent the energy flux components. Lines 10 to 15 represent the diagonal and off-diagonal terms of the pressure tensor. The last sixteenth to nineteenth rows represent higher-order stress terms and other fluid properties; through linear combinations of these rows and columns, the matrix... It can map the velocity distribution function of the fluid to the above physical quantities and process them according to the characteristics of the moment; Step 205, improve the computational efficiency of multiple relaxations: by transforming the collision process from velocity space to moment space; matrix It can be represented as:

[0023] ; In the above formula, For fluid density, For energy, The second moment of energy, They represent The momentum component in the direction. Represent directional energy flux component These are the diagonal elements of the pressure tensor. These are the off-diagonal elements of the pressure tensor; Equilibrium Moment Function Related to flux density, energy, and pressure, it can be represented by the following function:

[0024] ; In the above formula, and It is a constant, representing an adjustable parameter, while Let be the velocity vector of the fluid. By processing it in the above formula, we can avoid using the transformation matrix to calculate the velocity vector of discrete particles in each collision.

[0025] Furthermore, step 3 specifically includes: Step 301, the equilibrium equations for the structural nodes are:

[0026] ; In the above formula, Represents the node quality matrix. Indicates acceleration. Indicates external force. Indicates the internal forces within the element; Step 302: Discretize the continuum into a finite number of elements. The displacement within each element is represented by nodal displacements, thus transforming the infinite-degree-of-freedom problem into a finite-degree-of-freedom problem. Then, integrate the velocity to calculate the displacement equation:

[0027] ; In the above formula, Indicates the increment step.

[0028] Furthermore, step 4 specifically includes: Step 401: Use Co-simulation Engine to manage the coupling and data exchange between the fluid solver LBM and the structural solid solver FEM. The Co-simulation Engine program manages the coupling between the fluid solver LBM and the structural solid solver FEM, and is responsible for establishing a data exchange link between the fluid solver LBM and the structural solid solver FEM, while synchronizing the calculation step size of each solver in time. Step 402, execute explicit coupling algorithm and bidirectional data transfer: The solvers use the explicit Gauss-Seidel algorithm for serial coupling calculation, follow the master solver negotiation method, and consider using the time step of the fluid solver LBM as the coupling synchronization benchmark, so that the fluid solver LBM undertakes the master solution task, so as to balance computational efficiency and coupling accuracy. During the coupling process, the fluid solver LBM transfers the stress load at the fluid-structure interaction node to the structural solid solver FEM. After calculating the displacement information after being subjected to the corresponding load, the structural solver FEM transfers the displacement information back to the fluid solver LBM, realizing seamless linkage and bidirectional fluid-structure interaction between the flow field and the structure.

[0029] As can be seen from the above description of the present invention, compared with the prior art, the novel device of the present invention has at least one of the following beneficial effects: (1) Numerical calculations were performed on flapping wings with different Young's moduli using a two-way fluid-structure interaction method combining LBM and FEM. The effects of flow field on flexible flapping wings and the effects of flexible flapping wing deformation on flow field under different Young's moduli were compared and analyzed. The LBM-FEM method can establish data exchange between the fluid solver and the structural mechanics solver and perform two-way fluid-structure interaction solution. The LBM method has simple mesh generation, high parallel efficiency, and can handle complex boundary problems, while the accuracy of FEM in simulating the mechanical deformation of wings has also been verified. (2) The solver uses the explicit Gauss-Seidel algorithm for coupling. The explicit Gauss-Seidel algorithm is a serial process, so during the coupling process, the solution of one physics field will lag behind that of another. The negotiation method of CSE follows the master solver and considers giving priority to using the time step of LBM for data exchange, so that LBM undertakes the master solution task. LBM transfers the force load at the fluid-structure interaction node to FEM and transfers the displacement information after being subjected to the corresponding load from FEM to LBM, realizing bidirectional fluid-structure interaction simulation. Attached Figure Description

[0030] Figure 1 This is a flowchart illustrating the steps in an embodiment of the present invention; Figure 2 This is a schematic diagram of the computational domain and boundary of the butterfly-shaped flexible flapping wing in an embodiment of the present invention; Figure 3 This is a graph showing the deformation displacement of the flexible flapping wingtip in an embodiment of the present invention. Figure 4 This is a schematic diagram comparing the deformation of the rigid flapping wing and the flexible flapping wing during a complete flapping cycle in an embodiment of the present invention; Figure 5 The graph shows the instantaneous lift coefficient curves of the flexible flapping wing in the embodiment of the invention under different Young's moduli. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0032] refer to Figure 1 As shown in the preferred embodiment of the present invention, a solution method for a bidirectional fluid-structure interaction flexible flapping wing based on LBM and FEM includes the following steps: Step 1, Establish a dynamic model: By establishing a quantitative relationship between force and motion, the dynamic behavior of the system under different conditions can be predicted; Step 2: Use the fluid solver LBM method to simulate the dynamic model: Based on the Boltzmann equation, by discretizing the space and velocity space, the macroscopic motion of the continuous medium is described as the motion of discrete particles, thereby simulating the evolution of the fluid in space and time. Step 3: Solve the highly nonlinear transient problem in the dynamic model using the structural solid solver FEM method: Discretize the time into a series of small increment steps and solve using the central difference method; Step 4: Establish a bidirectional data exchange mechanism between the fluid solver LBM and the structural solid solver FEM and perform fluid-structure interaction (FSI) solution: Through the coupling interface and bidirectional iteration, the fluid load is applied to the structure and the structural deformation is fed back to the fluid domain; the fluid solver LBM is set as the master solver, and the synchronization time step of the coupled calculation is negotiated and determined.

[0033] As a preferred embodiment of the present invention, it may also have the following additional technical features: In this embodiment, the process of establishing the model in step 1 specifically includes: Step 101, Establishing the flapping wing model: Due to the intricate and complex structures of the butterfly wing body, such as wing veins, wing membranes, and stigmas, the deformation of the body itself has a very low degree of interference with aerodynamic distribution and unsteady flow states. To broaden the applicability of the simulation conclusions and improve their universality, the original wing body structure needs to be reasonably simplified when establishing the butterfly flapping wing simulation model. An isotropic flapping wing configuration is selected for this study. The body section and the wing root junction are both set to rigid body properties using an isotropic flapping wing configuration. The structure of the butterfly-shaped flexible flapping wing consists of a pair of wings and a body. The structural parameters of the butterfly-shaped flexible flapping wing can be found in the detailed parameter information listed in Table 1.

[0034] Table 1. Parameters of the butterfly-shaped flapping wing structure

[0035]

[0036] Step 102, Mesh generation and boundary conditions setting: In this embodiment, the computational region is divided into two regions. The buffer region has a size of 20L×10L×10L, and the free flight region has a size of 20L×6L×8L. The grid size of the buffer region is set to... The coarse grid, the grid size of the flight area is set to... Fine mesh. The Lagrange points representing the flapping wing boundary are related to the size of the flapping wing body and the left and right wings. Indicates the time step. Lattice spacing.

[0037] The butterfly-shaped flexible flapping wing adopts a hexahedral structure mesh, with periodic boundary conditions selected and wall boundary conditions selected for the remaining lateral and upper / lower end face boundaries. Periodic boundary conditions are primarily used to simulate the infinite extension characteristics of the x-axis transverse space, enabling iterative flow field iteration in the x-axis direction, avoiding lateral boundary interference, and eliminating the influence of transverse boundary effects on the flapping fin flow field. Wall boundary conditions can simulate the closed flow field environment around the flapping fin, constraining the fluid overflow path, ensuring the internal flow field is not disturbed by external flow, and simultaneously conforming to the flow constraint characteristics of real natural and experimental environments, constraining the longitudinal vortex diffusion range. Specific boundary conditions and computational domains are as follows: Figure 2 The computational domain and boundary of the butterfly-shaped flexible flapping wing are shown.

[0038] The butterfly-shaped flexible flapping wing employs a hexahedral mesh with C3D8R elements, each node having three translational degrees of freedom. Volumetric strain is calculated only at a single integration point at the element center. Its core mechanism is to significantly improve computational efficiency while maintaining displacement accuracy and avoiding the numerical limitations of fully integrated elements by reducing the number of integration points. Its advantages include high shear locking efficiency, strong robustness against large deformations, low computational cost, and good convergence. This paper considers isotropic linear elastic materials, discretized using 2372 elements.

[0039] The body and wings are connected in a Cartesian + Cardan configuration to allow for sufficient translational and rotational degrees of freedom. The butterfly-shaped flexible flapping wing torso only participates in the translation of the entire fuselage. While the wings pitch around the center point of the body, they also flap around the wing roots on both sides of the body. Therefore, both the torso and the wing roots are designed as rigid bodies.

[0040] Step 103, setting the motion of the butterfly-shaped flexible flapping wings: using Euler angles to describe... coordinate system relative to Rotational motion of the coordinate system, The three unit vectors along the X, Y, and Z axes in a coordinate system are defined as follows: , and , and The unit vector in a coordinate system is defined as , and vector group , , A vector group can be obtained through continuous orthogonal transformations. , , as follows:

[0041] ; in , , The definition is as follows:

[0042] ;

[0043] ;

[0044] ; In the above formula, This represents the angle of attack of the flapping wing at time t. The flapping amplitude of the wing at time t is defined as follows:

[0045] ;

[0046] ; In the above formula, The maximum angle of attack for flapping wings, The amplitude of wing flapping. The flapping cycle, in which Its magnitude is equal to half of the full-stroke amplitude.

[0047] The initial motion state of free flight in this paper is set according to the research of Kosuke Suzuki et al., and flapping wings are selected. The trajectory of motion when =16.

[0048] In this embodiment, step 2, the process of simplifying the model specifically includes: The basic idea of ​​LBM is to establish a simplified dynamic model. Based on the Boltzmann equation, by discretizing the space and velocity space, the macroscopic motion of the continuous medium is described as the motion of discrete particles, thereby simulating the evolution process of the fluid in space and time.

[0049] Step 201: Based on the Boltzmann equation, construct an evolution model of the particle distribution function; The evolution formula of the particle distribution function is expressed as:

[0050] ; In the above formula, It is a collision term, velocity distribution function. It is a function of the velocity, space, and time of the object being considered. Represents the acceleration term. Represents spatial coordinates The partial derivatives, Represents the velocity component The partial derivatives, Indicates time The partial derivatives; in Representing the acceleration term, a first-order Taylor expansion of the above equation yields:

[0051] ; After the collision, the above equation becomes:

[0052] ; Step 202: Apply a linear approximation to the Boltzmann equation, i.e., the Boltzmann-BGK approximation equation:

[0053] ;

[0054] ; In the above formula, It is a relaxation factor It is the equilibrium distribution function; This method is essentially a simulation process in which statistical mesoscopic particles collide and migrate on a discrete grid, and the flow effects of the fluid are described accordingly. During the simulation, the focus is primarily on the evolution of the particle velocity distribution function in various directions. Through this simulation and analysis, the flow characteristics of particles during collisions and migrations can be obtained, further enhancing our understanding of the fluid's motion.

[0055] According to the definition of gas dynamics, the macroscopic density of the fluid can be obtained from the evolved distribution function. velocity distribution And three macroeconomic quantities:

[0056] ;

[0057] ;

[0058] ; Step 203, to enable the model to correctly recover macroscopic fluid dynamics, uses the following equilibrium distribution function equation: The collision term in the equation is:

[0059] ;

[0060] ; In the above formula, Represented as weighting coefficients in different discrete velocity directions. Represented as lattice speed of sound, macroscopic density velocity distribution , Indicates different discrete velocity directions. Indicates the first A discrete velocity in directional components, Indicates the first A discrete velocity in A vector of direction; The Kronecker symbol is defined as follows: ={1, α=β0, α≠β}; Step 204: Define an additional matrix that includes macroscopic quantities in the fluid system, and assign a relaxation time to each macroscopic quantity. The multi-relaxation-time model formula is as follows:

[0061] ; The collision process formula is as follows:

[0062] ; In the above formula, This refers to the spatial position vector and time, where r represents the space of the grid node. In two dimensions, r = (x, y), and in three dimensions, r = (x, y, z). The formula for the coordinate migration process is as follows:

[0063] ; In the above formula, This refers to the time step. Representing the local collision process in vector form:

[0064] ; In the above formula, It is a distribution function. It is a relaxation matrix that controls the relaxation rate of various macroscopic quantities; It is a collision item; Representation of the collision matrix in moment space:

[0065] ; In the above equation, the middle rectangle A rectangle is formed after the collision. , This is the quantity that the velocity distribution function is transformed and mapped to the space matrix. Representing the equilibrium moment function in moment space, as shown in the above equation, the particle distribution function... By transforming the matrix Mapping to a matrix Obtain the transformed matrix , is represented as:

[0066] ; In the above formula, The definition is as follows:

[0067] ; In the above formula, This is called the diagonal collision relaxation matrix. It is a relaxation factor; ; Transformation matrix The matrix has a dimension of 19×19, where each row represents a different physical quantity. The rows of the matrix mainly correspond to the following physical quantities: The first row corresponds to the density. The second row corresponds to the kinetic energy term; the third row corresponds to the higher-order energy term; the fourth to sixth rows respectively represent... The seventh to ninth rows represent the energy flux components. Lines 10 to 15 represent the diagonal and off-diagonal terms of the pressure tensor. The last sixteenth to nineteenth rows represent higher-order stress terms and other fluid properties; through linear combinations of these rows and columns, the matrix... It can map the velocity distribution function of the fluid to the above physical quantities and process them according to the characteristics of the moment; Step 205, improve the computational efficiency of multiple relaxations: by transforming the collision process from velocity space to moment space; matrix It can be represented as:

[0068] ; In the above formula, For fluid density, For energy, The second moment of energy, They represent The momentum component in the direction. Represent directional energy flux component These are the diagonal elements of the pressure tensor. These are the off-diagonal elements of the pressure tensor; Equilibrium Moment Function Related to flux density, energy, and pressure, it can be represented by the following function:

[0069] ; In the above formula, and It is a constant, representing an adjustable parameter, while Let be the velocity vector of the fluid. By processing it in the above formula, we can avoid using the transformation matrix to calculate the velocity vector of discrete particles in each collision.

[0070] In this embodiment, step 3 specifically includes: Step 301, the equilibrium equations for the structural nodes are:

[0071] ; In the above formula, Represents the node quality matrix. Indicates acceleration. Indicates external force. Indicates the internal forces within the element; According to the above equation, in the increment step The nodal acceleration at time is:

[0072] ; The solver uses a diagonal mass matrix, so there is no need to solve a system of equations. If the acceleration is constant in a very short time, the central difference method can be used to integrate the acceleration and obtain the nodal velocity.

[0073] ; Step 302: Discretize the continuum into a finite number of elements. The displacement within each element is represented by nodal displacements, thus transforming the infinite-degree-of-freedom problem into a finite-degree-of-freedom problem. Then, integrate the velocity to calculate the displacement equation:

[0074] ; In the above formula, Indicates the increment step.

[0075] In this embodiment, step 4 specifically includes: The bidirectional coupling between the fluid solver LBM and the structural solid solver FEM has the strongest compatibility, especially suitable for fluid-structure interaction simulation of flexible flapping wings. The fluid solver LBM uses a fixed Eulerian grid to simulate the flow field, without the need for grid reconstruction as the solid structure deforms. It describes the geometry and deformation of the flexible flapping wing through Lagrangian, which can accurately capture the microscopic vortex structure (such as leading edge vortex and wake vortex) of unsteady flow field at low Reynolds number, and can flexibly adapt to the large-amplitude motion and large-scale passive deformation of the flexible flapping wing. FEM can accurately solve the mechanical response of flexible structures and capture the passive torsion and camber deformation of the airfoil. The two-way coupling can achieve seamless linkage between the flow field and the structure.

[0076] Step 401: Use Co-simulation Engine to manage the coupling and data exchange between the fluid solver LBM and the structural solid solver FEM. The Co-simulation Engine program manages the coupling between the fluid solver LBM and the structural solid solver FEM, and is responsible for establishing a data exchange link between the fluid solver LBM and the structural solid solver FEM, while synchronizing the calculation step size of each solver in time. Step 402, execute explicit coupling algorithm and bidirectional data transfer: The solvers use the explicit Gauss-Seidel algorithm for serial coupling calculation, follow the master solver negotiation method, and consider using the time step of the fluid solver LBM as the coupling synchronization benchmark, so that the fluid solver LBM undertakes the master solution task, so as to balance computational efficiency and coupling accuracy. During the coupling process, the fluid solver LBM transfers the stress load at the fluid-structure interaction node to the structural solid solver FEM. After calculating the displacement information after being subjected to the corresponding load, the structural solver FEM transfers the displacement information back to the fluid solver LBM, realizing seamless linkage and bidirectional fluid-structure interaction between the flow field and the structure.

[0077] In this embodiment, the method further includes numerical calculation and verification: In 2010, Aono et al. verified the aeroelastic response of a hovering Zimmermann wing under single-DOF rotation, using an isotropic Zimmermann wing flapping at 10 Hz in incompressible, viscous flow. The numerical results of Aono et al. are compared with the numerical results obtained from two-way fluid-structure interaction solutions to verify the reliability of the method.

[0078] The functional expression for the airfoil is as follows:

[0079] ;

[0080] ; In the above formula, The chord length at the wing root is 0.025m. The maximum span of the airfoil is 0.075m. Airfoil thickness... It has a diameter of 0.0004m, is made of aluminum plate, and has an airfoil flapping angle of 0.0004m. The pattern of change over time is as follows ,in For amplitude, The flapping frequency is 10 Hz. Specific structural parameters of the Zimmermann airfoil are shown in Table 2, and relevant flow parameters are shown in Table 3, including the airfoil's structural parameters and related numerical simulation parameters. The Reynolds number is... Defined as ,in is the kinetic viscosity coefficient, and the reference velocity is the average flapping velocity at the wingtip.

[0081] Table 2 Zimmermann wing structure parameters

[0082]

[0083] Table 3 Relevant Flow Parameters

[0084]

[0085] By comparing the results calculated using the numerical method in this embodiment with those obtained by Aono et al., it was found that the differences between the two sets of results are small, and the overall trend is consistent. Within the error range, the bidirectional fluid-structure interaction numerical simulation method for flexible flapping wings in this embodiment can be considered reliable.

[0086] In this embodiment, specifically, a method for numerically monitoring the deformation of a flexible flapping wing is also included. A monitoring point A is placed at the wingtip to monitor wingtip deformation. Since the bending deformation of the forewing is more pronounced than that of the hindwing, two monitoring points B and D, and C and E, are placed near the middle and tip of the forewing span, respectively, to monitor the bending deformation. Substituting the coordinates of monitoring points B and D, C and E on the flexible flapping wing forewing with the coordinates of the rigid flapping wing forewing into the bending angle calculation formula yields the bending deformation at the middle and wingtip of the butterfly-shaped flapping wing. The bending angle β is calculated as follows:

[0087]

[0088] In the above formula, , , This represents the difference in coordinates between two points along the spanwise direction of the rigid flapping wing. , , The coordinate difference between two points corresponding to the spanwise position of the flexible flapping wing can be calculated to obtain the bending angle of the two parts at that moment, reflecting the overall bending state of the forewing.

[0089] Based on the monitored values, the deformation values ​​of the wingtip monitoring point A during one cycle in free flight under different Young's moduli were obtained. The variation curves are shown below. Figure 3 As shown, the curves exhibit a generally periodic pattern. The flexible wing reaches its minimum and maximum wingtip deformation at two moments near the top and bottom dead centers (around T=0.25 and T=0.75, respectively). 'a' represents the minimum value, and 'b' represents the maximum value; specific values ​​are shown in Table 4. As the Young's modulus increases, the difference between 'b' and 'a' gradually decreases. This is because the wing consistently generates lift during flight. Within a reasonable range, a smaller Young's modulus results in greater lift. However, the aerodynamic drag generated by flapping is in the opposite direction to the flapping motion. Furthermore, it is also affected by elastic and inertial forces, leading to changes in the forces acting on the wings during flapping.

[0090] Table 4. Deformation results and differences of the upper and lower dead centers at the tip of the flexible flapping wing under different Young's moduli.

[0091]

[0092] In this embodiment, an aerodynamic analysis method is also included. like Figure 4 As shown, green represents a flexible flapping wing; flexible flapping wings with Young's moduli of 90 MPa and 110 MPa were selected for comparison with rigid flapping wings. The chordal torsional deformation and spanwise bending deformation of the butterfly-shaped wings during flapping can be clearly seen in the figures; the flapping phase of the flexible flapping wing is delayed compared to the rigid flapping wing. Combined with... Figure 3 The wingtip deformation displacement analysis shown shows that during the time period T=0-0.5, the flexible flapping wing is positioned below the rigid flapping wing due to deformation; at T=0.75, the flapping wing flaps upward, and the flexible flapping wing is positioned above the rigid flapping wing.

[0093] In this embodiment, flexible flapping wings with Young's modulus of 70 MPa, 90 MPa, 110 MPa and 130 MPa are used as the analysis objects, and other parameters are set as follows: flapping angle of attack 40°, flapping amplitude 45° and flapping frequency 5 Hz.

[0094] like Figure 4 This is a comparison chart of lift coefficients under different Young's moduli. The gray shaded area represents the descent phase, and the white area represents the ascent phase. The entire motion process is as follows: starting from the horizontal plane, descent to the lowest point (descent phase); starting from the lowest point, ascent to the highest point (ascent phase); and descent from the highest point back to the horizontal plane (descent phase).

[0095] like Figure 5 As shown, when T=0 is in the middle of the descent phase, the lift coefficient of the flexible flapping wing is significantly higher than that of the rigid flapping wing, indicating that flexibility does indeed have a positive effect on the lift of the flapping wing. At this point, it can be observed that as the Young's modulus decreases, the more flexible the wing material, the greater the corresponding lift coefficient, with the flexible flapping wing with a Young's modulus of 70 MPa exhibiting the largest lift coefficient. As the flapping wing continues to descend, the lift coefficient of each flapping wing continuously decreases, and the smaller the Young's modulus, the faster the lift coefficient decreases. Around T=0.18, the lift coefficient curves intersect at a point, and negative lift begins to occur.

[0096] When 0.25 < T < 0.75 and the flapping wing is in the upstroke stage, the entire upstroke, whether it is a flexible flapping wing or a rigid flapping wing, is in a negative lift state. For the flexible flapping wing, when T is approximately 0.4, the lift coefficient drops to the lowest, showing a valley point; while the time for the rigid flapping wing to generate the valley point of the lift coefficient is a bit later, around T = 0.5, with a phase difference. And as the Young's modulus decreases, the valley point of the negative lift becomes lower. This is because the large deformation of the flexible structure exacerbates the flow separation, resulting in a greater decrease in lift; although the flexibility brings a greater positive lift, it also brings a greater negative lift. As the flapping wing continues to move up, the lift coefficient begins to rise, but it is still in the negative lift state.

[0097] The upstroke of the flapping wing ends and it enters the downstroke stage 0.75 < T < 1. The flapping wing begins to generate a positive lift coefficient, and the lift coefficient curve climbs rapidly. The smaller the Young's modulus, the faster the climbing speed. The flexible flapping wing with a Young's modulus of 70 mpa reaches the peak of the positive lift coefficient at approximately T = 0.9, and the peak value is much larger than the peak values generated by other flapping wings. For the flapping wings with other Young's moduli, the lift coefficient reaches the peak at T = 1. This shows that the large deformation of the flapping wing will cause a phase difference, and the flow separation will also occur earlier accordingly.

[0098] As described above, it is only the preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention, according to the technical solution of the present invention and its improved concept, making equivalent substitutions or changes, should be covered by the protection scope of the present invention.

Claims

1. A two-way fluid-structure interaction (LSM) solution method for flexible flapping wings based on LBM and FEM, used to solve the LSM problem of flexible flapping wings, characterized in that, Includes the following steps: Step 1, Establish a dynamic model: By establishing a quantitative relationship between force and motion, the dynamic behavior of the system under different conditions can be predicted; Step 2: Use the fluid solver LBM method to simulate the dynamic model: Based on the Boltzmann equation, by discretizing the space and velocity space, the macroscopic motion of the continuous medium is described as the motion of discrete particles, thereby simulating the evolution of the fluid in space and time. Step 3: Solve the highly nonlinear transient problem in the dynamic model using the structural solid solver FEM method: Discretize time into a series of small increment steps using the explicit time integration method, and solve using the central difference method; Step 4: Establish a bidirectional data exchange mechanism between the fluid solver LBM and the structural solid solver FEM and perform fluid-structure interaction (FSI) solution: Through the coupling interface and bidirectional iteration, the fluid load is applied to the structure and the structural deformation is fed back to the fluid domain; the fluid solver LBM is set as the master solver, and the synchronization time step of the coupled calculation is negotiated and determined.

2. The solution method for a bidirectional fluid-structure interaction flexible flapping wing based on LBM and FEM according to claim 1, characterized in that, In step 1, the process of building the model specifically includes: Step 101, establish a flexible flapping wing model: select a butterfly-shaped flexible flapping wing. The butterfly-shaped flexible flapping wing adopts an isotropic flapping wing configuration. The body and the wing root connection are both set as rigid bodies. The structure of the butterfly-shaped flexible flapping wing consists of a pair of wings and a body. Step 102, mesh generation and boundary conditions setting: periodic boundary conditions are selected, and wall boundary conditions are selected for the remaining lateral and upper and lower end face boundaries. The butterfly-shaped flexible flapping wing adopts a hexahedral structure mesh. Step 103, setting the motion of the butterfly-shaped flexible flapping wings: using Euler angles to describe... coordinate system relative to Rotational motion of the coordinate system, The three unit vectors along the X, Y, and Z axes in a coordinate system are defined as follows: , and ,exist The unit vector in a coordinate system is defined as , and vector group , , A vector group can be obtained through continuous orthogonal transformations. , , as follows: ； In the above formula, , , The definition is as follows: ； ； ； In the above formula, This represents the angle of attack of the flapping wing at time t. The flapping amplitude of the wing at time t is defined as follows: ； ； In the above formula, The maximum angle of attack for flapping wings, The amplitude of wing flapping. The flapping cycle, in which Its magnitude is equal to half of the full-stroke amplitude.

3. The solution method for a bidirectional fluid-structure interaction flexible flapping wing based on LBM and FEM according to claim 1, characterized in that, In step 2, the process of simplifying the model specifically includes: Step 201: Based on the Boltzmann equation, construct an evolution model of the particle distribution function; The evolution formula of the particle distribution function is expressed as: ； In the above formula, It is a collision term, velocity distribution function. It is a function of the velocity, space, and time of the object being considered. Represents the acceleration term. Represents spatial coordinates The partial derivatives, Represents the velocity component The partial derivatives, Indicates time The partial derivatives; Step 202, perform a linear approximation of the Boltzmann equation, i.e., the Boltzmann-BGK approximation equation: ； In the above formula, It is a relaxation factor It is the equilibrium distribution function; Step 203, to enable the model to correctly recover macroscopic fluid dynamics, uses the following equilibrium distribution function equation: ； In the above formula, Represented as weighting coefficients in different discrete velocity directions. Represented as lattice speed of sound, macroscopic density velocity distribution , Indicates different discrete velocity directions. Indicates the i-th discrete velocity at directional components, Indicates the first A discrete velocity in A vector of direction; The Kronecker symbol is defined as follows: ={1, α=β0, α≠β}; Step 204: Define an additional matrix that includes macroscopic quantities in the fluid system, and assign a relaxation time to each macroscopic quantity. The multi-relaxation-time model formula is as follows: ； In the above formula, It is a distribution function. It is a relaxation matrix that controls the relaxation rate of various macroscopic quantities; It is a collision item; Representation of the collision matrix in moment space: ； In the above equation, the middle rectangle A rectangle is formed after the collision. , This is the quantity that the velocity distribution function is transformed and mapped to the space matrix. Representing the equilibrium moment function in moment space, as shown in the above equation, the particle distribution function... By transforming the matrix Mapping to a matrix Obtain the transformed matrix , represented as: ； In the above formula, The definition is as follows: ； In the above formula, This is called the diagonal collision relaxation matrix. It is a relaxation factor; ; Transformation matrix The matrix has a dimension of 19×19, where each row represents a different physical quantity. The rows of the matrix mainly correspond to the following physical quantities: The first row corresponds to the density. The second row corresponds to the kinetic energy term; the third row corresponds to the higher-order energy term; the fourth to sixth rows respectively represent... The seventh to ninth rows represent the energy flux components. Lines 10 to 15 represent the diagonal and off-diagonal terms of the pressure tensor. The last sixteenth to nineteenth rows represent higher-order stress terms and other fluid properties; through linear combinations of these rows and columns, the matrix... It can map the velocity distribution function of the fluid to the above physical quantities and process them according to the characteristics of the moment; Step 205, improve the computational efficiency of multiple relaxations: by transforming the collision process from velocity space to moment space; matrix It can be represented as: ； In the above formula, For fluid density, For energy, The second moment of energy, They represent The momentum component in the direction. Represent Energy flux component in the direction These are the diagonal elements of the pressure tensor. These are the off-diagonal elements of the pressure tensor; Equilibrium Moment Function Related to flux density, energy, and pressure, it can be represented by the following function: ； In the above formula, and It is a constant, representing an adjustable parameter, while Let be the velocity vector of the fluid. By processing it in the above formula, we can avoid using the transformation matrix to calculate the velocity vector of discrete particles in each collision.

4. The solution method for a bidirectional fluid-structure interaction flexible flapping wing based on LBM and FEM according to claim 1, characterized in that, Step 3 specifically includes: Step 301, the equilibrium equations for the structural nodes are: ； In the above formula, Represents the node quality matrix. Indicates acceleration. Indicates external force. Indicates the internal forces within the element; Step 302: Discretize the continuum into a finite number of elements. The displacement within each element is represented by nodal displacements, thus transforming the infinite-degree-of-freedom problem into a finite-degree-of-freedom problem. Then, integrate the velocity to calculate the displacement equation: ； In the above formula, Indicates the increment step.

5. The solution method for a bidirectional fluid-structure interaction flexible flapping wing based on LBM and FEM according to claim 1, characterized in that, Step 4 specifically includes: Step 401: Use Co-simulation Engine to manage the coupling and data exchange between the fluid solver LBM and the structural solid solver FEM. The Co-simulation Engine program manages the coupling between the fluid solver LBM and the structural solid solver FEM, and is responsible for establishing a data exchange link between the fluid solver LBM and the structural solid solver FEM, while synchronizing the calculation step size of each solver in time. Step 402, execute explicit coupling algorithm and bidirectional data transfer: The solvers use the explicit Gauss-Seidel algorithm for serial coupling calculation, follow the master solver negotiation method, and consider using the time step of the fluid solver LBM as the coupling synchronization benchmark, so that the fluid solver LBM undertakes the master solution task, so as to balance computational efficiency and coupling accuracy. During the coupling process, the fluid solver LBM transfers the stress load at the fluid-structure interaction node to the structural solid solver FEM. After calculating the displacement information after being subjected to the corresponding load, the structural solver FEM transfers the displacement information back to the fluid solver LBM, realizing seamless linkage and bidirectional fluid-structure interaction between the flow field and the structure.

Images