A method and device for simulating two-phase flow in an aero-engine bearing cavity

Abstract

This application provides a two-phase flow simulation method and apparatus for an aero-engine bearing cavity, relating to the field of aero-engine simulation technology. The method is applied to the geometric model of an aero-engine bearing cavity, which includes a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet. The method includes: meshing the geometric model of the aero-engine bearing cavity to obtain a discretized Cartesian mesh; setting initial and boundary conditions for the computational domain of the discretized Cartesian mesh; and performing flow field simulation calculations on the geometric model of the aero-engine bearing cavity based on the set initial and boundary conditions to obtain the physical quantity distribution of the mixed flow field at each mesh point at the current moment. This application can accurately obtain the distribution pattern of lubricating oil in the engine bearing cavity with fewer meshes and less computation in simulation calculations, saving computation time.

Description

Technical Field

[0001] This application relates to the field of aero-engine simulation technology, and in particular to a method and apparatus for simulating two-phase flow in aero-engine bearing cavities. Background Technology

[0002] As a component of aero engines, the lubrication system primarily performs the functions of lubrication and cooling, thereby ensuring the engine's long-term, highly reliable, and stable operation. Currently, the aero engine industry is still developing rapidly, with increasingly harsher operating environments and higher performance requirements for engines, which also places more stringent demands on the lubrication system.

[0003] During the operation of an aero-engine, the lubricating oil in the bearing cavity interacts with the sealing gas inside the cavity under the high-speed rotation of the bearing, forming a complex oil-gas two-phase flow distribution. Since this has a significant impact on the lubrication design and heat transfer analysis of the bearing cavity, studying the oil-gas two-phase flow inside the cavity can improve the reliability and accuracy of the research.

[0004] The simulation of two-phase flow in the bearing cavity of aero-engines involves deformation problems such as interface migration and fragmentation. Traditional computational fluid dynamics methods, such as VOF (Volume of Fluid) and LevelSet, require tracking or capturing the interface, which increases the complexity of the computational model for two-phase flow problems and cannot consider the interactions between molecules from a microscopic perspective. Summary of the Invention

[0005] In view of this, this application provides a two-phase flow simulation method and apparatus for aero-engine bearing cavities to solve the above-mentioned technical problems.

[0006] In a first aspect, embodiments of this application provide a two-phase flow simulation method for an aero-engine bearing cavity, applied to a geometric model of the aero-engine bearing cavity. The geometric model includes a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet. The method includes:

[0007] The geometric model of the aero-engine bearing cavity is meshed to obtain a discretized Cartesian mesh;

[0008] Set the initial and boundary conditions for the computational domain of the discretized Cartesian mesh;

[0009] Based on the set initial and boundary conditions, flow field simulation calculations are performed on the geometric model of the aero-engine bearing cavity to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment.

[0010] Furthermore, the initial and boundary conditions of the computational domain of the discretized Cartesian mesh are set, including:

[0011] Set the distribution of lubricating oil and air in the computational domain at the initial moment;

[0012] Obtain the density and velocity at the lubricating oil inlet, and set the boundary conditions for the distribution function at the lubricating oil inlet:

[0013]

[0014] Where t represents the current time, x r ρ represents the coordinates at the boundary of the lubricating oil inlet. r U represents the density at the lubricating oil inlet boundary. r The velocity vector at the boundary of the lubricating oil inlet, x rf The coordinates of the fluid domain immediately adjacent to the lubricating oil inlet boundary, u rf ρ represents the velocity vector of the fluid domain immediately adjacent to the lubricating oil inlet boundary. rf f represents the density of the fluid domain immediately adjacent to the lubricating oil inlet boundary. r (x r (t) represents x at time t. r Density distribution function at that location, This represents the equilibrium state of the density distribution function;

[0015] Obtain the density and pressure at the lubricating oil outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions of the distribution function at the lubricating oil outlet.

[0016] Obtain the density and velocity at the air inlet, and set the boundary conditions for the density distribution function at the air inlet:

[0017]

[0018] Where, x a ρ represents the coordinates at the air inlet boundary. a U represents the density at the air inlet boundary. a The velocity vector at the air inlet boundary, x af The coordinates of the fluid domain immediately adjacent to the air inlet boundary, u af ρ represents the velocity vector of the fluid domain immediately adjacent to the air inlet boundary. af f represents the density of the fluid domain immediately adjacent to the air inlet boundary. a (x r (t) represents x at time t. r The air density distribution function at that location, This represents the equilibrium state of the air density distribution function;

[0019] Obtain the density and pressure at the air outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions for the distribution function at the air outlet.

[0020] The boundary condition for the density distribution function of the mixed fluid at the wall boundary is:

[0021]

[0022] Where, x b The coordinates at the wall boundary are represented by Δt, where Δt represents the time step, and u is the coordinates at the wall boundary. w This represents the velocity vector at the wall boundary. When the wall is a stationary wall of the cavity, u w =0, ρ w w represents the density at the wall boundary. k e represents the weighting coefficient. k Let c be the vector representing the direction of the k-th discrete velocity. s For the speed of sound, f k (x b (t) represents x at time t. b The density distribution function at point k, where k represents the index of the direction of the distribution function. It represents the distribution function of the velocity direction opposite to the k-th velocity direction after the collision at time t-Δt.

[0023] Furthermore, based on the set initial and boundary conditions, flow field simulation calculations are performed on the geometric model of the aero-engine bearing cavity to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment; including:

[0024] Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, calculate the vector f of the lubricating oil density distribution function at grid point x at the current time t. r , where f r =(f r，1 f r，2， …f r，K K is the number of discrete velocity directions;

[0025] Calculate the lubricating oil density ρ at grid point x at the current time t. r (t, x):

[0026]

[0027] Calculate the lubricating oil velocity u at grid point x at the current time t. r (t, x):

[0028]

[0029] Among them, e k It is the vector representing the direction of the k-th discrete velocity;

[0030] Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, the vector f of the air density distribution function at grid point x at the current time t is calculated. a , where f a =(f a，1 f a，2， …f a，K );

[0031] Calculate the air density ρ at grid point x at time t. a (t, x):

[0032]

[0033] Calculate the air velocity u at grid point x at the current time t. a (t, x):

[0034]

[0035] Then the mixed fluid density ρ(t, x) at grid point x at the current time t is:

[0036] ρ(t, x) = ρ r (t, x) + ρ a (t, x)

[0037] Then the velocity u(t, x) of the mixed fluid at grid point x at the current time t is:

[0038]

[0039] Among them, F r，total For the combined force of the lubricating oil, F a，total It is the resultant force of air.

[0040] Furthermore, based on the evolution equation of the LBM distribution function of the multi-relaxation collision operator, the vector f of the density distribution function of the lubricating oil at grid point x at the current time t is calculated. r ;include:

[0041] Vector f r The k-th component f r，k The calculation formula is:

[0042]

[0043] Among them, f r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Density distribution function at Δt; M -1 S is the inverse of the transformation matrix M; r Let m be the relaxation matrix. r(xe k Δt, t-Δt) are the density distribution function f r (xe k The moments of Δt, t-Δt), Let f be the density distribution function r (xe k The equilibrium moment F of Δt, t-Δt) r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external force of the lubricating oil at point Δt;

[0044] Lubricating oil external force F r (xe k The formula for calculating Δt, t-Δt) is:

[0045]

[0046] Where, ρ r (t-Δt, x) represents the lubricating oil density at grid point x at the previous time t-Δt; u eq Let x be the equilibrium velocity of the mixed fluid at grid point x at the previous time t-Δt:

[0047]

[0048] Lubricating oil drift speed Δu r The calculation formula is:

[0049] Δu r =F r,total Δt / ρ r (t-Δt, x)

[0050] Among them, the combined force of lubricating oil F r,total for:

[0051]

[0052] The force F between lubricating oil and air ra for:

[0053]

[0054] Among them, g ra w represents the strength of the interaction between the lubricating oil and air. k The weighting coefficient for the k-th discrete velocity direction;

[0055] The force between the lubricating oil and the wall surface for:

[0056]

[0057] in, The strength of the interaction between the lubricating oil and the wall surface, s(xe) k Δt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0;

[0058] external forces of gravity for:

[0059]

[0060] Where g is the gravity coefficient.

[0061] Furthermore, based on the LBM distribution function evolution equation of the multiple relaxation collision operator, the vector f of the air density distribution function at grid point x at the current time t is calculated. a ;include:

[0062] Vector f a The k-th component f a，k The calculation formula is:

[0063]

[0064] Among them, f a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Air density distribution function at Δt; S a Let m be the relaxation matrix. a (xe k Δt, t-Δt) are the density distribution function f a (xe k The moments of Δt, t-Δt), Let f be the density distribution function a (xe k The equilibrium moment F of Δt, t-Δt) a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external air force at point Δt:

[0065] External air force F a (xe k The formula for calculating Δt, t-Δt) is:

[0066]

[0067] Where, ρ a (t-Δt, x) represents the air density at grid point x at the previous time t-Δt.

[0068] air drift velocity Δu a The calculation formula is:

[0069] Δu a =F a，total Δt / ρ a (t-Δt, x)

[0070] Among them, the resultant force of air F a，total for:

[0071]

[0072] The force F between air and lubricating oil ar for:

[0073]

[0074] Among them, g ar The strength of the interaction between air and lubricating oil, w k The weighting coefficient for the k-th discrete velocity direction;

[0075] Forces between air and wall for:

[0076]

[0077] in, The intensity of the interaction between the air and the wall, s(xe k Δt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0;

[0078] external forces of gravity for:

[0079]

[0080] Furthermore, the method also includes:

[0081] Based on the relationship between the relaxation factor τ and the time step Δt, spatial step Δx, and kinematic viscosity v: The range of the relaxation factor is 0.5 < τ < 2, which determines the values ​​of step size Δt, spatial step size Δx, and kinematic viscosity v.

[0082] Furthermore, after meshing the geometric model of the aero-engine bearing cavity to obtain the discretized Cartesian mesh, the following steps are also included:

[0083] Calculate the magnitude of the phase volume fraction gradient at the phase interface

[0084] Using the magnitude of gradient Determine the encryption level factor A:

[0085]

[0086] Where A0 is the encryption level factor constant, [] represents integer operations, tanh(·) is the hyperbolic tangent function, and C α The reference value for the magnitude of the gradient is set to determine the degree of encryption based on the magnitude of different gradients;

[0087] For a square grid in the computational domain, perform 2 A This process involves multiple divisions to achieve mesh encryption at the interface.

[0088] Furthermore, after meshing the geometric model of the aero-engine bearing cavity to obtain the discretized Cartesian mesh, the following steps are also included:

[0089] The modulus for calculating the gradient of phase volume fraction within the gas or liquid phase. When the magnitude of the gradient When the value is less than the set coarsening threshold, the mesh is merged into two layers, thereby achieving mesh coarsening within the gas or liquid phase.

[0090] Secondly, embodiments of this application provide a two-phase flow simulation device for an aero-engine bearing cavity, applied to a geometric model of an aero-engine bearing cavity. The geometric model includes a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet. The device includes:

[0091] Mesh generation unit is used to mesh the geometric model of the bearing cavity of an aero-engine to obtain a discretized Cartesian mesh;

[0092] The cell is used to set the initial and boundary conditions of the computational domain of the discretized Cartesian mesh:

[0093] The computing unit is used to perform flow field simulation calculations on the geometric model of the aero-engine bearing cavity based on the set initial and boundary conditions, and to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment.

[0094] Thirdly, embodiments of this application provide an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method of embodiments of this application.

[0095] Fourthly, embodiments of this application provide a computer-readable storage medium storing computer instructions that, when executed by a processor, implement the methods of embodiments of this application.

[0096] This application can accurately obtain the distribution pattern of lubricating oil in the engine bearing cavity with less grid and computational cost in simulation calculations, saving computation time. Attached Figure Description

[0097] To more clearly illustrate the technical solutions in the specific embodiments of this application or the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0098] Figure 1 A flowchart of a two-phase flow simulation method for an aero-engine bearing cavity provided in this application embodiment;

[0099] Figure 2 A schematic diagram of the geometric model of the bearing cavity of an aero-engine provided in the embodiments of this application;

[0100] Figure 3 A schematic diagram of the Cartesian mesh used to divide the geometric model of the aero-engine bearing cavity provided in the embodiments of this application;

[0101] Figure 4 A manifold diagram of two-phase flow in the bearing cavity of an aero-engine provided for an embodiment of this application;

[0102] Figure 5 Functional structure diagram of the two-phase flow simulation device for aero-engine bearing cavity provided in the embodiments of this application;

[0103] Figure 6 This is a structural diagram of an electronic device provided in an embodiment of this application. Detailed Implementation

[0104] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0105] Therefore, the following detailed description of the embodiments of this application provided in the accompanying drawings is not intended to limit the scope of the claimed application, but merely to illustrate selected embodiments of the application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application.

[0106] First, a brief introduction to the design concept of the embodiments of this application will be given.

[0107] The simulation of two-phase flow in the bearing cavity of aero-engines involves deformation problems such as interface migration and fragmentation. Traditional computational fluid dynamics methods, such as VOF (Volume of Fluid) and Level Set, require tracking or capturing the interface, which increases the complexity of the computational model for two-phase flow problems and cannot consider the interactions between molecules from a microscopic perspective.

[0108] To this end, this application provides a simulation method for two-phase flow in the bearing cavity of an aero-engine. This method uses the Lattice Boltzmann Method (LBM) to calculate the two-phase flow in the bearing cavity of an aero-engine. As a diffusion interface method, it does not require tracking, capturing or reconstructing complex phase interfaces, thus reducing the computational complexity. At the same time, a refinement and coarsening level factor related to the phase volume fraction gradient is introduced during the calculation process. The mesh can be automatically refined and coarsened as the phase interface moves, ensuring accurate and fast calculation of the two-phase flow morphology.

[0109] This application can obtain highly accurate physical quantity results using simulation calculations. In the simulation calculations, the distribution pattern of lubricating oil in the engine bearing cavity can be accurately obtained with fewer grids and less computation, saving computation time and helping to explain the flow mechanism of flow phenomena in the aero-engine bearing cavity.

[0110] After introducing the application scenarios and design concepts of the embodiments of this application, the technical solutions provided by the embodiments of this application will be described below.

[0111] like Figure 1 As shown in the figure, this application provides a two-phase flow simulation method for an aero-engine bearing cavity, applied to the geometric model of the aero-engine bearing cavity. The bearing cavity model includes a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet, as shown in the figure. Figure 2As shown, the geometric model can be created using the 3D modeling software UG.

[0112] The method includes:

[0113] Step 101: Mesh the geometric model of the aero-engine bearing cavity to obtain a discretized Cartesian mesh;

[0114] To quickly and accurately calculate the shape changes of phase interfaces, a dynamic adaptive meshing method is used, which dynamically refines and coarsens the computational mesh based on changes in phase volume fraction. The main idea is to refine the mesh near the phase interface to improve the resolution of the interface, and coarse the mesh inside the gas or liquid phase to improve computational speed.

[0115] The specific mesh refinement method is as follows: calculate the gradient of the phase volume fraction, determine the refinement level factor, and dynamically refine the computational mesh. The gradient of the volume fraction and the refinement level factor satisfy the following relationship:

[0116]

[0117] Where A is the encryption level factor, and A0 is the encryption level factor constant, preferably set to 2. Let C be the magnitude of the phase volume fraction gradient, [] represents integer operations, tanh(·) is the hyperbolic tangent function, and C is the value of C. α The reference value for the magnitude of the gradient is set to 100, which is used to determine the degree of mesh refinement based on the magnitude of different gradients. A minimum mesh size limit is introduced here: when the size of the refined mesh reaches the minimum size limit, even if the refinement conditions are met, further mesh subdivision will not be performed. Using a dynamic adaptive mesh refinement method, the computational mesh at the phase interface can be automatically refined as the phase interface moves, improving computational accuracy. Mesh refinement is performed at the wall to improve computational accuracy near the wall, such as... Figure 3 As shown.

[0118] The specific mesh coarsening method involves calculating the gradient of the phase volume fraction and determining the coarsening level factor to dynamically coarse the computational mesh. The magnitude of the gradient of the volume fraction... The coarsening level factor B0 satisfies the following judgment relationship:

[0119]

[0120] Where ∈ represents a given coarsening threshold, with a value of 0.05. The idea is that when the magnitude of the gradient is less than the set coarsening threshold, two layers of mesh are merged; when the magnitude of the gradient is greater than the coarsening threshold, no mesh merging operation is performed. A maximum mesh size limit is introduced here: when the coarsened mesh size reaches the maximum size limit, even if the coarsening condition is met, no further mesh merging operations are performed. Using a dynamic adaptive mesh coarsening method, the computational mesh within the phase can be automatically coarsened as the phase interface moves, improving computational speed.

[0121] Step 102: Set the initial and boundary conditions for the computational domain of the discretized Cartesian mesh;

[0122] The initial conditions for two-phase simulation of the bearing cavity of an aero-engine are the initial distribution of lubricating oil and air in a given computational domain, specifically the volume fraction of the lubricating oil.

[0123] The boundary at the lubricating oil inlet can be configured as a macroscopic velocity inlet, flow inlet, or pressure inlet, depending on the specific physical scenario. For the flow inlet condition, since there is no direct flow quantity in the lattice Boltzmann method, it needs to be indirectly converted to an equivalent velocity inlet condition. This is based on the lubricating oil mass flow rate. The following relationship is satisfied by density ρ, velocity U, and inlet cross-sectional area A. The inlet velocity value is calculated to apply the boundary conditions.

[0124] Obtain the density and velocity at the lubricating oil inlet, and set the boundary conditions for the distribution function at the lubricating oil inlet:

[0125]

[0126] Where t represents the current time, x r ρ represents the coordinates at the boundary of the lubricating oil inlet. r u represents the density of the lubricating oil at the lubricating oil inlet boundary. r x represents the lubricating oil velocity vector at the lubricating oil inlet boundary. rf The coordinates of the fluid domain immediately adjacent to the lubricating oil inlet boundary, u rf ρ represents the velocity vector of the fluid domain immediately adjacent to the lubricating oil inlet boundary. rf f represents the density of the fluid domain immediately adjacent to the lubricating oil inlet boundary. r (x r (t) represents x at time t. r Density distribution function at that location, This represents the equilibrium state of the density distribution function.

[0127] The lubricating oil outlet of the bearing cavity can be designated as either a pressure outlet or a free outflow boundary, depending on the specific physical scenario. If the lubricating oil outlet of the bearing cavity is a pressure outlet, the specific implementation method for this pressure outlet is based on the state equation satisfied by incompressible fluids. The density is calculated given the pressure, where p is the pressure, ρ is the density, and c is the density. s The velocity of sound is used, and then the boundary conditions are implemented using the distribution function at the lubricating oil inlet. Specifically, after each iteration of the calculation, the distribution function value of the grid points at the lubricating oil outlet boundary is set to the distribution function value of the internal grid points immediately adjacent to the boundary.

[0128] The sealed air inlet is a velocity inlet, and it is treated in the same way as the lubricating oil velocity inlet. The boundary conditions are specifically implemented in the lattice Boltzmann method.

[0129] Obtain the density and velocity at the air inlet, and set the boundary conditions for the density distribution function at the air inlet:

[0130]

[0131] Where, x a ρ represents the coordinates at the air inlet boundary. a U represents the density at the air inlet boundary. a The velocity vector at the air inlet boundary, x af The coordinates of the fluid domain immediately adjacent to the air inlet boundary, u af ρ represents the velocity vector of the fluid domain immediately adjacent to the air inlet boundary. af f represents the density of the fluid domain immediately adjacent to the air inlet boundary. a (x r (t) represents x at time t. r The air density distribution function at that location, This represents the equilibrium state of the air density distribution function.

[0132] The air outlet of the bearing cavity can be designated as a pressure outlet or a free outflow boundary, depending on the specific physical scenario. The boundary conditions are implemented using the same treatment method as for the lubricating oil outlet, within the lattice Boltzmann method.

[0133] The rotating wall of the aero-engine bearing cavity is given a wall condition with a fixed rotation speed. The specific implementation method is as follows: the distribution function of the grid points at the rotating wall needs to be corrected for velocity after the collision and migration process, depending on whether momentum is gained or lost during the collision with the rotating wall.

[0134] The boundary condition for the density distribution function of the mixed fluid at the wall boundary is:

[0135]

[0136] Where, x b The coordinates at the wall boundary are represented by Δt, where Δt represents the time step, and u is the coordinates at the wall boundary. w This represents the velocity vector at the wall boundary. When the wall is a stationary wall of the cavity, u w =0, ρ w w represents the density at the wall boundary. k e represents the weighting coefficient. k Let c be the vector representing the direction of the k-th discrete velocity. s For the speed of sound, f k (x b (t) represents x at time t. b The density distribution function at point k, where k represents the index of the direction of the distribution function. It represents the distribution function of the velocity direction opposite to the k-th velocity direction after the collision at time t-Δt.

[0137] The stationary wall of the aero-engine bearing cavity is given as a non-slip wall condition. The specific implementation method is as follows: the LBM grid distribution function at the wall is bounced after the collision and migration process, and the value of the distribution function flowing into the computational domain from the wall boundary is taken as the value of the distribution function in the opposite direction.

[0138] Step 103: Based on the set initial and boundary conditions, perform flow field simulation calculations on the geometric model of the aero-engine bearing cavity to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment.

[0139] The specific solution process for the simulation calculation involves calculating the equations for the collision and migration processes according to a given calculation duration and a time-progression method until the set number of time steps is reached. The relaxation factor is τ = 0.8, the maximum time step is Δt = 0.0001s, and the kinematic viscosity of air is 1.48 × 10⁻⁶. -5 The kinematic viscosity of the lubricating oil is taken as 17.7390 m² / s. The local spatial step is determined based on the relationship between the relaxation factor, time step, spatial step, and kinematic viscosity.

[0140] In this embodiment, the step includes:

[0141] Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, calculate the vector f of the lubricating oil density distribution function at grid point x at the current time t. r , where f r =(f r，1 f r，2 , ...f r，KK is the number of discrete velocity directions;

[0142] Calculate the lubricating oil density ρ at grid point x at the current time t. r (t, x):

[0143]

[0144] Calculate the lubricating oil velocity u at grid point x at the current time t. r (t, x):

[0145]

[0146] Among them, e k It is the vector representing the direction of the k-th discrete velocity;

[0147] Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, the vector f of the air density distribution function at grid point x at the current time t is calculated. a , where f a =(f a，1 f a，2 , ...f a，K );

[0148] Calculate the air density ρ at grid point x at time t. a (t, x):

[0149]

[0150] Calculate the air velocity u at grid point x at the current time t. a (t, x):

[0151]

[0152] Then the mixed fluid density ρ(t, x) at grid point x at the current time t is:

[0153] ρ(t, x) = ρ r (t, x) + ρ a (t, x)

[0154] Then the velocity u(t, x) of the mixed fluid at grid point x at the current time t is:

[0155]

[0156] Among them, F r，total For the combined force of the lubricating oil, F a，total It is the resultant force of air.

[0157] The density and velocity distribution of the mixed fluid are imported into the post-processing software Tecplot to display the morphological distribution of the two-phase flow, such as... Figure 4 As shown.

[0158] Specifically, based on the evolution equation of the LBM distribution function of the multi-relaxation collision operator, the vector f of the lubricating oil density distribution function at grid point x at the current time t is calculated. r ;include:

[0159] Vector f r The k-th component f r，k The calculation formula is:

[0160]

[0161] Among them, f r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Density distribution function at Δt; M -1 S is the inverse of the transformation matrix M; r Let m be the relaxation matrix. r (xe k Δt, t-Δt) are the density distribution function f r (xe k The moments of Δt, t-Δt), Let f be the density distribution function r (xe k The equilibrium moment F of Δt, t-Δt) r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external force of the lubricating oil at point Δt;

[0162] Lubricating oil external force F r (xe k The formula for calculating Δt, t-Δt) is:

[0163]

[0164] Where, ρ r (t-Δt, x) represents the lubricating oil density at grid point x at the previous time t-Δt; u eq Let x be the equilibrium velocity of the mixed fluid at grid point x at the previous time t-Δt:

[0165]

[0166] Lubricating oil drift speed Δu r The calculation formula is:

[0167] Δu r =F r,total Δt / ρ r (t-Δt, x)

[0168] Among them, the combined force of lubricating oil F r,total for:

[0169]

[0170] The force F between lubricating oil and air ra for:

[0171]

[0172] Among them, g ra w represents the strength of the interaction between the lubricating oil and air. k The weighting coefficient for the k-th discrete velocity direction;

[0173] The force between the lubricating oil and the wall surface for:

[0174]

[0175] in, The strength of the interaction between the lubricating oil and the wall surface, s(xe) k Δt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0;

[0176] external forces of gravity for:

[0177]

[0178] Where g is the gravity coefficient.

[0179] Specifically, based on the LBM distribution function evolution equation of the multiple relaxation collision operator, the vector f of the air density distribution function at grid point x at the current time t is calculated. a ;include:

[0180] Vector f a The k-th component f a，k The calculation formula is:

[0181]

[0182] Among them, f a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Air density distribution function at Δt; S a Let m be the relaxation matrix.a (xe k Δt, t-Δt) are the density distribution function f a (xe k The moments of Δt, t-Δt), Let f be the density distribution function a (xe k The equilibrium moment F of Δt, t-Δt) a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external air force acting at point Δt;

[0183] External air force F a (xe k The formula for calculating Δt, t-Δt) is:

[0184]

[0185] Where, ρ a (t-Δt, x) represents the air density at grid point x at the previous time t-Δt.

[0186] air drift velocity Δu a The calculation formula is:

[0187] Δu a =F a，total Δt / ρ a (t-Δt, x)

[0188] Among them, the resultant force of air F a，total for:

[0189]

[0190] The force F between air and lubricating oil ar for:

[0191]

[0192] Among them, g ar The strength of the interaction between air and lubricating oil, w k The weighting coefficient for the k-th discrete velocity direction;

[0193] Forces between air and wall for:

[0194]

[0195] in, The intensity of the interaction between the air and the wall, s(xe kΔt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0;

[0196] external forces of gravity for:

[0197]

[0198] Based on the above embodiments, this application provides a two-phase flow simulation device for an aero-engine bearing cavity, see reference. Figure 5 As shown, the two-phase flow simulation device 200 for aero-engine bearing cavities provided in this application embodiment includes at least:

[0199] Mesh generation unit 201 is used to mesh the geometric model of the aero-engine bearing cavity to obtain a discretized Cartesian mesh;

[0200] Setting element 202 is used to set the initial conditions and boundary conditions of the computational domain of the discretized Cartesian mesh;

[0201] The calculation unit 203 is used to perform flow field simulation calculations on the geometric model of the aero-engine bearing cavity according to the set initial conditions and boundary conditions, and to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment.

[0202] It should be noted that the principle of the two-phase flow simulation device 200 for the aero-engine bearing cavity provided in this application embodiment to solve the technical problem is similar to the method provided in this application embodiment. Therefore, the implementation of the two-phase flow simulation device 200 for the aero-engine bearing cavity provided in this application embodiment can refer to the implementation of the method provided in this application embodiment, and the repeated parts will not be described again.

[0203] like Figure 6 As shown, the electronic device 300 provided in this application embodiment includes at least: a processor 301, a memory 302, and a computer program stored in the memory 302 and capable of running on the processor 301. When the processor 301 executes the computer program, it implements the two-phase flow simulation method for the bearing cavity of an aero-engine provided in this application embodiment.

[0204] The electronic device 300 provided in this application embodiment may further include a bus 303 connecting different components (including processor 301 and memory 302). The bus 303 represents one or more types of bus structures, including memory bus, peripheral bus, local area bus, etc.

[0205] The memory 302 may include a readable medium in the form of volatile memory, such as random access memory (RAM) 3021 and / or cache memory 3022, and may further include read-only memory (ROM) 3023.

[0206] The memory 302 may also include a program tool 3024 having a set (at least one) of program modules 3025, including but not limited to: an operating subsystem, one or more application programs, other program modules, and program data, each or some combination of these examples may include an implementation of a network environment.

[0207] Electronic device 300 can also communicate with one or more external devices 304 (e.g., keyboard, remote control, etc.), and with one or more devices that enable a user to interact with electronic device 300 (e.g., mobile phone, computer, etc.), and / or with any device that enables electronic device 300 to communicate with one or more other electronic devices 300 (e.g., router, modem, etc.). This communication can be performed through input / output (I / O) interface 305. Furthermore, electronic device 300 can also communicate with one or more networks (e.g., local area network (LAN), wide area network (WAN), and / or public networks, such as the Internet) through network adapter 306. Figure 6 As shown, network adapter 306 communicates with other modules of electronic device 300 via bus 303. It should be understood that, although... Figure 6 As not shown, other hardware and / or software modules may be used in conjunction with electronic device 300, including but not limited to: microcode, device drivers, redundant processors, external disk drive arrays, Redundant Arrays of Independent Disks (RAID) subsystems, tape drives, and data backup storage subsystems.

[0208] It should be noted that, Figure 6 The electronic device 300 shown is merely an example and should not impose any limitations on the functionality and scope of use of the embodiments of this application.

[0209] This application also provides a computer-readable storage medium storing computer instructions that, when executed by a processor, implement the two-phase flow simulation method for aero-engine bearing cavity provided in this application.

[0210] Furthermore, although the operations of the method of this application are described in a specific order in the accompanying drawings, this does not require or imply that these operations must be performed in that specific order, or that all the operations shown must be performed to achieve the desired result. Additionally or alternatively, certain steps may be omitted, multiple steps may be combined into one step, and / or one step may be broken down into multiple steps.

[0211] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0212] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

Claims

1. A two-phase flow simulation method for an aero-engine bearing cavity, applied to the geometric model of the aero-engine bearing cavity, wherein the geometric model includes a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet; characterized in that, The method includes: The geometric model of the aero-engine bearing cavity is meshed to obtain a discretized Cartesian mesh; Set the initial and boundary conditions for the computational domain of the discretized Cartesian mesh; Based on the set initial and boundary conditions, flow field simulation calculations are performed on the geometric model of the aero-engine bearing cavity to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment; Set the initial and boundary conditions for the computational domain of the discretized Cartesian mesh, including: Set the distribution of lubricating oil and air in the computational domain at the initial moment; Obtain the density and velocity at the lubricating oil inlet, and set the boundary conditions for the distribution function at the lubricating oil inlet: Where t represents the current time, x r ρ represents the coordinates at the boundary of the lubricating oil inlet. r u represents the density of the lubricating oil at the lubricating oil inlet boundary. r x represents the lubricating oil velocity vector at the lubricating oil inlet boundary. rf The coordinates of the fluid domain immediately adjacent to the lubricating oil inlet boundary, u rf ρ represents the velocity vector of the fluid domain immediately adjacent to the lubricating oil inlet boundary. rf f represents the density of the fluid domain immediately adjacent to the lubricating oil inlet boundary. r (x r (,t) represents x at time t r Density distribution function at that location, This represents the equilibrium state of the density distribution function; Obtain the density and pressure at the lubricating oil outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions of the distribution function at the lubricating oil outlet. Obtain the density and velocity at the air inlet, and set the boundary conditions for the density distribution function at the air inlet: Where, x a ρ represents the coordinates at the air inlet boundary. a U represents the density at the air inlet boundary. a The velocity vector at the air inlet boundary, x af The coordinates of the fluid domain immediately adjacent to the air inlet boundary, u af ρ represents the velocity vector of the fluid domain immediately adjacent to the air inlet boundary. af f represents the density of the fluid domain immediately adjacent to the air inlet boundary. a (x r (,t) represents x at time t r The air density distribution function at that location, This represents the equilibrium state of the air density distribution function; Obtain the density and pressure at the air outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions for the distribution function at the air outlet. The boundary condition for the density distribution function of the mixed fluid at the wall boundary is: Where, x b The coordinates at the wall boundary are represented by Δt, where Δt represents the time step, and u is the coordinates at the wall boundary. w This represents the velocity vector at the wall boundary. When the wall is a stationary wall of the cavity, u w =0, ρ w w represents the density at the wall boundary. k e represents the weighting coefficient. k Let c be the vector representing the direction of the k-th discrete velocity. s For the speed of sound, f k (x b (,t) represents x at time t b The density distribution function at point k, where k represents the index of the direction of the distribution function. It represents the distribution function of the velocity direction opposite to the k-th velocity direction after the collision at time t-Δt.

2. The method according to claim 1, characterized in that, Based on the set initial and boundary conditions, flow field simulation calculations are performed on the geometric model of the aero-engine bearing cavity to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment; including: Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, calculate the vector f of the lubricating oil density distribution function at grid point x at the current time t. r , where f r =(f r,1 ,f r,2 ,…f r,K K is the number of discrete velocity directions; Calculate the lubricating oil density ρ at grid point x at the current time t. r (t,x): Calculate the lubricating oil velocity u at grid point x at the current time t. r (t,x): Among them, e k It is the vector representing the direction of the k-th discrete velocity; Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, the vector f of the air density distribution function at grid point x at the current time t is calculated. a , where f a =(f a,1 ,f a,2 ,…f a,K ); Calculate the air density ρ at grid point x at time t. a (t,x): Calculate the air velocity u at grid point x at the current time t. a (t,x): Then the mixed fluid density ρ(t,x) at grid point x at the current time t is: p(t,x)=p r (t,x)+ρ a (t,x) Then the velocity u(t,x) of the mixed fluid at grid point x at the current time t is: Among them, F r,total For the combined force of the lubricating oil, F a,total It is the resultant force of air.

3. The method according to claim 2, characterized in that, Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, calculate the vector f of the lubricating oil density distribution function at grid point x at the current time t. r ; include: Vector f r The k-th component f r,k The calculation formula is: Among them, f r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Density distribution function at Δt; M -1 S is the inverse of the transformation matrix M; r Let m be the relaxation matrix. r (xe k Δt, t-Δt) is the density distribution function f r (xe k The moments of Δt, t-Δt), Let f be the density distribution function r (xe k The equilibrium moment F of Δt, t-Δt) r (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external force of the lubricating oil at point Δt; Lubricating oil external force F r (xe k The formula for calculating Δt, t-Δt) is: Where, ρ r (t-Δt,x) represents the lubricating oil density at grid point x at the previous time t-Δt; u eq Let x be the equilibrium velocity of the mixed fluid at grid point x at the previous time t-Δt: Lubricating oil drift speed Δu r The calculation formula is: Thu r =F r,total Δt / ρ r (t-Δt,x) Among them, the combined force of lubricating oil F r,total for: The force F between lubricating oil and air ra for: Among them, g ra w represents the strength of the interaction between the lubricating oil and air. k The weighting coefficient for the k-th discrete velocity direction; The force between the lubricating oil and the wall surface for: in, The strength of the interaction between the lubricating oil and the wall surface, s(xe) k Δt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0; external forces of gravity for: Where g is the gravity coefficient.

4. The method according to claim 3, characterized in that, Based on the LBM distribution function evolution equation of the multiple relaxation collision operator, calculate the vector f of the air density distribution function at grid point x at the current time t. a ; include: Vector f a The k-th component f a,k The calculation formula is: Among them, f a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k Air density distribution function at Δt; S a Let m be the relaxation matrix. a (xe k Δt, t-Δt) is the density distribution function f a (xe k The moments of Δt, t-Δt), Let f be the density distribution function a (xe k The equilibrium moment F of Δt, t-Δt) a (xe k Δt, t-Δt) represents the previous time t-Δt at xe k The external air force acting at point Δt; External air force F a (xe k The formula for calculating Δt, t-Δt) is: Where, ρ a (t-Δt,x) represents the air density at grid point x at the previous time t-Δt. air drift velocity Δu a The calculation formula is: Thu a =F a,total Δt / ρ a (t-Δt,x) Among them, the resultant force of air F a,total for: The force F between air and lubricating oil ar for: Among them, g ar The strength of the interaction between air and lubricating oil, w k The weighting coefficient for the k-th discrete velocity direction; Forces between air and wall for: in, The intensity of the interaction between the air and the wall, s(xe k Δt) is the density function near the wall, when the coordinates xe k When Δt is in the solid domain, s(xe) k Δt)=1, when coordinates xe k When Δt is in the fluid domain, s(x+e) k Δt)=0; external forces of gravity for:

5. The method according to claim 1, characterized in that, The method further includes: Based on the relationship between the relaxation factor τ and the time step Δt, spatial step Δx, and kinematic viscosity v: The range of the relaxation factor is 0.5 < τ < 2, which determines the values ​​of the step size Δt, spatial step size Δx, and kinematic viscosity v.

6. The method according to claim 1, characterized in that, After meshing the geometric model of the aero-engine bearing cavity to obtain the discretized Cartesian mesh, the following steps are also included: Calculate the magnitude of the phase volume fraction gradient at the phase interface Using the magnitude of gradient Determine the encryption level factor A: Where A0 is the encryption level factor constant, [] represents integer operation, tanh(·) is the hyperbolic tangent function, and C α The reference value for the magnitude of the gradient is set to determine the degree of encryption based on the magnitude of different gradients; For a square grid in the computational domain, perform 2 A This process involves multiple divisions to achieve mesh encryption at the interface.

7. The method according to claim 1, characterized in that, After meshing the geometric model of the aero-engine bearing cavity to obtain the discretized Cartesian mesh, the following steps are also included: The modulus for calculating the gradient of phase volume fraction within the gas or liquid phase. When the magnitude of the gradient When the value is less than the set coarsening threshold, the mesh is merged into two layers, thereby achieving mesh coarsening within the gas or liquid phase.

8. A two-phase flow simulation device for an aero-engine bearing cavity, applied to the geometric model of an aero-engine bearing cavity, the geometric model including a rotating cavity wall, a stationary cavity wall, a lubricating oil inlet, an air inlet, a lubricating oil outlet, and an air outlet; characterized in that, The device includes: Mesh generation unit is used to mesh the geometric model of the bearing cavity of an aero-engine to obtain a discretized Cartesian mesh; The setting unit is used to set the initial conditions and boundary conditions of the computational domain of the discretized Cartesian mesh; The computing unit is used to perform flow field simulation calculations on the geometric model of the aero-engine bearing cavity based on the set initial and boundary conditions, and to obtain the physical quantity distribution of the mixed flow field at each grid point at the current moment; The setting unit is specifically used for: Set the distribution of lubricating oil and air in the computational domain at the initial moment; Obtain the density and velocity at the lubricating oil inlet, and set the boundary conditions for the distribution function at the lubricating oil inlet: Where t represents the current time, x r ρ represents the coordinates at the boundary of the lubricating oil inlet. r u represents the density of the lubricating oil at the lubricating oil inlet boundary. r x represents the lubricating oil velocity vector at the lubricating oil inlet boundary. rf The coordinates of the fluid domain immediately adjacent to the lubricating oil inlet boundary, u rf ρ represents the velocity vector of the fluid domain immediately adjacent to the lubricating oil inlet boundary. rf f represents the density of the fluid domain immediately adjacent to the lubricating oil inlet boundary. r (x r (,t) represents x at time t r Density distribution function at that location, This represents the equilibrium state of the density distribution function; Obtain the density and pressure at the lubricating oil outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions of the distribution function at the lubricating oil outlet. Obtain the density and velocity at the air inlet, and set the boundary conditions for the density distribution function at the air inlet: Where, x a ρ represents the coordinates at the air inlet boundary. a U represents the density at the air inlet boundary. r The velocity vector at the air inlet boundary, x af The coordinates of the fluid domain immediately adjacent to the air inlet boundary, u af ρ represents the velocity vector of the fluid domain immediately adjacent to the air inlet boundary. af f represents the density of the fluid domain immediately adjacent to the air inlet boundary. a (x r (,t) represents x at time t r The air density distribution function at that location, This represents the equilibrium state of the air density distribution function; Obtain the density and pressure at the air outlet, calculate the density according to the state equation, and then calculate the corresponding velocity. Based on this, set the boundary conditions for the distribution function at the air outlet. The boundary condition for the density distribution function of the mixed fluid at the wall boundary is: Where, x b The coordinates at the wall boundary are represented by Δt, where Δt represents the time step, and u is the coordinates at the wall boundary. w This represents the velocity vector at the wall boundary. When the wall is a stationary wall of the cavity, u w =0, ρ w w represents the density at the wall boundary. k e represents the weighting coefficient. k Let c be the vector representing the direction of the k-th discrete velocity. s For the speed of sound, f k (x b (,t) represents x at time t b The density distribution function at point k, where k represents the index of the direction of the distribution function. It represents the distribution function of the velocity direction opposite to the k-th velocity direction after the collision at time t-Δt.

9. An electronic device, characterized in that, include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the method as claimed in any one of claims 1-7.

Images