A numerical simulation method for fluid-structure interaction based on lattice Boltzmann flux algorithm

Abstract

The present invention discloses a fluid-solid coupling numerical simulation method based on the lattice Boltzmann flux algorithm, and the method steps are as follows: S1, determine the calculation area, read the grid information, physical parameters and control parameters; S2, solve the Navier-Stokes equations of the two-dimensional unsteady flow field, wherein the convection term is solved using the lattice Boltzmann method, and then the flow field parameters are updated; S3, solve the motion equation; S4, update the flow field boundary conditions and grid coordinates according to the result of solving the motion equation in step S3; S5, judge whether the simulation time is reached, if yes, then enter step S7, otherwise, then enter step S6; S6, return to step S2 to perform the next time step calculation; S7, simulation ends. The method of the present invention applies the lattice Boltzmann flux method to the inviscid flux solution of the unsteady flow field, which can improve the capture of the detailed characteristics of the flow field from a physical level, thereby effectively improving the aerodynamic simulation capability.

Description

Technical Field

[0001] The invention relates to the technical field of fluid mechanics, in particular to a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm. Background Art

[0002] Fluid-Structure Interaction (FSI) is an independent branch of fluid mechanics that intersects solid mechanics. It primarily studies the deformation and dynamic properties of solids under the influence of fluids, as well as the effects of motion on the fluid. This research began with the aeroelasticity of aircraft wings, initially focusing on calculating the flutter boundaries of wings. It has since expanded to include broader fields such as water conservancy, oceanography, aerospace, and medical treatment.

[0003] Common flow field numerical simulation methods, such as the finite difference method (FDM) and the finite volume method (FVM), currently divide the cell boundary flux into two components: inviscid and viscous. The calculation of the inviscid flux component often uses numerical schemes based on the characteristics of the Euler equation and often treats the cell interface as a discontinuity (e.g., Riemann solver). The solution to the viscous flux generally uses central differences based on the continuity assumption of interface physical quantities. Furthermore, due to the inherent inadequacy of the Euler equation in describing the actual physical evolution of the interface, both finite difference and finite volume method-based schemes have specific numerical flaws and often require manual "correction."

[0004] The aforementioned flow field solution technology presents at least one problem: the calculation of flow field parameters is often constrained by the continuity assumptions inherent in the equations. Because these assumptions deviate from the actual situation, the simulated flow field can deviate significantly from the actual flow field, leading to significant discrepancies between the calculated flow field parameters and the true values, which limits its practical application. Summary of the Invention

[0005] The purpose of the present invention is to provide a fluid-solid coupling numerical simulation method based on the lattice Boltzmann flux algorithm to address the problems existing in the prior art.

[0006] The purpose of the present invention is to be solved by the following technical solutions:

[0007] A fluid-structure coupling numerical simulation method based on a lattice Boltzmann flux algorithm is characterized in that the method comprises the following steps:

[0008] S1, determine the calculation area, read the grid information, physical parameters and control parameters, and enter step S2;

[0009] S2. Solve the Navier-Stokes equations for the two-dimensional unsteady flow field, where the convection term is solved using the lattice Boltzmann method and the viscous term is solved using the SA turbulence model, then update the flow field parameters and proceed to step S3;

[0010] S3, solve the motion equation and go to step S4;

[0011] S4. Update the flow field boundary conditions and grid coordinates based on the results of solving the motion equations in step S3 (the displacement caused by the object movement requires updating the grid coordinates, and the surface velocity and acceleration caused by the object movement will affect the boundary conditions of the flow field calculation), and then proceed to step S5;

[0012] S5, determine whether the simulation time has been reached, if yes, proceed to step S7, if not, proceed to step S6;

[0013] S6, return to step S2 to perform the next time step calculation;

[0014] S7. Simulation ends.

[0015] The grid information in step S1 includes: the number of grid blocks, the number of grids in each grid block, the connection information of the grid blocks and the grid coordinates; the physical property parameters include: the density, velocity, temperature, pressure and viscosity coefficient of each grid point in the flow field; the control parameters include: the time step, the CFL number and the number of inner iteration steps, as well as the incoming flow Mach number, Reynolds number and incoming flow angle of attack; the flow field parameters in step S2 are the lift coefficient and the torque coefficient; the results of solving the motion equation in step S4 include the displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object.

[0016] The Navier-Stokes equation in step S2 is:

[0017]

[0018] In formula (1), Ω is the control volume; S is the boundary surface of the control volume unit; is the outer normal area vector of the S element; Re is the Reynolds number; W is the macroscopic conservation quantity; F is the convective flux; F v is the viscous flux;

[0019] The macroscopic conservation quantity W, convective flux F, and viscous flux F in formula (1) are v The specific expressions are as follows:

[0020]

[0021]

[0022]

[0023] Viscous flux F v The expressions in are as follows:

[0024]

[0025]

[0026]

[0027]

[0028]

[0029]

[0030]

[0031]

[0032] In equations (2) to (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; and μ is the dynamic viscosity coefficient. Based on the eddy viscosity assumption, the viscosity coefficient μ = μ l +μ t , μ l and μ t They are the laminar viscosity coefficient and turbulent viscosity coefficient respectively. The laminar viscosity coefficient can be obtained from the Sutherland formula in is the dimensionless incoming flow temperature, the turbulent viscosity coefficient μ t Given by the turbulence model; λ is given by the Stokes assumption Pr l is the laminar Prandtl number; Pr t is the turbulent Prandtl number, Pr l =0.72, Pr t =0.9; γ is the specific heat ratio, for air, γ = 1.4; M ∞ is the incoming flow Mach number.

[0033] Furthermore, by expressing the convective flux F in the Navier-Stokes equation as the inviscid flux F, Equation (4) can be rewritten as Equation (13):

[0034]

[0035] In formula (13), F1 is the inviscid flux at steady state, which is solved by the lattice Boltzmann equation with BGK approximation; F2 is the flux change caused by grid motion, which is solved directly by taking the average value of the grid physical quantities on both sides of the grid boundary as the physical quantity at the boundary.

[0036] First, let's explain the solution of F2. The physical quantity ρ at the grid boundary i-1 / 2 on the left side of the grid cell (i, j) is i-1 / 2,j ,u i-1 / 2,j , v i-1 / 2,j , T i-1 / 2,j , p i-1 / 2,j , E i-1 / 2,j It can be obtained by averaging the physical quantities of the grids on both sides of the boundary:

[0037]

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] The mesh movement speed u at the boundary b,i-1 / 2,j ,v b,i-1 / 2,j The average velocity of the grid points at both ends of the boundary is obtained:

[0044]

[0045]

[0046] So the value of F2 in formula (13) is:

[0047]

[0048] The lattice Boltzmann equation of the BGK approximation is:

[0049]

[0050] In formula (14), r represents the physical position; τ represents the distribution time for the particle distribution to tend to the equilibrium state through collision; f α is the density distribution function along the α direction; is the corresponding equilibrium state; δ t is the flow time step; e α is the velocity of the particle in the α direction; N is the velocity number of discrete particles;

[0051] Through Chapman-Enskog multiscale expansion [1] It can be proved that: Equation (14) can successfully restore the macroscopic NS equation, establish the macroscopic conservation quantity W and density distribution function fα The relationship between the microscopic physical quantities and the macroscopic physical quantities; using the LB model to solve equation (14), we can get the density distribution function f α , then according to the macroscopic conservation quantity W and density distribution function f α The connection between them gives the macroscopic conservation quantity W, and thus the inviscid flux F1 is obtained.

[0052] Since most existing one-dimensional LB models contain a large number of user-specified parameters, which have a significant impact on the performance of LBFS. To eliminate these deficiencies, Yang et al. [1，2] The non-free parameter D1Q4 model is proposed to simulate inviscid compressible flow with strong shock waves. The density distribution function f is solved by the non-free parameter D1Q4 model. α In the non-free parameter D1Q4 model, the density distribution function f α Discrete into four directions for solution, the results are the density distribution functions g1, g2, g3, g4 in four directions, and the non-free parameter D1Q4 model is written as:

[0053]

[0054]

[0055]

[0056]

[0057]

[0058]

[0059] In formulas (15) to (20), d1 and d2 are the grid velocities; u is the average flow velocity in one dimension; c is defined as The special velocity of the particle, D represents the spatial dimension of the lattice Boltzmann model. For the non-free parameter D1Q4 model, D = 1;

[0060] After obtaining the density distribution functions g1, g2, g3, g4 and the velocities d1 and d2 in the four directions, the following relationship between macroscopic and microscopic physical quantities can be obtained according to the Chapman-Enskog multiscale expansion:

[0061]

[0062]

[0063]

[0064]

[0065]

[0066] In formulas (21) to (25), ξ i is the particle velocity in direction i, ξ1=d1, ξ2=-d1, ξ3=d2, ξ4=-d2; e p is the potential energy of the particle,

[0067] After obtaining the macroscopic physical quantities, the inviscid flux can be solved. The inviscid flux of the unit interface is:

[0068]

[0069] In formula (26), is the inviscid flux at the unit interface; is the inviscid flux calculated without considering numerical dissipation; is the inviscid flux obtained by introducing numerical dissipation calculation; α * is a switching function. In areas with small numerical dissipation such as boundary layers, α * The value of tends to zero, ensuring the introduction of smaller numerical dissipation. In the area where the numerical dissipation is larger, * tends to 1, which can accurately capture strong shock waves, α * The specific way to obtain the value of is as follows:

[0070]

[0071]

[0072]

[0073] α * =max{α L ,α R} (30)

[0074] In formulas (27) to (30), tanh(x) is the hyperbolic tangent function; p L and p R is the pressure on the left and right sides of the unit interface; C is the amplification coefficient; From formula (26), we can see that the range of α is between 0 and 1; α L and α R are the maximum values of the left and right control volume switching functions respectively; N L and N R are the number of faces of the left and right control volumes respectively.

[0075] Inviscid flux calculated without considering numerical dissipation and the inviscid flux obtained by introducing numerical dissipation calculation The solution formulas are:

[0076]

[0077]

[0078] In formula (31)-formula (32), u=n x U n -n y U τ , v = n x U τ +n y U n , U n is the normal velocity at the element interface, U τ is the tangential velocity at the cell interface, n x , n y are the components of the normal vector at the element interface in the x and y directions.

[0079] In the above method, since the density distribution function f α It consists of a balanced part and an unbalanced part. The balanced part contributes to the inviscid flux, and the unbalanced part contributes to the viscous flux. Since LBM is only used to calculate the inviscid flux, the unbalanced part can be regarded as numerical dissipation. When calculating the inviscid flux with small numerical dissipation, the unbalanced part is not considered, and it can provide very accurate results for the boundary layer flow with very little numerical dissipation. However, when simulating hypersonic flows with strong shock waves, this format often exhibits oscillations or even divergence. In order to accurately capture strong shock waves, numerical dissipation needs to be introduced, but large numerical dissipation will affect the solution in the smooth area, resulting in inaccurate solutions in the boundary layer. Therefore, in order to accurately capture strong shock waves and thin boundary layers, we need to carefully control the numerical dissipation. In areas with small numerical dissipation, such as near the boundary layer, numerical dissipation is ignored when calculating the inviscid flux. At this time, the inviscid flux It can be solved by the following formula (31); in the area far away from the boundary layer and prone to strong shock waves, where the numerical dissipation is large, it is necessary to introduce numerical dissipation when calculating the inviscid flux. At this time, the inviscid flux It can be solved by the following formula (32); in order to combine the advantages of the two solutions, the switch function α is introduced * control and The inviscid flux is solved by the proportion of

[0080] The motion equation in step S3 refers to the motion equation of the elastic system, specifically:

[0081]

[0082] In formula (33), h and α are the lifting displacement and pitch angle respectively; m is the unit of mass; S α is the static moment about the elastic axis; I α is the moment of inertia; K h and K α are the lift and pitch spring constants, respectively; L is the lift force; M is the moment about the elastic axis;

[0083] The motion equation of formula (33) takes the half-chord length b as the length dimension and the natural frequency ω of the uncoupled pitching motion as α Non-dimensionalizing the time dimension yields the following results:

[0084]

[0085] In formula (34), x α is the static imbalance; ω h is the natural frequency of the uncoupled lifting motion; is the square of the radius of rotation; U * is the dimensionless velocity, given by Definition; C l and C m are the lift coefficient and moment coefficient respectively;

[0086] Due to the difference in the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the dimensionless time dimension of the structure Need to be readjusted to ensure the consistency of the calculation time scale of the entire system, in is the dimensionless time of the flow field, L is the length scale;

[0087] The lift coefficient C obtained by solving the two-dimensional unsteady flow field Navier-Stokes equation in step S2 is l and moment coefficient C m Update to the dimensionless motion equation of formula (34), decouple the dimensionless motion equation (34), and obtain the following form:

[0088]

[0089] Among them, a1, a2, b1, b2, c1 and c2 are the coefficients of the decoupled equations, all of which are constants;

[0090] make Formula (35) is rewritten as:

[0091]

[0092] By performing second-order discretization on Equation (36), we can obtain:

[0093]

[0094] Arranging the equation (37) into a matrix form yields:

[0095]

[0096] The solution to the equation of motion can be obtained by directly performing matrix inversion.

[0097] The results of the motion equation in step S4 include the values of displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object. Because the dimensionless reference quantities used by the Navier-Stokes equations of the two-dimensional unsteady flow field and the motion equation are different, they need to be transformed. The corresponding relationship is as follows:

[0098]

[0099] Simplifying formula (39) yields formula (40):

[0100]

[0101] In formula (39)-formula (40), h, and are the vertical displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object respectively; h physical 、h f and h s They are the vertical displacement of the object, the dimensionless displacement in the flow field equation, and the dimensionless displacement in the motion equation; and are the velocity of the object in the vertical direction, the dimensionless velocity in the flow field equation, and the dimensionless velocity in the motion equation; and They are the acceleration of the object in the vertical direction, the dimensionless acceleration in the flow field equation, and the dimensionless acceleration in the motion equation; and They are the object's rotational angular velocity in the vertical direction, the dimensionless rotational angular velocity in the flow field equation, and the dimensionless rotational angular velocity in the motion equation; and They are the angular acceleration of the object in the vertical direction, the dimensionless angular acceleration in the flow field equation, and the dimensionless angular acceleration in the motion equation.

[0102] The flow field boundary conditions in step S4 include two types: ideal fluid and viscous fluid:

[0103] For an ideal fluid, the surface adopts a non-penetration boundary condition, that is, the normal velocity of the fluid on the surface is equal to the normal velocity of the surface, such as:

[0104]

[0105] For viscous fluids, the surface adopts the no-slip boundary condition, that is, the velocity of the fluid on the surface is equal to the velocity of the surface, such as:

[0106]

[0107] In formula (41)-formula (42), x t and y t is the object's speed; is the acceleration of the object, which is determined by the laws of motion of the object; is the normal vector to the object surface.

[0108] The grid coordinates in step S4 are updated using the infinite interpolation method (TFI), and the specific method is as follows:

[0109] S41. First, the grid points in the flow field are parameterized based on the arc length. Taking an edge in the i direction as an example, the arc length calculation formula is written as: S i =S i-1 +|r i+1,j -r i,j |,i=1,2,...,imax, then S k / S kmax , k=1,2,…,kmax, which is the parameterized value of each grid point on this edge. The parameterized variables in the i and j directions are respectively denoted as SI i,j ,SJ i,j ;

[0110] S42, using one-dimensional TFI technology to calculate the deformation of each block vertex and edge, when the mesh deformation dP at the edge vertex i=1 and i=imax is known 1,j and dP imax,j After that, the mesh deformation of any point on the edge is determined by formula (43):

[0111] dP i,j =(1-SI i,j )dP 1,j +SI i,j dP imax,j (43)

[0112] The mesh deformation on other edges is calculated according to step S42;

[0113] S43. Use the two-dimensional TFI technique to interpolate the mesh deformation of the internal points. For a mesh surface with imax mesh points in the i direction and jmax mesh points in the j direction, after the mesh deformation on the four edges of the surface is obtained, the deformation of any mesh point on the surface can be written as:

[0114]

[0115] The bending function A i,j 、B i,j 、C i,j and D i,j The definition is as follows:

[0116]

[0117] In formulas (44) and (45), dP 1,1 is the displacement at the grid block corner (1, 1), dP 1,jmax is the displacement at the grid block corner point (1, jmax), dP imax,1 is the displacement at the grid block corner (imax,1), dP imax,jmax is the displacement at the grid block corner (imax,jmax), dP i,1 is the displacement of the grid point (i, j) at the grid line vertex (i, 1) in the j direction, dP i,jmax is the displacement of the grid point (i, j) at the grid line vertex (i, jmax) in the j direction, dP 1,j is the displacement of the grid point (i, j) at the grid line vertex (1, j) in the i direction, dP imax,j is the displacement of the grid point (i, j) at the grid line vertex (imax, j) in the i direction; SI i,1 is the parameterized value of the grid line vertex (i,1) in the j direction where the grid point (i,j) is located on the edge with points (1,1) and (imax,1) as vertices, SI i,jmax SJ is the parameterized value of the grid line vertex (i, jmax) in the j direction where the grid point (i, j) is located on the edge with points (1, jmax) and (imax, jmax) as vertices. 1,j is the parameterized value of the grid line vertex (1, j) in the i direction where the grid point (i, j) is located on the edge with points (1, 1) and (1, jmax) as vertices, SJ imax,j ξ is the parameterized value of the grid line vertex (imax,j) in the i-direction where the grid point (i,j) is located on the edge with points (imax,1) and (imax,jmax) as vertices. i,j ,η i,j is a constant, an intermediate value in the calculation, and has no physical meaning;

[0118] S44. Load the grid point displacement onto the initial grid to generate a new flow field grid: P new =dP i,j +P original .

[0119] In the above-mentioned fluid-solid coupling numerical simulation method, the fluid equation in step S2 includes not only the solution of the convection term, but also the solution of the viscous term. The present invention uses the SA turbulence model to calculate the turbulent viscosity coefficient. The SA turbulence model is a one-equation model proposed by Spalart and Allmaras. The model directly assumes that the turbulent viscosity coefficient satisfies the scalar equation in the flow field (the convection-diffusion equation with the source term), and gives the equation and coefficient based on experience; the model contains a large number of empirical constants and empirical functions. Since Spalart et al. have rich experience in aviation calculations and have mastered rich experimental (calculation) data, these empirical constants and empirical functions have been well adjusted. The SA turbulence model has a good calculation effect in aviation, so the present invention directly uses the SA turbulence model to solve the viscous term.

[0120] The turbulent viscosity coefficient μ given by the SA turbulence model t Use the following formula to calculate:

[0121]

[0122]

[0123]

[0124]

[0125]

[0126] in right Limit it to prevent it from taking negative values, which would cause calculation instability; ||Ω|| is the modulus of the vorticity, d is the distance from the grid point to the wall; C b1 =0.1355,σ=2 / 3,C b2 =0.622, k=0.41, C w2 =0.3, C w3 =2, C v1 =7.1, C t3 =1.2, C t4 =0.5; the boundary conditions are set as follows, wall conditions: Flow conditions: The exit condition is obtained by extrapolation, that is, the value of the virtual grid point outside the exit is equal to the value of the inner point.

[0127] The Lattice Boltzmann Method (LBM) is an emerging mesoscopic numerical method. Its principle considers the migration and collision processes of particles near interfaces at the mesoscopic level and reconstructs macroscopic physical quantities at the interface using the obtained particle velocity distribution function. This approach offers richer and clearer physical implications than traditional numerical methods, and has advantages such as suitability for parallel computing, simple boundary handling, and ease of programming.

[0128] The present invention has the following advantages over the prior art:

[0129] The fluid-solid coupling numerical simulation method of the present invention introduces the lattice Boltzmann method to solve the inviscid flux of the flow field. The lattice Boltzmann method is a numerical simulation method based on contemporary statistical physics. It simulates complex physical phenomena at the macroscopic level by tracking the interaction of a large number of discrete particles at the microscopic scale of the medium. The simulation method has good stability, high calculation accuracy, clear physical meaning, and does not need to consider nonlinear terms. It can improve the capture of detailed characteristics of the flow field from a physical level, thereby effectively improving the aerodynamic simulation capability.

[0130] The innovation of the fluid-solid coupling numerical simulation method of the present invention lies in applying the lattice Boltzmann flux method to solve the unsteady flow field, and verifying the accuracy of the lattice Boltzmann flux method in solving the unsteady flow field. BRIEF DESCRIPTION OF THE DRAWINGS

[0131] Attachment Figure 1 Flowchart of the fluid-solid coupling numerical simulation method based on the lattice Boltzmann flux algorithm of the present invention;

[0132] Attachment Figure 2 Schematic diagram of the elastically mounted airfoil in the present invention;

[0133] Attachment Figure 3 Shows M = 0.825, V * = 0.5, the vertical displacement and rotation angle of the airfoil when the damping response occurs are at any time

[0134] ∞

[0135] the process of change between

[0136] Attachment Figure 4 Shows M = 0.825, V * = 0.51, the vertical displacement and rotation angle of the airfoil when the neutral response occurs are at any time

[0137] ∞

[0138] the process of change between

[0139] Attachment Figure 5Shows M = 0.825, V * = 0.65, the vertical displacement and rotation angle of the airfoil when the divergence response occurs are

[0140] ∞

[0141] the changing process of time;

[0142] Attachment Figure 6 The speed index V at which neutral response occurs at different Mach numbers is shown. * .

[0143] Attachment Figure 7 The flutter frequency ratio ω / ω for neutral response is shown at different Mach numbers. α . DETAILED DESCRIPTION

[0144] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0145] like Figure 1 As shown: A fluid-structure interaction numerical simulation method based on the lattice Boltzmann flux algorithm. The steps of this method are as follows:

[0146] S1, determine the calculation area, read the grid information, physical parameters and control parameters, and enter step S2;

[0147] S2. Solve the Navier-Stokes equations for the two-dimensional unsteady flow field, where the convection term is solved using the lattice Boltzmann method and the viscous term is solved using the SA turbulence model, then update the flow field parameters and proceed to step S3;

[0148] S3, solve the motion equation and go to step S4;

[0149] S4. Update the flow field boundary conditions and grid coordinates based on the results of solving the motion equations in step S3 (the displacement caused by the object movement requires updating the grid coordinates, and the surface velocity and acceleration caused by the object movement will affect the boundary conditions of the flow field calculation), and then proceed to step S5;

[0150] S5, determine whether the simulation time has been reached, if yes, proceed to step S7, if not, proceed to step S6;

[0151] S6, return to step S2 to perform the next time step calculation;

[0152] S7. Simulation ends.

[0153] The grid information in the above step S1 includes: the number of grid blocks, the number of grids in each grid block, the connection information of the grid blocks and the grid coordinates; the physical parameters include: the density, velocity, temperature, pressure and viscosity coefficient of each grid point in the flow field; the control parameters include: the time step, the CFL number and the number of inner iteration steps, as well as the incoming flow Mach number, Reynolds number and incoming flow angle of attack; the flow field parameters in the above step S2 are the lift coefficient and the torque coefficient; the results of solving the motion equation in the above step S4 include the displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object.

[0154] The specific steps for solving the Navier-Stokes equations for the two-dimensional unsteady flow field in step S2 are as follows.

[0155] The Navier-Stokes equations are:

[0156]

[0157] In formula (1), Ω is the control volume; S is the boundary surface of the control volume unit; is the outer normal area vector of the S element; Re is the Reynolds number; W is the macroscopic conservation quantity; F is the convective flux; F v is the viscous flux;

[0158] The macroscopic conservation quantity W, convective flux F, and viscous flux F in formula (1) are v The specific expressions are as follows:

[0159]

[0160]

[0161]

[0162] Viscous flux F v The expressions in are as follows:

[0163]

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170]

[0171] In equations (2) to (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; and μ is the dynamic viscosity coefficient. Based on the eddy viscosity assumption, the viscosity coefficient μ = μ l +μ t , μ l and μ t They are the laminar viscosity coefficient and turbulent viscosity coefficient respectively. The laminar viscosity coefficient can be obtained from the Sutherland formula in is the dimensionless incoming flow temperature, the turbulent viscosity coefficient μ t Given by the turbulence model; λ is given by the Stokes assumption Pr l is the laminar Prandtl number; Pr t is the turbulent Prandtl number, Pr l =0.72, Pr t =0.9; γ is the specific heat ratio, for air, γ = 1.4; M ∞ is the incoming flow Mach number.

[0172] By expressing the convective flux F in the Navier-Stokes equation as the inviscid flux F, Equation (4) can be rewritten as Equation (13):

[0173]

[0174] In formula (13), F1 is the inviscid flux at steady state, which is solved by the lattice Boltzmann equation with BGK approximation; F2 is the flux change caused by grid motion, which is solved directly by taking the average value of the grid physical quantities on both sides of the grid boundary as the physical quantity at the boundary.

[0175] The following is a detailed explanation of the solution to the steady-state inviscid flux F1. The lattice Boltzmann equation of the BGK approximation is:

[0176]

[0177] In formula (14), r represents the physical position; τ represents the distribution time for the particle distribution to tend to the equilibrium state through collision; f α is the density distribution function along the α direction; is the corresponding equilibrium state; δ t is the flow time step; e α is the velocity of the particle in the α direction; N is the velocity number of discrete particles;

[0178] Through Chapman-Enskog multiscale expansion [1]It can be proved that: Equation (14) can successfully restore the macroscopic NS equation, establish the macroscopic conservation quantity W and density distribution function f α The relationship between the microscopic physical quantities and the macroscopic physical quantities; using the LB model to solve equation (14), we can get the density distribution function f α , then according to the macroscopic conservation quantity W and density distribution function f α The connection between them gives the macroscopic conservation quantity W, and thus the inviscid flux F1 is obtained.

[0179] The following non-free parameter D1Q4 model is used to solve the density distribution function f α In the non-free parameter D1Q4 model, the density distribution function f α Discrete into four directions for solution, the results are the density distribution functions g1, g2, g3, g4 in four directions, and the non-free parameter D1Q4 model is written as:

[0180]

[0181]

[0182]

[0183]

[0184]

[0185]

[0186] In formulas (15) to (20), d1 and d2 are the grid velocities; u is the average flow velocity in one dimension; c is defined as The special velocity of the particle, D represents the spatial dimension of the lattice Boltzmann model. For the non-free parameter D1Q4 model, D = 1;

[0187] After obtaining the density distribution functions g1, g2, g3, g4 and the velocities d1 and d2 in the four directions, the following relationship between macroscopic and microscopic physical quantities can be obtained according to the Chapman-Enskog multiscale expansion:

[0188]

[0189]

[0190]

[0191]

[0192]

[0193] In formulas (21) to (25), ξ i is the particle velocity in direction i, ξ1=d1, ξ2=-d1, ξ3=d2, ξ4=-d2; e p is the potential energy of the particle,

[0194] After obtaining the macroscopic physical quantities, the inviscid flux can be solved. The inviscid flux of the unit interface is:

[0195]

[0196] In formula (26), is the inviscid flux at the unit interface; is the inviscid flux calculated without considering numerical dissipation; is the inviscid flux obtained by introducing numerical dissipation calculation; α * is a switching function. In areas with small numerical dissipation such as boundary layers, α * The value of tends to zero, ensuring the introduction of smaller numerical dissipation. In the area where the numerical dissipation is larger, * tends to 1, which can accurately capture strong shock waves, α * The specific way to obtain the value of is as follows:

[0197]

[0198]

[0199]

[0200] α * =max{α L ,α R} (30)

[0201] In formulas (27) to (30), tanh(x) is the hyperbolic tangent function; p L and p R is the pressure on the left and right sides of the unit interface; C is the amplification factor, and C = 100 is used in this paper; From formula (26), we can see that the range of α is between 0 and 1; α L and α R are the maximum values of the left and right control volume switching functions respectively; N L and N R are the number of faces of the left and right control volumes respectively.

[0202] Furthermore, the inviscid flux calculated without considering numerical dissipation is and the inviscid flux obtained by introducing numerical dissipation calculation The solution formulas are:

[0203]

[0204]

[0205] In formula (31)-formula (32), u=n x U n -n y U τ , v = n x U τ +n y U n , U n is the normal velocity at the element interface, U τ is the tangential velocity at the cell interface, n x , n y are the components of the normal vector at the element interface in the x and y directions.

[0206] The lift coefficient and torque coefficient obtained by solving the flow field equation are transferred to the solution part of the object's motion equation.

[0207] Next, we proceed to solve the motion equation. The motion equation in step S3 refers to the motion equation of the elastic system, specifically:

[0208]

[0209] In formula (33), h and α are the lifting displacement and pitch angle respectively; m is the unit of mass; S α is the static moment about the elastic axis; I α is the moment of inertia; K h and K α are the lift and pitch spring constants, respectively; L is the lift force; M is the moment about the elastic axis;

[0210] The motion equation of formula (33) takes the half-chord length b as the length dimension and the natural frequency ω of the uncoupled pitching motion as α Non-dimensionalizing the time dimension yields the following results:

[0211]

[0212] In formula (34), x α is the static imbalance; ω h is the natural frequency of the uncoupled lifting motion; is the square of the radius of rotation; U * is the dimensionless velocity, given by Definition; C l and C m are the lift coefficient and moment coefficient respectively;

[0213] Due to the difference in the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the dimensionless time dimension of the structure Need to be readjusted to ensure the consistency of the calculation time scale of the entire system, in is the dimensionless time of the flow field, L is the length scale;

[0214] The lift coefficient C obtained by solving the two-dimensional unsteady flow field Navier-Stokes equation in step S2 is l and moment coefficient C m Update to the dimensionless motion equation of formula (34), decouple the dimensionless motion equation (34), and obtain the following form:

[0215]

[0216] Among them, a1, a2, b1, b2, c1 and c2 are the coefficients of the decoupled equations, all of which are constants;

[0217] make Formula (35) is rewritten as:

[0218]

[0219] By performing second-order discretization on Equation (36), we can obtain:

[0220]

[0221] Arranging the equation (37) into a matrix form yields:

[0222]

[0223] The solution to the equation of motion can be obtained by directly performing a matrix inversion. After solving the equation of motion, the motion parameters such as the displacement, velocity, acceleration, angular velocity, and angular acceleration of the object are updated. Based on the updated motion parameters, the flow field boundary conditions and grid coordinates are updated, and the solution for the next time step is started or ended.

[0224] Because the dimensionless reference quantities used in the Navier-Stokes equations and the equation of motion for two-dimensional unsteady flow are different, a transformation is required. The corresponding relationship is as follows:

[0225]

[0226] Simplifying formula (39) yields formula (40):

[0227]

[0228] In formula (39)-formula (40), h, and are the vertical displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object respectively; h physical 、h f and h s They are the vertical displacement of the object, the dimensionless displacement in the flow field equation, and the dimensionless displacement in the motion equation; and are the velocity of the object in the vertical direction, the dimensionless velocity in the flow field equation, and the dimensionless velocity in the motion equation; and They are the acceleration of the object in the vertical direction, the dimensionless acceleration in the flow field equation, and the dimensionless acceleration in the motion equation; and They are the object's rotational angular velocity in the vertical direction, the dimensionless rotational angular velocity in the flow field equation, and the dimensionless rotational angular velocity in the motion equation; and They are the angular acceleration of the object in the vertical direction, the dimensionless angular acceleration in the flow field equation, and the dimensionless angular acceleration in the motion equation.

[0229] In the above method, the flow field boundary conditions in step S4 include two types: ideal fluid and viscous fluid. For the ideal fluid, the object surface adopts the non-penetration boundary condition, that is, the normal velocity of the fluid on the object surface is equal to the normal velocity of the object surface, such as:

[0230]

[0231] For viscous fluids, the surface adopts the no-slip boundary condition, that is, the velocity of the fluid on the surface is equal to the velocity of the surface, such as:

[0232]

[0233] In formula (41)-formula (42), x t and y t is the object's speed; is the acceleration of the object, which is determined by the laws of motion of the object; is the normal vector to the object surface.

[0234] In the above method, the grid coordinates in step S4 are updated using the infinite interpolation method (TFI), and the specific method is as follows:

[0235] S41. First, the grid points in the flow field are parameterized based on the arc length. Taking an edge in the i direction as an example, the arc length calculation formula is written as: S i =S i-1 +|r i+1,j -r i,j |,i=1,2,...,imax, then Sk / S kmax , k=1,2,…,kmax, which is the parameterized value of each grid point on this edge. The parameterized variables in the i and j directions are respectively denoted as SI i,j ,SJ i,j ;

[0236] S42, using one-dimensional TFI technology to calculate the deformation of each block vertex and edge, when the mesh deformation dP at the edge vertex i=1 and i=imax is known 1,j and dP imax,j After that, the mesh deformation of any point on the edge is determined by formula (43):

[0237] dP i,j =(1-SI i,j )dP 1,j +SI i,j dP imax,j (43)

[0238] The mesh deformation on other edges is calculated according to step S42;

[0239] S43. Use the two-dimensional TFI technique to interpolate the mesh deformation of the internal points. For a mesh surface with imax mesh points in the i direction and jmax mesh points in the j direction, after the mesh deformation on the four edges of the surface is obtained, the deformation of any mesh point on the surface can be written as:

[0240]

[0241] The bending function A i,j 、B i,j 、C i,j and D i,j The definition is as follows:

[0242]

[0243] In formulas (44) and (45), dP 1,1 is the displacement at the grid block corner (1, 1), dP 1,jmax is the displacement at the grid block corner point (1, jmax), dP imax,1 is the displacement at the grid block corner (imax,1), dP imax,jmax is the displacement at the grid block corner (imax,jmax), dP i,1 is the displacement of the grid point (i, j) at the grid line vertex (i, 1) in the j direction, dP i,jmax is the displacement of the grid point (i, j) at the grid line vertex (i, jmax) in the j direction, dP 1,j is the displacement of the grid point (i, j) at the grid line vertex (1, j) in the i direction, dPimax,j is the displacement of the grid point (i, j) at the grid line vertex (imax, j) in the i direction; SI i,1 is the parameterized value of the grid line vertex (i,1) in the j direction where the grid point (i,j) is located on the edge with points (1,1) and (imax,1) as vertices, SI i,jmax SJ is the parameterized value of the grid line vertex (i, jmax) in the j direction where the grid point (i, j) is located on the edge with points (1, jmax) and (imax, jmax) as vertices. 1,j is the parameterized value of the grid line vertex (1, j) in the i direction where the grid point (i, j) is located on the edge with points (1, 1) and (1, jmax) as vertices, SJ imax,j ξ is the parameterized value of the grid line vertex (imax,j) in the i-direction where the grid point (i,j) is located on the edge with points (imax,1) and (imax,jmax) as vertices. i,j ,η i,j is a constant, an intermediate value in the calculation, and has no physical meaning;

[0244] S44. Load the grid point displacement onto the initial grid to generate a new flow field grid: P new =dP i,j +P original .

[0245] Verification Example

[0246] In order to verify the correctness of the algorithm of the present invention, the following example is carried out, and the model information used is as follows: Figure 2 As shown in Figure 1, a structural model of the vibration of a two-dimensional swept wing with a NACA64A010 cross section is used. The airfoil undergoes pitch and vertical movement about a given elastic axis. The pitch axis is defined by the distance a, which is the distance between the elastic axis and the midchord as a percentage of the half-chord length. A positive a indicates that the axis is downstream of the midchord, while a negative a indicates that it is upstream of the midchord. The structural parameters used in this model are shown in Table 1.

[0247] Table 1

[0248]

[0249] Due to the symmetrical shape of the NACA 64A010 airfoil, an initial perturbation is applied to trigger the oscillatory motion. The airfoil moves sinusoidally at the natural pitch frequency ω α The airfoil is pitched about its elastic axis with an amplitude of α0 = 1°. The forced pitch mode typically lasts for 1-3 cycles. The elastically mounted airfoil is then set to free motion in both the vertical and pitch directions, and the dynamic response is recorded.

[0250] Figure 3-5 Showing M∞ =0.825, three different V * The vertical displacement and rotation angle of the airfoil vary with time under the condition of . These figures correspond to the damping response, neutral response and divergent response respectively. The main task of calculating the flutter boundary is to determine the critical position when the neutral response occurs by observing these images. When V * When the value of is less than the flutter boundary, a damping response occurs, and the changes in vertical displacement and pitch angle decay over time, such as Figure 3 As shown; when V * When the value of is close to the flutter boundary, a neutral response occurs, and the vertical displacement and pitch angle basically do not change with time, such as Figure 4 As shown; when V * When the value of is greater than the flutter boundary, a divergent response occurs, and the changes in vertical displacement and pitch angle increase with time, such as Figure 5 shown.

[0251] Then calculate the speed index V at which neutral response occurs at different Mach numbers. * and the flutter frequency ratio ω / ω α , the results are as follows Figure 6 and Figure 7 shown. Figure 6 and Figure 7 The calculated flutter velocity boundary and flutter frequency ratio boundary are shown, which are in good agreement with the previous research results. The flutter boundary is a "double step" shape. ∞ =0.83 and M ∞ = 0.91, the flutter velocity jumps twice. The flutter boundary can be roughly divided into four regions: Subsonic region: M ∞ <0.80; Transonic pit: M ∞ =0.80~0.83; the first vibration speed jump area: M ∞ =0.84~0.91; Locking area: M ∞ >0.91. The results of this paper differ from those of previous studies, mainly due to the different turbulence models and Reynolds numbers used in the simulation. In this study, all calculations were performed with a fixed Re=1.256×10 6 .

[0252] The innovation of the present invention lies in applying the lattice Boltzmann flux method to solve unsteady flow fields. The above examples prove the accuracy of the lattice Boltzmann flux method in solving unsteady flow fields.

[0253] The fluid-structure interaction numerical simulation method of this invention incorporates the lattice Boltzmann method to solve the inviscid flux of the flow field. The lattice Boltzmann method is a numerical simulation method based on contemporary statistical physics. By tracking the interactions of a large number of discrete particles in the medium at the microscopic scale, it simulates complex physical phenomena at the macroscopic level. This simulation method has good stability, high computational accuracy, clear physical meaning, and does not consider nonlinear terms. It can improve the capture of detailed flow field characteristics at the physical level, thereby effectively enhancing the ability to simulate aerodynamic forces.

[0254] The above embodiments are only for illustrating the technical ideas of the present invention and cannot be used to limit the scope of protection of the present invention. Any changes made on the basis of the technical solutions in accordance with the technical ideas proposed by the present invention fall within the scope of protection of the present invention; any technologies not involved in the present invention can be implemented by existing technologies.

[0255] References:

[0256] [1]He

[0257] [2]LMYang,C.Shu,J.Wu.Development and comparative studies of threenon-free parameter lattice Boltzmann.

Claims

1. A fluid-structure interaction numerical simulation method based on the lattice Boltzmann flux algorithm, characterized by: The steps of this method are as follows: S1, determine the calculation area, read the grid information, physical parameters and control parameters, and enter step S2; S2, solving the Navier-Stokes equations for the two-dimensional unsteady flow field, wherein the convection term is solved using the lattice Boltzmann method, and then the flow field parameters are updated, and then proceeding to step S3; S3, solving the motion equation of the elastic system and proceeding to step S4; S4, updating the flow field boundary conditions and grid coordinates according to the results of solving the motion equation in step S3, and proceeding to step S5; S5, determine whether the simulation time has been reached, if yes, proceed to step S7, if not, proceed to step S6; S6, return to step S2 to perform the next time step calculation; S7, simulation ends; The Navier-Stokes equation in step S2 is: In formula (1), Ω is the control volume; S is the boundary surface of the control volume element; is the outer normal area vector of the S element; Re is the Reynolds number; W is the macroscopic conservation quantity; F is the convective flux; F v is the viscous flux; The macroscopic conservation quantity W, convective flux F, and viscous flux F in formula (1) are v The specific expressions are as follows: Viscous flux F v The expressions in are as follows: In equations (2) to (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; and μ is the dynamic viscosity coefficient. Based on the eddy viscosity assumption, the viscosity coefficient μ = μ l +μ t , μ l and μ t They are the laminar viscosity coefficient and turbulent viscosity coefficient respectively. The laminar viscosity coefficient can be obtained from the Sutherland formula in is the dimensionless incoming flow temperature, the turbulent viscosity coefficient μ t Given by the turbulence model; λ is given by the Stokes assumption Pr l is the laminar Prandtl number; Pr t is the turbulent Prandtl number, Pr l =0.72, Pr t =0.9; γ is the specific heat ratio, for air, γ = 1.4; M ∞ is the incoming flow Mach number; By expressing the convective flux F in the Navier-Stokes equation as the inviscid flux F, Equation (4) can be rewritten as Equation (13): In formula (13), F1 is the inviscid flux at steady state, which is solved by the lattice Boltzmann equation with BGK approximation; F2 is the flux change caused by grid motion, which is solved directly by taking the average value of the grid physical quantities on both sides of the grid boundary as the physical quantity at the boundary; The lattice Boltzmann equation of the BGK approximation is: In formula (14), r represents the physical position; τ represents the distribution time for the particle distribution to tend to the equilibrium state through collision; f α is the density distribution function along the α direction; is the corresponding equilibrium state; δ t is the flow time step; e α is the velocity of the particle in the α direction; N is the velocity number of discrete particles; Through the Chapman-Enskog multi-scale expansion, it can be proved that Equation (14) can be successfully restored to the macroscopic NS equation, and the macroscopic conservation quantity W and density distribution function f can be established. α The relationship between the microscopic physical quantities and the macroscopic physical quantities; using the LB model to solve equation (14), we can get the density distribution function f α , then according to the macroscopic conservation quantity W and density distribution function f α The connection between them gives the macroscopic conservation quantity W, and thus the inviscid flux F1 is obtained.

2. The fluid-structure coupling numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 1 is characterized in that: The grid information in step S1 includes: the number of grid blocks, the number of grids in each grid block, the connection information of the grid blocks and the grid coordinates; the physical property parameters include: the density, velocity, temperature, pressure and viscosity coefficient of each grid point in the flow field; the control parameters include: the time step, the CFL number and the number of inner iteration steps, as well as the incoming flow Mach number, Reynolds number and incoming flow angle of attack; the flow field parameters in step S2 are the lift coefficient and the torque coefficient.

3. The fluid-structure coupling numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 1 is characterized in that: The density distribution function f is solved using the non-free parameter D1Q4 model α In the non-free parameter D1Q4 model, the density distribution function f α Discrete into four directions for solution, the results are the density distribution functions g1, g2, g3, g4 in four directions, and the non-free parameter D1Q4 model is written as: In formulas (15) to (20), d1 and d2 are the grid velocities; u is the average flow velocity in one dimension; c is defined as The special velocity of the particle, D represents the spatial dimension of the lattice Boltzmann model. For the non-free parameter D1Q4 model, D = 1; After obtaining the density distribution functions g1, g2, g3, g4 and the velocities d1 and d2 in the four directions, the following relationship between macroscopic and microscopic physical quantities can be obtained according to the Chapman-Enskog multiscale expansion: In formulas (21) to (25), ξ i is the particle velocity in direction i, ξ1=d1, ξ2=-d1, ξ3=d2, ξ4=-d2; e p is the potential energy of the particle, After obtaining the macroscopic physical quantities, the inviscid flux can be solved. The inviscid flux of the unit interface is: In formula (26), is the inviscid flux at the unit interface; is the inviscid flux calculated without considering numerical dissipation; is the inviscid flux obtained by introducing numerical dissipation calculation; α * is the switching function. In the region where the numerical dissipation of the boundary layer is small, α * The value of tends to zero, ensuring the introduction of smaller numerical dissipation. In the area where the numerical dissipation is larger, * tends to 1, which can accurately capture strong shock waves, α * The specific way to obtain the value of is as follows: a * =max{α L ,a R } (30) In formulas (27) to (30), tanh(x) is the hyperbolic tangent function; p L and p R is the pressure on the left and right sides of the unit interface; C is the amplification coefficient; From formula (26), we can see that the range of α is between 0 and 1; α L and α R are the maximum values of the left and right control volume switching functions respectively; N L and N R are the number of faces of the left and right control volumes respectively.

4. The fluid-structure coupling numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 3 is characterized in that: Inviscid flux calculated without considering numerical dissipation and the inviscid flux obtained by introducing numerical dissipation calculation The solution formulas are: In formula (31)-formula (32), u=n x U n -n y U τ , v = n x U τ +n y U n , U n is the normal velocity at the element interface, U τ is the tangential velocity at the cell interface, n x , n y are the components of the normal vector at the element interface in the x and y directions.

5. The fluid-structure interaction numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 1, characterized in that: The motion equation of the elastic system in step S3 is specifically: In formula (33), h and α are the lifting displacement and pitch angle respectively; m is the unit of mass; S α is the static moment about the elastic axis; I α is the moment of inertia; K h and K α are the lift and pitch spring constants, respectively; L is the lift force; M is the moment about the elastic axis; The motion equation of formula (33) takes the half-chord length b as the length dimension and the natural frequency ω of the uncoupled pitching motion as α Non-dimensionalizing the time dimension yields the following results: In formula (34), x α is the static imbalance; ω h is the natural frequency of the uncoupled lifting motion; is the square of the radius of rotation; U * is the dimensionless velocity, given by Definition; C l and C m are the lift coefficient and moment coefficient respectively; Due to the difference in the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the dimensionless time dimension of the structure Need to be readjusted to ensure the consistency of the calculation time scale of the entire system, in is the dimensionless time of the flow field, L is the length scale; The lift coefficient C obtained by solving the two-dimensional unsteady flow field Navier-Stokes equation in step S2 is l and moment coefficient C m Update to the dimensionless motion equation of formula (34), Decoupling the dimensionless motion equation (34) yields the following form: Among them, a1, a2, b1, b2, c1 and c2 are the coefficients of the decoupled equations, all of which are constants; make Formula (35) is rewritten as: By performing second-order discretization on Equation (36), we can obtain: Arranging the equation (37) into a matrix form yields: The solution to the equation of motion can be obtained by directly performing matrix inversion.

6. The fluid-structure interaction numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 1, characterized in that: The results of the motion equation in step S4 include the values of displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object. Because the dimensionless reference quantities used by the Navier-Stokes equations of the two-dimensional unsteady flow field and the motion equation are different, they need to be transformed. The corresponding relationship is as follows: Simplifying formula (39) yields formula (40): In formula (39)-formula (40), h, and are the vertical displacement, velocity, acceleration, rotational angular velocity and angular acceleration of the object respectively; h physical 、h f and h s They are the vertical displacement of the object, the dimensionless displacement in the flow field equation, and the dimensionless displacement in the motion equation; and are the velocity of the object in the vertical direction, the dimensionless velocity in the flow field equation, and the dimensionless velocity in the motion equation; and They are the acceleration of the object in the vertical direction, the dimensionless acceleration in the flow field equation, and the dimensionless acceleration in the motion equation; and They are the object's rotational angular velocity in the vertical direction, the dimensionless rotational angular velocity in the flow field equation, and the dimensionless rotational angular velocity in the motion equation; and They are the angular acceleration of the object in the vertical direction, the dimensionless angular acceleration in the flow field equation, and the dimensionless angular acceleration in the motion equation.

7. The fluid-structure interaction numerical simulation method based on the lattice Boltzmann flux algorithm according to claim 1 or 6, characterized in that: The flow field boundary conditions in step S4 include two types: ideal fluid and viscous fluid: For an ideal fluid, the surface adopts a non-penetration boundary condition, that is, the normal velocity of the fluid on the surface is equal to the normal velocity of the surface, such as: For viscous fluids, the surface adopts the no-slip boundary condition, that is, the velocity of the fluid on the surface is equal to the velocity of the surface, such as: In formula (41)-formula (42), x t and y t is the object's speed; is the acceleration of the object, which is determined by the laws of motion of the object; is the normal vector to the object surface; The grid coordinates in step S4 are updated using the infinite interpolation method (TFI), and the specific method is as follows: S41. First, the grid points in the flow field are parameterized based on the arc length. Taking an edge in the i direction as an example, the arc length calculation formula is written as: S i =S i-1 +|r i+1,j -r i,j |,i=1,2,...,imax, then S k / S kmax , k=1,2,…,kmax, which is the parameterized value of each grid point on this edge. The parameterized variables in the i and j directions are respectively denoted as SI i,j ,SJ i,j ; S42, using one-dimensional TFI technology to calculate the deformation of each block vertex and edge, when the mesh deformation dP at the edge vertex i=1 and i=imax is known 1,j and dP imax,j After that, the mesh deformation of any point on the edge is determined by formula (43): dP i,j =(1-SI i,j )dP 1,j +SI i,j dP imax,j (43) The mesh deformation on other edges is calculated according to step S42; S43. Use the two-dimensional TFI technique to interpolate the mesh deformation of the internal points. For a mesh surface with imax mesh points in the i direction and jmax mesh points in the j direction, after the mesh deformation on the four edges of the surface is obtained, the deformation of any mesh point on the surface can be written as: The bending function A i,j 、B i,j 、C i,j and D i,j The definition is as follows: In formulas (44) and (45), dP 1,1 is the displacement at the grid block corner (1, 1), dP 1,jmax is the displacement at the grid block corner point (1, jmax), dP imax,1 is the displacement at the grid block corner (imax,1), dP imax,jmax is the displacement at the grid block corner (imax,jmax), dP i,1 is the displacement of the grid point (i, j) at the grid line vertex (i, 1) in the j direction, dP i,jmax is the displacement of the grid point (i, j) at the grid line vertex (i, jmax) in the j direction, dP 1,j is the displacement of the grid point (i, j) at the grid line vertex (1, j) in the i direction, dP imax,j is the displacement of the grid point (i, j) at the grid line vertex (imax, j) in the i direction; SI i,1 is the parameterized value of the grid line vertex (i,1) in the j direction where the grid point (i,j) is located on the edge with points (1,1) and (imax,1) as vertices, SI i,jmax SJ is the parameterized value of the grid line vertex (i, jmax) in the j direction where the grid point (i, j) is located on the edge with points (1, jmax) and (imax, jmax) as vertices. 1,j is the parameterized value of the grid line vertex (1, j) in the i direction where the grid point (i, j) is located on the edge with points (1, 1) and (1, jmax) as vertices, SJ imax,j is the parameterized value of the grid line vertex (imax,j) in the i-direction where the grid point (i,j) is located on the edge with points (imax,1) and (imax,jmax) as vertices; i,j ,η i,j is a constant, an intermediate value in the calculation, and has no physical meaning; S44. Load the grid point displacement onto the initial grid to generate a new flow field grid: P new =dP i,j +P original .

Images