CN115828782A - Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm - Google Patents
Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm Download PDFInfo
- Publication number
- CN115828782A CN115828782A CN202211546774.5A CN202211546774A CN115828782A CN 115828782 A CN115828782 A CN 115828782A CN 202211546774 A CN202211546774 A CN 202211546774A CN 115828782 A CN115828782 A CN 115828782A
- Authority
- CN
- China
- Prior art keywords
- grid
- equation
- flux
- formula
- motion
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm, which comprises the following steps: s1, determining a calculation area, and reading grid information, physical property parameters and control parameters; s2, solving a Navier-Stokes equation of the two-dimensional unsteady flow field, wherein a convection term is solved by using a lattice Boltzmann method, and then updating flow field parameters; s3, solving a motion equation; s4, updating the boundary conditions of the flow field and grid coordinates according to the result of solving the motion equation in the step S3; s5, judging whether the simulation time is up, if so, entering a step S7, and if not, entering a step S6; s6, returning to the step S2 to calculate the next time step; and S7, ending the simulation. The method applies the lattice boltzmann flux method to the inviscid flux solution of the unsteady flow field, and can improve the capture of the flow field detail characteristics from the physical level, thereby effectively improving the aerodynamic simulation capability.
Description
Technical Field
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
Fluid-structure interaction (FSI) is an independent branch formed by intersecting Fluid mechanics and solid mechanics, and is mainly used for researching the actions of deformation, dynamic performance and the like of a solid under the action of a Fluid and the influence of motion on the Fluid. The research starts from the aeroelasticity problem of airplane wings, mainly calculates flutter boundaries of the wings at first, and gradually develops into wider fields such as water conservancy, ocean, aerospace, medical treatment and the like.
In the prior art, commonly used flow field numerical simulation methods, such as a Finite Difference Method (FDM), a Finite Volume Method (FVM), and the like, solve a unit boundary flux by dividing the unit boundary flux into two parts, i.e., a non-viscous part and a viscous part, the calculation of the non-viscous part mostly adopts a numerical format constructed based on Euler equation characteristics and often treats a unit interface as a discontinuity (e.g., riemann solver), and the solution of the viscous flux generally adopts a central difference based on continuous assumption of an interface physical quantity. In addition to this contradiction, both the finite difference and finite volume methods have specific numerical deficiencies and most require some manual "correction" due to the "inherent deficiency" of Euler's equations in describing the true physical evolution of the interface.
The flow field solving technology at least has the following problems: when the flow field parameter values are calculated, the flow field parameter values are limited by some continuous assumed conditions of the equation. Due to the deviation between the assumed condition and the actual condition, the flow field obtained by simulation has larger deviation with the actual flow field, so that larger error exists between the obtained flow field parameter value and the actual value, and the application of the flow field parameter value in actual production is limited.
Disclosure of Invention
The invention aims to provide a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm, aiming at the problems in the prior art.
The invention aims to solve the problems by the following technical scheme:
a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm is characterized in that: the method comprises the following steps:
s1, determining a calculation area, reading grid information, physical property parameters and control parameters, and entering a step S2;
s2, solving a Navier-Stokes equation of the two-dimensional unsteady flow field, wherein a convection term is solved by using a lattice Boltzmann method, a viscosity term is solved by using an SA turbulence model, then updating flow field parameters, and entering a step S3;
s3, solving a motion equation, and entering the step S4;
s4, updating the boundary conditions of the flow field and grid coordinates (the grid coordinates need to be updated due to the displacement caused by the motion of the object, and the boundary conditions of the flow field calculation are influenced by the speed and the acceleration of the object plane caused by the motion of the object) according to the result of solving the motion equation in the step S3, and entering a step S5;
s5, judging whether the simulation time is up, if so, entering a step S7, and if not, entering a step S6;
s6, returning to the step S2 to calculate the next time step;
and S7, ending the simulation.
The mesh information in step S1 includes: the number of grid blocks, the grid number of each grid block, the connection information of the grid blocks and grid coordinates, and the physical property parameters comprise: density, speed, temperature, pressure and viscosity coefficient of each grid point of the flow field, and control parameters comprise: the method comprises the following steps of time step length, CFL (computational fluid dynamics) number, inner iteration step number, incoming flow Mach number, reynolds number and incoming flow incidence angle; the flow field parameters in the step S2 are a lift coefficient and a moment coefficient; the result of solving the motion equation in step S4 includes the displacement, velocity, acceleration, angular velocity of rotation, and angular acceleration of the object.
The Navier-Stokes equation in the step S2 is as follows:
in the formula (1), Ω is a controller; s is a control volume element boundary surface;
is the external normal area vector of the S infinitesimal; re is Reynolds number; w is macroscopic conservation quantity; f is convection flux; f v Is the viscous flux;
macroscopic conservation quantity W, convection flux F, and viscous flux F in formula (1) v The specific expressions of (a) are respectively as follows:
viscous flux F v The expressions in (1) are as follows:
in the formula (2) -formula (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; μ is the kinetic viscosity coefficient, based on the vortex viscosity assumption, viscosity coefficient μ = μ l +μ t ，μ l And mu t The laminar flow viscosity coefficient and the turbulent flow viscosity coefficient are respectively obtained by Sutherland formula
Wherein
The coefficient of turbulence viscosity mu for the dimensional incoming flow temperature t Given by a turbulence model; λ is given by Stokes assumption by Pr l The number of Prandtl is laminar flow; pr (Pr) of 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.
Further, by expressing the convection flux F in the Navier-Stokes equation as the non-viscous flux F, the formula (4) can be rewritten into the form of the formula (13):
in the formula (13), F 1 Solving by using a lattice boltzmann equation similar to BGK for the constant non-viscous flux; f 2 And taking the average value of grid physical quantities at two sides of the grid boundary as the physical quantity at the boundary for the flux change caused by the grid motion, and directly solving.
First, F is explained 2 Is solved by the physical quantity p at the grid boundary i-1/2 on the left side of the grid cell (i, j) 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 The physical quantity of the grids on the two sides of the boundary can be averaged to obtain:
velocity u of the grid motion at the boundary b,i-1/2,j ,v b,i-1/2,j The motion speed of the grid points at the two ends of the boundary is averaged to obtain:
so F in the formula (13) 2 The values of (A) are:
the lattice boltzmann equation approximated by BGK is:
in formula (14), r represents a physical position; τ represents the distribution time for the particle distribution to approach equilibrium by collision; f. of α Is the density distribution function along the alpha direction;
is a corresponding equilibrium state; delta t Is a flow time step; e.g. of the type α Is the velocity of the particle in the alpha direction; n is the velocity number of the discrete particles;
multiscale expansion by Chapman-Enskog [1] It can be demonstrated that: the formula (14) can be successfully restored to a macroscopic NS equation to establish a macroscopic conservation quantity W and a density distribution function f α The relationship between the micro physical quantity and the macro physical quantity; by solving equation (14) using the LB model, the density distribution function f can be obtained α Then according to the macroscopic conservation quantity W and the density distribution function f α The macroscopic conservation quantity W is obtained through the relation between the two, so that the inviscid flux F is obtained 1 。
Since most existing one-dimensional LB models contain a large number of user-specified parameters, these parameters have a large impact on the performance of LBFS. To overcome these deficiencies, yang et al [1，2] A non-free parameter D1Q4 model is provided for simulating non-viscous compressible flow with strong shock waves; solving a density distribution function f by adopting a non-free parameter D1Q4 model α In the non-free parameter D1Q4 model, the density distribution function f α The dispersion is solved in four directions, and the obtained result is a density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 The non-free parametric D1Q4 model is written as:
in formula (15) to formula (20), d 1 ，d 2 Is the grid speed; u is the average flow velocity in the one-dimensional case; c represents definition of
D represents a spatial dimension of the lattice boltzmann model, and for the non-free parameter D1Q4 model, D =1;
obtaining density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 And velocity d 1 ，d 2 Then, according to Chapman-Enskog multi-scale expansion, the following relationship between macroscopic physical quantities and microscopic physical quantities can be obtained:
in the formulae (21) to (25), ξ i Particle velocity in the i direction, xi 1 ＝d 1 ，ξ 2 ＝-d 1 ，ξ 3 ＝d 2 ，ξ 4 ＝-d 2 ；e p Is the potential energy of the particles and,
and (3) solving the inviscid flux after obtaining the macroscopic physical quantity, wherein the inviscid flux of the unit interface is as follows:
in the formula (26), the reaction mixture is,
is the cell interface at no viscous flux;
calculating the calculated non-viscous flux for the dissipation without taking into account the value;
calculating the obtained non-viscous flux for introducing numerical dissipation; alpha is alpha * In regions of boundary layer with low numerical dissipation, alpha, as a function of switching * The value of (a) tends to zero, so that a smaller numerical dissipation is introduced, and the area alpha with a larger numerical dissipation is ensured * Tends to 1, can accurately capture strong shock wave alpha * The specific way of calculating (c) is as follows:
α * ＝max{α L ,α R } (30)
formula (27) formula (30), tanh (x) is a hyperbolic tangent function; p is a radical of L And p R Is the pressure on the left and right sides of the cell interface; c is the amplification factor; as can be seen from formula (26), α ranges between 0 and 1; alpha is alpha L And alpha R Maximum values of left and right control volume switch functions respectively; n is a radical of L And N R The number of surfaces of the left and right control bodies.
Calculated inviscid flux without consideration of numerical dissipation
And introducing a non-viscous flux calculated from the numerical dissipation
The solving formulas of (1) are respectively as follows:
in formula (31) to formula (32), u = n x U n -n y U τ ，v＝n x U τ +n y U n ，U n Is the normal velocity, U, at the cell interface τ Is the tangential velocity at the cell interface, n x ，n y Is the component of the normal vector at the cell interface in the x, y direction.
In the above method, the function f is distributed due to density α The flux-increasing agent consists of a balance part and an unbalance part, wherein the balance part contributes to the non-viscous flux, and the unbalance part contributes to the viscous flux. Since LBM is only used to calculate the inviscid flux, it is not flatThe equilibrium portion can be regarded as numerical dissipation, and when calculating the non-viscous flux with smaller numerical dissipation, the non-equilibrium portion is not considered, so that very accurate results can be provided for boundary layer flow under the condition of introducing less numerical dissipation; however, when simulating hypersonic flow with strong shock waves, this format often shows oscillations or even divergence, and in order to accurately capture strong shock waves, numerical dissipation needs to be introduced, but large numerical dissipation affects the solution in smooth areas, resulting in inaccurate solution in the boundary layer, so that in order to accurately capture strong shock waves and thin boundary layers, we need to carefully control the numerical dissipation. In the region with small numerical dissipation near the boundary layer, the numerical dissipation is ignored when calculating the non-viscous flux, and the non-viscous flux is calculated at the moment
Can be solved by the following equation (31); in the region which is far from the boundary layer and is easy to generate strong shock waves and the region with larger numerical dissipation, the numerical dissipation needs to be introduced when calculating the inviscid flux, and the inviscid flux at the moment
Can be solved by the following equation (32); to integrate the advantages of both solutions, by introducing a switching function α * Control of
And
the specific gravity of (A) is solved by the inviscid flux, and the inviscid flux obtained at the moment is
The motion equation in step S3 refers to a motion equation of an elastic system, and specifically includes:
in formula (33), h and α are respectively a heave displacement and a pitch angle; m is a mass unit; s. the α Is the static moment about the elastic axis; i is α Is the moment of inertia; k h And K α Respectively heave and pitch spring constants; l is a lift force; m is the moment about the elastic axis;
the equation of motion of equation (33) has a length dimension of half chord length b and a natural frequency ω of uncoupled pitching motion α Dimensionless for the time dimension, the results are as follows:
in the formula (34), x α Is the static unbalance; omega h A non-coupled natural frequency of the elevating movement;
is the square of the radius of rotation; u shape * Is a dimensionless speed composed of
Defining; c l And C m Respectively, a lift coefficient and a moment coefficient;
due to the difference between the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the structure dimensionless time
Readjustment is required to ensure consistency of the overall system computation time scale,
wherein
Is flow field dimensionless time, L is length scale;
obtaining a lift coefficient C by solving the Navier-Stokes equation of the two-dimensional unsteady flow field in the step S2 l Coefficient of sum moment C m Updating the equation of motion in the dimensionless equation of motion in the formula (34), and decoupling the equation of motion in the dimensionless equation of motion (34) to obtain the following form:
wherein a is 1 、a 2 、b 1 、b 2 、c 1 And c 2 The equation coefficients after decoupling are all constants;
order to
Formula (35) is rewritten as:
second order discretization of equation (36) yields:
the equation of equation (37) is arranged in a matrix form to obtain:
the solution of the equation of motion can be obtained by directly carrying out matrix inversion.
The result of the motion equation in step S4 includes values of displacement, velocity, acceleration, angular velocity and angular acceleration of the motion of the object, and because the Navier-Stokes equation of the two-dimensional unsteady flow field and the motion equation adopt different dimensionless reference quantities, transformation is required, and the correspondence relationship is as follows:
reducing formula (39) to obtain formula (40):
in the formulae (39) to (40), h,
And
respectively displacement, speed, acceleration, rotation angular velocity and angular acceleration of an object in the vertical direction; h is physical 、h f And h s Respectively representing the displacement of an object in the vertical direction, dimensionless displacement in a flow field equation and dimensionless displacement in a motion equation;
and
respectively representing the velocity of the object in the vertical direction, the dimensionless velocity in a flow field equation and the dimensionless velocity in a motion equation;
and
acceleration of an object in the vertical direction, dimensionless acceleration in a flow field equation and dimensionless acceleration in a motion equation are respectively measured;
and
the rotation angular velocity of the object in the vertical direction, the dimensionless rotation angular velocity in the flow field equation and the dimensionless rotation angular velocity in the motion equation are respectively;
and
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 are respectively.
The flow field boundary conditions in step S4 include both ideal fluid and viscous fluid:
for an ideal fluid, the object plane adopts the condition of no penetration boundary, that is, the normal velocity of the fluid at the object plane is equal to the normal velocity of the object plane, such as:
for viscous fluids, the object plane adopts the non-slip boundary condition, that is, the velocity of the fluid at the object plane is equal to that of the object plane, such as:
in the formulae (41) to (42), x t And y t The speed of the object motion;
the acceleration of the object is determined by the motion rule of the object;
is the normal vector of the object plane.
The grid coordinates in step S4 are updated by using an infinite interpolation method (TFI), which specifically includes:
s41, firstly, parameterizing grid points in the flow field based on the arc length, taking an edge in the i direction as an example, and writing an arc length calculation formula as follows: s i ＝S i-1 +|r i+1,j -r i,j I, i =1,2, imax, then S k /S kmax K =1,2, \ 8230;, kmax, i.e. the parameterized values of each grid point on the edge, and the parameterized variables in the i, j directions are respectively denoted as SI i,j ,SJ i,j ；
S42, calculating the deformation quantity on each block vertex and each edge by utilizing a one-dimensional TFI technology, and knowing the grid deformation quantity dP of i =1 and i = imax at the vertex of each edge 1,j And dP imax,j Then, the grid deformation amount of any point on the edge is determined by equation (43):
dP i,j ＝(1-SI i,j )dP 1,j +SI i,j dP imax,j (43)
the grid deformation on other edges is calculated according to step S42;
s43, interpolating by using a two-dimensional TFI technology to obtain grid deformation of an internal point, and for a grid surface with imax grid points in the i direction and jmax grid points in the j direction, writing the deformation of any grid point on the surface as follows after the grid deformation on 4 edges on the surface is obtained:
wherein the bending function A i,j 、B i,j 、C i,j And D i,j The definition is as follows:
in formulae (44) to (45), dP 1,1 Is the displacement, dP, at the corner (1, 1) of the grid block 1,jmax Is the displacement, dP, at the grid block corner (1, jmax) imax,1 Is the displacement, dP, at the corner point (imax, 1) of the grid block imax,jmax Is the displacement, dP, at the grid block corner (imax, jmax) i,1 dP is the displacement at the vertex (i, 1) of the j-direction grid line where the grid point (i, j) is located i,jmax dP is the displacement at the vertex (i, jmax) of the j-direction grid line where the grid point (i, j) is located 1,j Is the displacement, dP, at the i-direction grid line vertex (1, j) where the grid point (i, j) is located imax,j Is the displacement at the i-direction grid line vertex (imax, j) where the grid point (i, j) is located; SI (Standard interface) i,1 Parameterization values of the j-direction grid line vertex (i, 1) where the grid point (i, j) is located on an edge whose vertex is the point (1, 1) and the point (imax, 1), SI i,jmax Parameterized values, SJ, of a j-direction grid line vertex (i, jmax) where the grid point (i, j) is located on an edge having the point (1, jmax) and the point (imax, jmax) as vertices 1,j Parameterized values of i-direction grid line vertex (1, j) where grid point (i, j) is located on an edge whose vertex is point (1, 1) and point (1, jmax), SJ imax,j And (ii) a parameterized value of an i-direction grid line vertex (imax, j) where the grid point (i, j) is located on an edge with the point (imax, 1) and the point (imax, jmax) as vertexes. Xi i,j ,η i,j Is constant, isCalculating an intermediate value without physical significance;
s44, loading the grid point displacement to the initial grid to generate a new flow field grid: p new ＝dP i,j +P original 。
In the fluid-solid coupling numerical simulation method, the fluid equation in the step S2 not only comprises the solution of a convection term, but also comprises the solution of a viscous flow term, the invention uses an SA turbulence model to calculate the turbulence viscosity coefficient, the SA turbulence model is a process model proposed by spaart and Allimaras, the model directly assumes that the turbulence viscosity coefficient meets a scalar equation (a convection-diffusion equation with a source term) in a flow field, and gives an equation and a coefficient based on experience; the model comprises a large number of empirical constants and empirical functions, and Spalart and the like have abundant experience in the aspect of aviation calculation and grasp abundant experimental (calculation) data, and the empirical constants and the empirical functions are well calibrated, so that the SA turbulence model has a good calculation effect in the aspect of aviation, and the solution of the viscous flow term is directly carried out by adopting the SA turbulence model.
Coefficient of turbulence viscosity mu given by SA turbulence model t Calculated using the following formula:
wherein
To pair
Limiting to prevent negative values from occurring to cause unstable calculation; l Ω is the modulus of vorticity, and d is the distance from the grid point to the wall surface; 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:
inflow conditions are as follows:
the exit conditions are extrapolated by making the values at the outer virtual grid points of the exit equal to the values at the inner points.
The Lattice Boltzmann Method (LBM) is an emerging mesoscopic numerical Method. In principle, LBM considers the migration and collision processes of particles near an interface at a mesoscopic level and reconstructs macroscopic physical quantities on the interface through an obtained particle velocity distribution function.
Compared with the prior art, the invention has the following advantages:
the fluid-solid coupling numerical simulation method introduces a lattice boltzmann method to solve the flow field without viscous flux, the lattice boltzmann method is a numerical simulation method taking contemporary statistical physics as a theoretical basis, and the complex physical phenomenon on a macroscopic level is simulated by tracking the interaction of a large number of discrete particles on a medium microscopic scale; the simulation method has the advantages of good stability, high calculation precision, clear physical significance and no consideration of nonlinear terms, and can improve the capture of flow field detail characteristics from the physical level, thereby effectively improving the aerodynamic simulation capability.
The fluid-structure interaction numerical simulation method has the innovation points that the lattice boltzmann flux method is applied to solving of the unsteady flow field, and the accuracy of the lattice boltzmann flux method in solving of the unsteady flow field is verified.
Drawings
FIG. 1 is a flow chart of the fluid-solid coupling numerical simulation method based on the lattice Boltzmann flux algorithm of the present invention;
FIG. 2 is a schematic view of a resiliently mounted airfoil of the present invention;
fig. 3 shows M =0.825, V * When =0.5, the airfoil profile vertical displacement and the rotation angle at any time when the damping response occurs
∞
The course of change in time;
fig. 4 shows M =0.825, V * When =0.51, the vertical displacement and the rotation angle of the airfoil profile at the time of neutral response are at any time
∞
The course of change in time;
fig. 5 shows M =0.825, V * When =0.65, the vertical displacement and rotation angle of the airfoil are dependent on the divergence response
∞
A process of change in time;
FIG. 6 shows the velocity index V at which neutral response occurs at different Mach numbers * 。
FIG. 7 shows the flutter frequency ratio ω/ω for neutral response at different Mach numbers α 。
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail below with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
As shown in fig. 1: a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm comprises the following steps:
s1, determining a calculation area, reading grid information, physical property parameters and control parameters, and entering a step S2;
s2, solving a Navier-Stokes equation of the two-dimensional unsteady flow field, wherein a convection term is solved by using a lattice Boltzmann method, a viscosity term is solved by using an SA turbulence model, then updating flow field parameters, and entering a step S3;
s3, solving a motion equation, and entering the step S4;
s4, updating the boundary conditions of the flow field and grid coordinates (the grid coordinates need to be updated due to the displacement caused by the motion of the object, and the boundary conditions of the flow field calculation are influenced by the speed and the acceleration of the object plane caused by the motion of the object) according to the result of solving the motion equation in the step S3, and entering a step S5;
s5, judging whether the simulation time is up, if so, entering a step S7, and if not, entering a step S6;
s6, returning to the step S2 to calculate the next time step;
and S7, ending the simulation.
The mesh information in step S1 includes: the number of grid blocks, the grid number of each grid block, the connection information of the grid blocks and grid coordinates, and the physical property parameters comprise: density, speed, temperature, pressure and viscosity coefficient of each grid point of the flow field, and control parameters comprise: the time step length, the CFL number, the inner iteration step number, the incoming flow Mach number, the Reynolds number and the incoming flow incidence angle; the flow field parameters in the step S2 are a lift coefficient and a moment coefficient; the result of solving the motion equation in the above step S4 includes displacement, velocity, acceleration, rotation angular velocity, and angular acceleration of the object.
The specific solving step of the Navier-Stokes equation of the two-dimensional unsteady flow field in the step S2 is as follows.
Wherein the Navier-Stokes equation is:
in the formula (1), Ω is a controller; s is a control volume element boundary surface;
external normal plane of S infinitesimalA product vector quantity; re is Reynolds number; w is macroscopic conservation quantity; f is the convection flux; f v Is the viscous flux;
macroscopic conservation quantity W, convection flux F, and viscous flux F in formula (1) v The specific expressions of (a) are respectively as follows:
viscous flux F v The expressions in (1) are as follows:
in the formula (2) -formula (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; μ is the kinetic viscosity coefficient, based on the vortex viscosity assumption, viscosity coefficient μ = μ l +μ t ，μ l And mu t The laminar flow viscosity coefficient and the turbulent flow viscosity coefficient are respectively obtained by Sutherland formula
Wherein
For dimensional incoming flow temperature, the turbulent viscosity coefficient mu t Given by a turbulence model; λ is given by Stokes assumption by Pr l The number of Prandtl is laminar flow; pr (Pr) of t Is the turbulent Prandtl number, pr l ＝0.72，Pr t =0.9; γ is the specific heat ratio, for air, γ =1.4; m ∞ The incoming flow mach number.
Expressing the convection flux F in the Navier-Stokes equation as a non-viscous flux F, equation (4) can be rewritten into the form of equation (13):
in formula (13), F 1 Solving by using a lattice boltzmann equation similar to BGK for the constant non-viscous flux; f 2 And taking the average value of grid physical quantities at two sides of the grid boundary as the physical quantity at the boundary for the flux change caused by the grid motion, and directly solving.
The non-viscous flux F at the time of setting will be described in detail below 1 Solving the aspect, the lattice boltzmann equation of the BGK approximation is:
in formula (14), r represents a physical position; τ represents the distribution time for the particle distribution to approach equilibrium by collision; f. of α Is the density distribution function along the alpha direction;
is a corresponding equilibrium state; delta t Is a flow time step; e.g. of the type α Is the velocity of the particle in the alpha direction; n is the velocity number of the discrete particles;
multiscale expansion by Chapman-Enskog [1] It can be demonstrated that: the formula (14) can be successfully restored to a macroscopic NS equation to establish a macroscopic conservation quantity W and a density distribution function f α The relationship between the micro physical quantity and the macro physical quantity; by solving equation (14) using the LB model, the density distribution function f can be obtained α Then according to the macroscopic conservation quantity W and the density distribution function f α The macroscopic conservation quantity W is obtained through the relation between the two, so that the inviscid flux F is obtained 1 。
The density distribution function f is solved by adopting a non-free parameter D1Q4 model α In the non-free parameter D1Q4 model, the density distribution function f α The dispersion is solved in four directions, and the obtained result is a density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 The non-free parametric D1Q4 model is written as:
in formula (15) -formula (20), d 1 ，d 2 Is the grid speed; u is the average flow velocity in the one-dimensional case; c represents definition as
D represents a spatial dimension of the lattice boltzmann model, and for the non-free parameter D1Q4 model, D =1;
obtaining density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 And velocity d 1 ，d 2 Then, according to Chapman-Enskog multi-scale expansion, the following relationship between macroscopic physical quantities and microscopic physical quantities can be obtained:
in the formulae (21) to (25), ξ i Particle velocity in the i direction, xi 1 ＝d 1 ，ξ 2 ＝-d 1 ，ξ 3 ＝d 2 ，ξ 4 ＝-d 2 ；e p Is the potential energy of the particles and,
and (3) solving the inviscid flux after obtaining the macroscopic physical quantity, wherein the inviscid flux of the unit interface is as follows:
in the formula (26), the reaction mixture is,
is the cell interface at no viscous flux;
calculating the calculated non-viscous flux for the dissipation without taking into account the value;
calculating the obtained non-viscous flux for introducing numerical dissipation; alpha is alpha * In regions of boundary layer with low numerical dissipation, alpha, as a function of switching * The value of (a) tends to zero, so that a smaller numerical dissipation is introduced, and the area alpha with a larger numerical dissipation is ensured * Tends to 1, can accurately capture strong shock wave alpha * The specific way of obtaining the value of (c) is as follows:
α * ＝max{α L ,α R } (30)
formula (27) formula (30), tanh (x) is a hyperbolic tangent function; p is a radical of formula L And p R Is the pressure on the left and right sides of the cell interface; c is the amplification factor, taken herein as C =100; as can be seen from formula (26), α ranges between 0 and 1; alpha is alpha L And alpha R Maximum values of left and right control volume switch functions respectively; n is a radical of L And N R The number of surfaces of the left and right control bodies.
Further, calculated inviscid flux without taking into account numerical dissipation
And introducing a non-viscous flux calculated from the numerical dissipation
The solving formulas of (1) are respectively as follows:
in formula (31) to formula (32), u = n x U n -n y U τ ，v＝n x U τ +n y U n ，U n Is the normal velocity, U, at the cell interface τ Is the tangential velocity at the cell interface, n x ，n y Is the component of the normal vector at the cell interface in the x, y direction.
And transferring the lift coefficient and the moment coefficient obtained by solving the flow field equation to a solving part of the object motion equation.
The following is a part of solving the motion equation, where the motion equation in step S3 refers to a motion equation of an elastic system, and specifically is:
in the formula (33), h and α are respectively a heave displacement and a pitch angle; m is a mass unit; s α Is the static moment about the elastic axis; i is α Is the moment of inertia; k h And K α Respectively heave and pitch spring constants; l is a lift force; m is the moment about the elastic axis;
the equation of motion of equation (33) has a half chord length b as a length dimension and a natural frequency ω of the uncoupled pitching motion α Dimensionless for the time dimension, the results are as follows:
in the formula (34), x α Is the static unbalance; omega h A non-coupled natural frequency of the elevating movement;
is the square of the radius of rotation; u shape * Is a dimensionless speed composed of
Defining; c l And C m Respectively, a lift coefficient and a moment coefficient;
due to the difference between the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the structure dimensionless time
Readjustment is required to ensure consistency of the overall system computation time scale,
wherein
Is flow field dimensionless time, L is length scale;
obtaining a lift coefficient C by solving the Navier-Stokes equation of the two-dimensional unsteady flow field in the step S2 l Coefficient of sum moment C m Update to equation (3)4) In the dimensionless equation of motion (34), the dimensionless equation of motion (34) is decoupled to obtain the following form:
wherein a is 1 、a 2 、b 1 、b 2 、c 1 And c 2 The equation coefficients after decoupling are all constants;
order to
Formula (35) is rewritten as:
second order discretization of equation (36) yields:
the equation of equation (37) is arranged in a matrix form to obtain:
the solution of the equation of motion can be obtained by directly carrying out matrix inversion. After solving the motion equation, updating motion parameters such as values of displacement, velocity, acceleration, rotation angular velocity and angular acceleration of the motion of the object, updating the boundary conditions of the flow field and grid coordinates based on the updated motion parameters, and starting or ending the solution of the next time step.
Because dimensionless reference quantities adopted by the Navier-Stokes equation and the motion equation of the two-dimensional unsteady flow field are different, transformation is required, and the corresponding relation is as follows:
reducing formula (39) to obtain formula (40):
in the formulae (39) to (40), h,
And
respectively displacement, speed, acceleration, rotation angular velocity and angular acceleration of an object in the vertical direction; h is physical 、h f And h s Respectively representing the displacement of an object in the vertical direction, dimensionless displacement in a flow field equation and dimensionless displacement in a motion equation;
and
respectively representing the velocity of the object in the vertical direction, the dimensionless velocity in a flow field equation and the dimensionless velocity in a motion equation;
and
acceleration of an object in the vertical direction, dimensionless acceleration in a flow field equation and dimensionless acceleration in a motion equation are respectively measured;
and
the rotation angular velocity of the object in the vertical direction, the dimensionless rotation angular velocity in the flow field equation and the dimensionless rotation angular velocity in the motion equation are respectively;
and
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 are respectively.
In the above method, the flow field boundary conditions in step S4 include both ideal fluid and viscous fluid: for an ideal fluid, the object plane adopts the condition of no penetration boundary, that is, the normal velocity of the fluid at the object plane is equal to the normal velocity of the object plane, such as:
for viscous fluids, the object plane adopts the non-slip boundary condition, that is, the velocity of the fluid at the object plane is equal to that of the object plane, such as:
in the formulae (41) to (42), x t And y t The speed of the object movement;
the acceleration of the object is determined by the motion rule of the object;
is the normal vector of the object plane.
In the above method, the grid coordinates in step S4 are updated by using an infinite interpolation method (TFI), which specifically includes the following steps:
s41, firstly, parameterizing grid points in the flow field based on the arc length, taking an edge in the i direction as an example, and writing an arc length calculation formula as follows: s i ＝S i-1 +|r i+1,j -r i,j I, i =1,2,. Imax, then S k /S kmax K =1,2, \ 8230;, kmax, i.e. the parameterized values of each grid point on the edge, and the parameterized variables in the i, j directions are respectively denoted as SI i,j ,SJ i,j ；
S42, calculating the deformation quantity of each block vertex and each block edge by utilizing a one-dimensional TFI technology, and knowing the deformation quantityMesh deformation dP at edge vertices i =1 and i = imax 1,j And dP imax,j Then, the grid deformation amount of any point on the edge is determined by equation (43):
dP i,j ＝(1-SI i,j )dP 1,j +SI i,j dP imax,j (43)
the grid deformation on other edges is calculated according to step S42;
s43, interpolating by using a two-dimensional TFI technology to obtain grid deformation of an internal point, and for a grid surface with imax grid points in the i direction and jmax grid points in the j direction, writing the deformation of any grid point on the surface as follows when the grid deformation on 4 edges on the surface is obtained:
wherein the bending function A i,j 、B i,j 、C i,j And D i,j The definition is as follows:
in formulae (44) to (45), dP 1,1 Is the displacement, dP, at the corner (1, 1) of the grid block 1,jmax Is the displacement, dP, at the grid block corner (1, jmax) imax,1 Is the displacement, dP, at the corner point (imax, 1) of the grid block imax,jmax Is the displacement at the corner point (imax, jmax) of the grid block, dP i,1 Is the displacement, dP, at the vertex (i, 1) of the grid line in the j direction at which the grid point (i, j) is located i,jmax Is the displacement, dP, at the j-direction grid line vertex (i, jmax) where the grid point (i, j) is located 1,j Is the displacement, dP, at the i-direction grid line vertex (1, j) where the grid point (i, j) is located imax,j Is the displacement at the i-direction grid line vertex (imax, j) where the grid point (i, j) is located; SI (Standard interface) i,1 Parameterization values of the j-direction grid line vertex (i, 1) where the grid point (i, j) is located on an edge whose vertex is the point (1, 1) and the point (imax, 1), SI i,jmax The parameter of the vertex (i, jmax) of the grid line in the j direction where the grid point (i, j) is located on the edge having the point (1, jmax) and the point (imax, jmax) as the verticesDigitized value, SJ 1,j Parameterization values of i-direction grid line vertexes (1, j) of grid points (i, j) on edges taking the points (1, 1) and (1, jmax) as vertexes, SJ imax,j Parameterised values of the i-direction grid line vertex (imax, j) where the grid point (i, j) is located on an edge whose vertex is the point (imax, 1) and the point (imax, jmax). Xi i,j ,η i,j The constant is a middle value, and has no physical significance;
s44, loading grid point displacement to the initial grid to generate a new flow field grid: p new ＝dP i,j +P original 。
Verification example
To verify the correctness of the algorithm of the present invention, the following example verification was performed, and the model information used is shown in fig. 2. The wing profile moves in pitching and vertical lifting around a given elastic axis by adopting a structural model of two-dimensional sweepback wing vibration with the section of NACA64A 010. The pitch axis is defined by the distance a, which is the distance between the elastic axis and the mid-chord line as a percentage of half the chord length. If a is a positive number, it indicates that the axis is downstream of the mid-chord point, and a negative number indicates that it is upstream of the mid-chord point. The structural parameters adopted by the model are shown in table 1.
TABLE 1
Due to the symmetrical profile of the NACA64a010 airfoil, an initial perturbation is applied to trigger the oscillatory motion. The wing profile makes sinusoidal motion at a natural pitch frequency omega α Pitching about its elastic axis, amplitude alpha 0 =1 °. The forced pitch mode typically lasts 1-3 cycles. On the basis, the elastic installation wing profile is set to freely move in the vertical direction and the pitching direction, and dynamic response is recorded.
FIGS. 3-5 show M ∞ When =0.825, three different V * The time course of the vertical displacement and the rotation angle of the wing profile. These figures correspond to the damped, neutral and divergent responses, respectively. The main task of computing the flutter boundaries is to determine the occurrence of neutral sounds by observing these imagesCritical position of stress, when V * When the value of (a) is less than the value of the flutter boundary, a damping response occurs, and the changes in vertical displacement and pitch angle attenuate with time, as shown in fig. 3; when V is * When the value of (a) is close to the value of the flutter boundary, a neutral response occurs, and the vertical displacement and the pitch angle are not changed substantially with time, as shown in fig. 4; when V is * When the value of (a) is greater than the value of the flutter boundary, a divergent response occurs, and the changes in vertical displacement and pitch angle increase with time, as shown in fig. 5.
Then calculating the speed index V of neutral response when different Mach numbers are calculated * And flutter frequency ratio omega/omega α The results are shown in FIGS. 6 and 7. Fig. 6 and 7 show the calculated flutter speed boundary and flutter frequency ratio boundary, which are well matched with the research results of the former. The flutter boundary appears as a "double step" shape at M ∞ =0.83 and M ∞ Two jumps in flutter speed occur at 0.91. The flutter boundary can be roughly divided into 4 regions: subsonic region: m ∞ Less than 0.80; transonic pits: m ∞ = 0.80-0.83; first flutter speed jump region: m ∞ =0.84 to 0.91; locking area: m ∞ Is greater than 0.91. The results herein differ from the previous results primarily due to the turbulence model and reynolds number used for the simulation. In this study, all calculations used fixed Re =1.256 × 10 6 。
The innovation point of the method is that the lattice boltzmann flux method is applied to solving the unsteady flow field, and the calculation example proves the accuracy of the lattice boltzmann flux method in the aspect of solving the unsteady flow field.
The fluid-solid coupling numerical simulation method introduces a lattice boltzmann method to solve the flow field without viscous flux, the lattice boltzmann method is a numerical simulation method which takes contemporary statistical physics as a theoretical basis, and the complex physical phenomenon on a macroscopic level is simulated by tracking the interaction of a large number of discrete particles on the microscopic scale of a medium. The simulation method has the advantages of good stability, high calculation precision, clear physical significance and no consideration of nonlinear terms, and can improve the capture of flow field detail characteristics from the physical level, thereby effectively improving the aerodynamic simulation capability.
The above embodiments are only for illustrating the technical idea of the present invention, and the protection scope of the present invention cannot be limited thereby, and any modification made on the basis of the technical scheme according to the technical idea proposed by the present invention falls within the protection scope of the present invention; the technology not related to the invention can be realized by the prior art.
Reference documents:
[1]He X Y,Luo L S.Lattice Boltzmann model for the incompressible Navier-Stokes equation[J].Journal of Statistical Physics,1997,88:927-944.
[2]L.M.Yang,C.Shu,J.Wu.Development and comparative studies of three non-free parameter lattice Boltzmann.
Claims (10)
1. a fluid-solid coupling numerical simulation method based on a lattice Boltzmann flux algorithm is characterized in that: the method comprises the following steps:
s1, determining a calculation area, reading grid information, physical property parameters and control parameters, and entering a step S2;
s2, solving a Navier-Stokes equation of the two-dimensional unsteady flow field, wherein a convection term is solved by using a lattice Boltzmann method, then updating flow field parameters, and entering a step S3;
s3, solving a motion equation, and entering the step S4;
s4, updating the boundary conditions of the flow field and grid coordinates according to the result of solving the motion equation in the step S3, and entering a step S5;
s5, judging whether the simulation time is up, if so, entering a step S7, and if not, entering a step S6;
s6, returning to the step S2 to calculate the next time step;
and S7, ending the simulation.
2. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 1, wherein: the mesh information in step S1 includes: the number of grid blocks, the grid number of each grid block, the connection information of the grid blocks and grid coordinates, and the physical property parameters comprise: density, speed, temperature, pressure and viscosity coefficient of each grid point of the flow field, and control parameters comprise: the time step length, the CFL number, the inner iteration step number, the incoming flow Mach number, the Reynolds number and the incoming flow incidence angle; and the flow field parameters in the step S2 are a lift coefficient and a moment coefficient.
3. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 1, wherein: the Navier-Stokes equation in the step S2 is as follows:
in the formula (1), Ω is a controller; s is a control volume element boundary surface;
is the external normal area vector of the S infinitesimal; re is Reynolds number; w is macroscopic conservation quantity; f is convection flux; f v Is the viscous flux;
macroscopic conservation quantity W, convection flux F, and viscous flux F in formula (1) v The specific expressions of (a) are respectively as follows:
viscous flux F v The expressions in (1) are as follows:
in the formula (2) -formula (12), ρ, E, p, and T represent density, total energy per unit mass, pressure, and temperature, respectively; τ is the viscous shear stress; μ is the kinetic viscosity coefficient, based on the vortex viscosity assumption, viscosity coefficient μ = μ l +μ t ，μ l And mu t The laminar flow viscosity coefficient and the turbulent flow viscosity coefficient are respectively obtained by Sutherland formula
Wherein
For dimensional incoming flow temperature, the turbulent viscosity coefficient mu t Given by a turbulence model; λ is given by Stokes assumption by Pr l Is laminar Prandtl number; pr (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.
4. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 2, wherein: expressing the convection flux F in the Navier-Stokes equation as a non-viscous flux F, equation (4) can be rewritten into the form of equation (13):
in the formula (13), F 1 Solving by using a lattice boltzmann equation similar to BGK for the constant non-viscous flux; f 2 And taking the average value of the grid physical quantities at two sides of the grid boundary as the physical quantity at the boundary for the flux change caused by the grid motion, and directly solving.
5. The lattice boltzmann flux algorithm-based fluid-structure interaction numerical simulation method of claim 4, wherein: the lattice boltzmann equation for BGK approximation is:
in formula (14), r represents a physical position; τ represents a radical passing throughThe distribution time at which the collision drives the particle distribution towards an equilibrium state; f. of α Is the density distribution function along the alpha direction;
is a corresponding equilibrium state; delta. For the preparation of a coating t Is a flow time step; e.g. of a cylinder α Is the velocity of the particle in the alpha direction; n is the velocity number of the discrete particles;
the following can be demonstrated by Chapman-Enskog multiscale expansion: the formula (14) can be successfully restored to a macroscopic NS equation to establish a macroscopic conservation quantity W and a density distribution function f α The relationship between the micro physical quantity and the macro physical quantity; by solving equation (14) using the LB model, the density distribution function f can be obtained α Then according to the macroscopic conservation quantity W and the density distribution function f α The macroscopic conservation quantity W is obtained through the relation between the two, so that the inviscid flux F is obtained 1 。
6. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 5, wherein: solving a density distribution function f by adopting a non-free parameter D1Q4 model α In the non-free parameter D1Q4 model, the density distribution function f α The solution is carried out when the dispersion is in four directions, and the obtained result is a density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 The non-free parametric D1Q4 model is written as:
in formula (15) -formula (20), d 1 ，d 2 Is the grid speed; u is the average flow velocity in the one-dimensional case; c represents definition of
D represents a spatial dimension of the lattice boltzmann model, and for the non-free parameter D1Q4 model, D =1;
obtaining density distribution function g in four directions 1 ，g 2 ，g 3 ，g 4 And velocity d 1 ，d 2 Then, according to Chapman-Enskog multi-scale expansion, the following relationship between macroscopic physical quantities and microscopic physical quantities can be obtained:
in the formulae (21) to (25), ξ i Velocity of the particle in the i direction, ξ 1 ＝d 1 ，ξ 2 ＝-d 1 ，ξ 3 ＝d 2 ，ξ 4 ＝-d 2 ；e p Is the potential energy of the particles and,
and (3) solving the inviscid flux after obtaining the macroscopic physical quantity, wherein the inviscid flux of the unit interface is as follows:
in the formula (26), the reaction mixture is,
is the cell interface at no viscous flux;
calculating the calculated non-viscous flux for the dissipation without taking into account the value;
calculating the obtained non-viscous flux for introducing numerical dissipation; alpha is alpha * In regions of boundary layer with low numerical dissipation, alpha, as a switching function * The value of (a) tends to zero, so that a smaller numerical dissipation is introduced, and the area alpha with a larger numerical dissipation is ensured * Tend to 1, can
Capturing strong shock wave, alpha, with sufficient accuracy * The specific way of obtaining the value of (c) is as follows:
α * ＝max{α L ,α R }(30)
formula (27) — formula (30), tanh (x) is a hyperbolic tangent function; p is a radical of L And p R Is the pressure on the left and right sides of the cell interface; c is the amplification factor; as can be seen from formula (26), α ranges between 0 and 1; alpha (alpha) ("alpha") L And alpha R Maximum values of left and right control volume switch functions respectively; n is a radical of L And N R The number of surfaces of the left and right control bodies.
7. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 6, wherein: calculated inviscid flux without consideration of numerical dissipation
And introducing a non-viscous flux calculated from the numerical dissipation
The solving formulas of (1) are respectively as follows:
in formula (31) to formula (32), u = n x U n -n y U τ ，v＝n x U τ +n y U n ，U n Is the normal velocity, U, at the cell interface τ Is the tangential velocity at the cell interface, n x ，n y Is the component of the normal vector at the cell interface in the x, y direction.
8. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 1, wherein: the motion equation in step S3 refers to a motion equation of an elastic system, and specifically includes:
in the formula (33), h and α are respectively a heave displacement and a pitch angle; m is a mass unit; s α Is the static moment about the elastic axis; i is α Is the moment of inertia; k h And K α Respectively heave and pitch spring constants; l is a lift force; m is the moment about the elastic axis;
the equation of motion of equation (33) has a half chord length b as a length dimension and a natural frequency ω of the uncoupled pitching motion α Dimensionless for the time dimension, the results are as follows:
in the formula (34), x α Is the static unbalance; omega h A non-coupled natural frequency of the elevating movement;
is the square of the radius of rotation; u shape * Is a dimensionless speed composed of
Defining; c l And C m Respectively, a lift coefficient and a moment coefficient;
due to the difference between the dimensionless time dimension of the equation of motion (34) and the fluid control equation, the structure dimensionless time
Readjustment is required to ensure consistency of the overall system computation time scale,
wherein
Is flow field dimensionless time, L is length scale;
obtaining a lift coefficient C by solving the Navier-Stokes equation of the two-dimensional unsteady flow field in the step S2 l Coefficient of sum moment C m Updated into the dimensionless equation of motion of equation (34),
decoupling the dimensionless equation of motion (34) to obtain the following form:
wherein a is 1 、a 2 、b 1 、b 2 、c 1 And c 2 The equation coefficients after decoupling are all constants;
order to
Formula (35) is rewritten as:
second order discretization of equation (36) yields:
the equation of equation (37) is arranged in a matrix form to obtain:
the solution of the equation of motion can be obtained by directly carrying out matrix inversion.
9. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 1, wherein: the result of the motion equation in step S4 includes values of displacement, velocity, acceleration, angular velocity and angular acceleration of the motion of the object, and because the Navier-Stokes equation of the two-dimensional unsteady flow field and the motion equation adopt different dimensionless reference quantities, transformation is required, and the correspondence relationship is as follows:
reducing formula (39) to obtain formula (40):
in the formulae (39) to (40), h,
And
respectively displacement, speed, acceleration, rotation angular velocity and angular acceleration of an object in the vertical direction; h is physical 、h f And h s Respectively the displacement of an object in the vertical direction, dimensionless displacement and motion in a flow field equationDimensionless displacement in the equation;
and
respectively representing the velocity of the object in the vertical direction, the dimensionless velocity in a flow field equation and the dimensionless velocity in a motion equation;
and
acceleration of an object in the vertical direction, dimensionless acceleration in a flow field equation and dimensionless acceleration in a motion equation are respectively measured;
and
the rotation angular velocity of the object in the vertical direction, the dimensionless rotation angular velocity in the flow field equation and the dimensionless rotation angular velocity in the motion equation are respectively;
and
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 are respectively.
10. The lattice boltzmann flux algorithm-based fluid-solid coupling numerical simulation method of claim 1 or 9, wherein: the flow field boundary conditions in step S4 include both ideal fluid and viscous fluid:
for an ideal fluid, the object plane adopts the condition of no penetration boundary, that is, the normal velocity of the fluid at the object plane is equal to the normal velocity of the object plane, such as:
for viscous fluids, the object plane adopts the non-slip boundary condition, that is, the velocity of the fluid at the object plane is equal to that of the object plane, such as:
in the formulae (41) to (42), x t And y t The speed of the object motion;
the acceleration of the object is determined by the motion rule of the object;
is the normal vector of the object plane;
the grid coordinates in step S4 are updated by using an infinite interpolation method (TFI), which includes the following steps:
s41, firstly, parameterizing grid points in the flow field based on the arc length, taking an edge in the i direction as an example, and writing an arc length calculation formula as follows: s i ＝S i-1 +|r i+1,j -r i,j I, i =1,2, imax, then S k /S kmax K =1,2, \ 8230;, kmax, i.e. the parameterized values of each grid point on the edge, and the parameterized variables in the i, j directions are respectively denoted as SI i,j ,SJ i,j ；
S42, calculating the deformation quantity on each block vertex and each edge by utilizing a one-dimensional TFI technology, and knowing the grid deformation quantity dP of i =1 and i = imax at the vertex of each edge 1,j And dP imax,j Then, the grid deformation amount of any point on the edge is determined by equation (43):
dP i,j ＝(1-SI i,j )dP 1,j +SI i,j dP imax,j (43)
the grid deformation on other edges is calculated according to step S42;
s43, interpolating by using a two-dimensional TFI technology to obtain grid deformation of an internal point, and for a grid surface with imax grid points in the i direction and jmax grid points in the j direction, writing the deformation of any grid point on the surface as follows after the grid deformation on 4 edges on the surface is obtained:
wherein the bending function A i,j 、B i,j 、C i,j And D i,j The definition is as follows:
in formulae (44) to (45), dP 1,1 Is the displacement, dP, at the corner (1, 1) of the grid block 1,jmax Is the displacement, dP, at the grid block corner (1, jmax) imax,1 Is the displacement, dP, at the corner point (imax, 1) of the grid block imax,jmax Is the displacement, dP, at the grid block corner (imax, jmax) i,1 dP is the displacement at the vertex (i, 1) of the j-direction grid line where the grid point (i, j) is located i,jmax Is the displacement, dP, at the j-direction grid line vertex (i, jmax) where the grid point (i, j) is located 1,j Is the displacement, dP, at the i-direction grid line vertex (1, j) where the grid point (i, j) is located imax,j Is the displacement at the i-direction grid line vertex (imax, j) where the grid point (i, j) is located; SI (Standard interface) i,1 Parameterization values of the j-direction grid line vertex (i, 1) where the grid point (i, j) is located on an edge whose vertex is the point (1, 1) and the point (imax, 1), SI i,jmax Parameterized values, SJ, of a j-direction grid line vertex (i, jmax) where the grid point (i, j) is located on an edge having the point (1, jmax) and the point (imax, jmax) as vertices 1,j Parameterization values of the vertices (1, j) of the i-direction grid lines where the grid points (i, j) are located on the edges with the vertices (1, 1) and (1, jmax),SJ imax,j parameterised values of the i-direction grid line vertex (imax, j) where the grid point (i, j) is located on an edge whose vertex is the point (imax, 1) and the point (imax, jmax). Xi i,j ,η i,j The constant is a middle value, and has no physical significance;
s44, loading the grid point displacement to the initial grid to generate a new flow field grid: p new ＝dP i,j +P original 。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211546774.5A CN115828782A (en) | 2022-12-05 | 2022-12-05 | Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211546774.5A CN115828782A (en) | 2022-12-05 | 2022-12-05 | Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115828782A true CN115828782A (en) | 2023-03-21 |
Family
ID=85545113
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211546774.5A Pending CN115828782A (en) | 2022-12-05 | 2022-12-05 | Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115828782A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116227388A (en) * | 2023-04-21 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
-
2022
- 2022-12-05 CN CN202211546774.5A patent/CN115828782A/en active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116227388A (en) * | 2023-04-21 | 2023-06-06 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
| CN116227388B (en) * | 2023-04-21 | 2023-07-07 | 中国空气动力研究与发展中心计算空气动力研究所 | Dynamic adjustment method, system, equipment and medium for CFL number of high-ultra-flow simulation |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yang et al. | Unstructured dynamic meshes with higher-order time integration schemes for the unsteady Navier-Stokes equations | |
| Sváček et al. | Numerical simulation of flow induced airfoil vibrations with large amplitudes | |
| CN104298869B (en) | A kind of fluid structurecoupling Numerical prediction method of elastic hydrofoil | |
| Sadeghi et al. | Application of three-dimensional interfaces for data transfer in aeroelastic computations | |
| Petinrin et al. | Computational study of aerodynamic flow over NACA 4412 airfoil | |
| CN115828782A (en) | Fluid-solid coupling numerical simulation method based on lattice Boltzmann flux algorithm | |
| Yang et al. | Higher-order time integration schemes for aeroelastic applications on unstructured meshes | |
| Qiao et al. | Nonlinear aeroelastic characteristics analysis of composite wing with high aspect ratio based on co-rotational method | |
| Ritter et al. | Dynamic aeroelastic simulations of the Pazy wing by UVLM with nonlinear viscous corrections | |
| Kuzmina et al. | Analysis of static and dynamic aeroelastic characteristics of airplane in transonic flow | |
| Dubcová et al. | Numerical simulation of interaction between turbulent flow and a vibrating airfoil | |
| Lin et al. | Discontinuous Galerkin finite element method for Euler and Navier-Stokes equations | |
| Cavagna et al. | Efficient application of CFD aeroelastic methods using commercial software | |
| Drofelnik et al. | Rapid calculation of unsteady aircraft loads | |
| Guruswamy et al. | Transonic aeroelastic analysis of the B-1 wing | |
| Geller et al. | An explicit model for three-dimensional fluid-structure interaction using LBM and p-FEM | |
| Streher | Large-eddy simulations of the flow around a NACA0012 airfoil at different angles of attack | |
| Morton et al. | Nonlinear analysis of airfoil flutter at transonic speeds | |
| Cavallaro et al. | Phenomenology of nonlinear aeroelastic responses of highly deformable joined-wings configurations | |
| Lahooti et al. | Thick strip method for efficient large-eddy simulations of flexible wings in stall | |
| Djayapertapa et al. | Two‐dimensional transonic aeroservoelastic computations in the time domain | |
| Athavale et al. | Application of an unstructured grid solution methodology to turbomachinery flows | |
| Rougeault-Sens et al. | Numerical unsteady aerodynamics for turbomachinery aeroelasticity | |
| Relvas et al. | Aeroelasticity of nonlinear structures using the corotational method | |
| Muffo et al. | Interface Velocity Consistency in time-accurate flow simulations on dynamic meshes |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |