Explicit volume sigma-Y atomization large vortex analog simulation method capable of being used for rotational flow cup atomization

Abstract

The invention discloses an explicit volume sigma-Y atomization large vortex simulation method for rotational flow cup atomization, which is used for calculation of primary atomization and secondary atomization under a coarse grid with the ratio of the average grid scale to the minimum droplet diameter being greater than 10, and the primary atomization and the secondary atomization contain multi-scale droplet flows. According to the method, a multi-scale droplet flow field containing primary atomization and secondary atomization at the same time can be calculated, a large vortex structure above a filtering scale can be captured by using LES + VOF under a coarse grid, and in addition, the influence of sub-grid scale pulsation surface tension on a mainstream can be calculated.

Description

Technical Field

[0001] The invention relates to the technical field of aero-engine atomization, and in particular to an explicit volume Σ-Y atomization large eddy simulation method that can be used for swirl cup atomization. Background Art

[0002] Atomization is a critical process in liquid combustion devices such as aircraft engines, and is related to factors such as fluid turbulence diffusion and the air resistance encountered by droplets when they pass through gas media. Atomization is an interfacial physical process, which is the result of the competition between inertia, viscosity and surface tension (or capillary) forces. The quality of atomization has an important influence on combustion performance and pollutant emissions. The droplet size distribution after atomization is wide, usually several orders of magnitude different from micron-sized droplets in the downstream spray to millimeter-sized large liquid clusters in the upstream continuous liquid core. The inherent complexity and multi-scale nature of atomization make both experimental research and numerical calculations very challenging.

[0003] At present, the quasi-direct numerical simulation of the atomization process widely uses gas-liquid interface capture methods, such as volume of fluid function (VOF), level set (LS) and diffusion interface (DI). These methods can automatically identify complex processes such as droplet deformation, breakup and aggregation. However, the LS method and some DI methods are not mass-conserving in nature and require additional numerical corrections. The VOF method itself satisfies mass conservation and is currently a hot topic of research.

[0004] The VOF method is an interface tracking method based on the Euler grid. In this method, the immiscible fluid components share a set of momentum equations, and the phase volume fraction α and the phase fraction transport equation are introduced to track the interface between phases in the computational domain. α It represents the ratio of the volume of the liquid phase to the volume of the grid. It contains a phase transport equation, which is a scalar transport equation: The concept of is based on the fact that two immiscible fluids can be modeled as a single fluid with singular values ​​of thermodynamic, transport and flow properties. The density ρ and dynamic viscosity μ at each point in space and time are: ρ = αρ h +(1-α)ρ l μ=αμ h +(1-α)μ l , the VOF equation is essentially the same as the single-phase flow equation, both of which contain a momentum equation and a continuity equation. The continuity equation is exactly the same. In the momentum equation, the VOF model adds a viscosity that varies with space, as well as an additional surface tension term. The continuous surface force model converts the influence of surface tension into a continuous volume in a region, σ is the tension coefficient, and κ is the curvature of the phase interface:

[0005] It has a wide range of applicability when dealing with macro-scale multiphase flow problems, and is suitable for multiphase flow problems in most industrial and engineering applications, such as wave motion, open channel flow, etc. In addition, the VOF method has been embedded in the commercial computational fluid dynamics software ANSYS Fluent and the open source computational fluid dynamics software OpenFOAM. It has been widely used and studied in academia, and has a relatively mature theoretical basis and engineering application experience.

[0006] A major drawback of the F method is that the numerical dissipation caused by the discretization of the phase fraction governing equations leads to excessive smearing of sharp interfaces at the boundaries, making it difficult to capture fine liquid structures using coarse grids. In order to reduce numerical diffusion and accurately capture the gas-liquid interface, very fine grids and high computational costs are usually required, so the traditional VOF method is difficult to apply to secondary atomization and subsequent fine droplet dispersion processes.

[0007] The VOF method is not suitable for large eddy simulation (LES) turbulence models. LES needs to analyze fluid structures such as large-scale vortices, and the grid required is finer than RANS. The phase fraction transport equation itself does not have a diffusion term. In the coarse-grid LES, the numerical dissipation of the discrete VOF equation has an enhanced impact on the flow field. The numerical dissipation grid is highly correlated, and the interface dynamics are highly sensitive to the grid, which will not only affect the small-scale pulsation quantity, but also further affect the large-scale statistics.

[0008] Large eddy simulation cannot correctly evaluate the subgrid stress caused by the presence of the interface in two-phase flow. Turbulence pulsation and interface inhomogeneity near the two-phase interface will produce subgrid stress, and the subgrid turbulence diffusion rate in the Smagorinsky turbulence model is not enough to simulate the subgrid pulsation caused by the inhomogeneous phase at the interface. In addition, due to the multi-scale characteristics of atomization, when the size of some droplets is smaller than the LES filter scale, the subgrid scale surface tension cannot be ignored, and the traditional VOF model cannot correctly analyze the two-phase flow.

[0009] In view of the above problems, the present invention provides an explicit volumetric ZY atomization large eddy simulation method that can be used for swirl cup atomization, which can use LES+VOF on a coarse grid to capture large eddy structures above the filter scale. Summary of the invention

[0010] The present invention provides an explicit volume Σ-Y atomization large eddy simulation method for swirl cup atomization, which can capture large eddy structures above the filter scale by using LES+VOF under a coarse grid.

[0011] The purpose of the present invention is to provide an explicit volume Σ-Y atomization large eddy simulation method for swirl cup atomization, which is used for the calculation of primary and secondary atomization under a coarser grid with a ratio of average grid size to minimum droplet diameter greater than 10. The primary and secondary atomization contain multi-scale droplet flows defined as:

[0012]

[0013] Furthermore, the specific steps of the method are as follows:

[0014] Step 1: Perform LES filtering on the VOF phase transport equation and introduce Farve averaging to obtain:

[0015]

[0016] in, is the volume average density, is the Favre average phase fraction, is the average speed of Favre, is the subgrid average phase fraction flux, α″ is the Favre average phase fraction pulsation, u″ is the Favre average velocity pulsation;

[0017] Step 2: Volume averaging will produce unenclosed subgrid α fluxes, by introducing an artificial diffusion coefficient D V Model it:

[0018]

[0019] Step 3: On the one hand, it measures the impact of sub-grid scale phase fraction fluctuations on large-scale flow, and on the other hand, it is equivalent to artificially adding diffusion terms. Appropriate artificial diffusion terms can suppress the impact of numerical dissipation caused by discreteness on the calculation results. V The calculation formula is as follows:

[0020]

[0021] Among them, l V is the grid-related characteristic length, taking the cube root of the grid volume, C αu is a constant, |Q / E| is the high vorticity characteristic factor;

[0022] Step 4: Limit D V It can only appear in unstable areas such as strong turbulence. Before the interface instability occurs, no interface curvature will be generated, and |Q / E| approaches 0, so it approaches 0. The formula is expressed as:

[0023]

[0024] in, is the strain rate tensor, is the vorticity tensor;

[0025] Step 5: Using the conservation of liquid mass, the phase fraction of the Farve average can be obtained Convert to

[0026]

[0027] Step 6: Apply LES filtering to the VOF momentum control equation and introduce Farve averaging to obtain: is the volume mean stress, is the sub-grid stress, is the volume average surface stress:

[0028]

[0029] Step 7: Stress of unclosed subgrid Use the following model:

[0030]

[0031] Among them, v V is the diffusion coefficient related to the inhomogeneity caused by the existence of the interface, and the calculation formula is:

[0032] v V =D V Sc V

[0033] Step 8: is the subgrid-scale stress mainly caused by isotropic turbulence, and its calculation formula is:

[0034]

[0035] Among them, the turbulent viscosity coefficient v t Using the conventional Smagorinsky model, C s is a dimensionless empirical coefficient. The turbulence intensity at the interface gradually decays to 0. The interface function Φ(α) is used to calculate v at the interface. t impose restrictions;

[0036] Step 9: When near the interface, the diffusion coefficient is mainly v V , and far away from the interface D V = 0, the diffusion coefficient degenerates to v t , based on the continuous surface tension model, the surface density is introduced to close the surface tension term, in which the sub-grid pulsation of the interface curvature is ignored. The surface density of each LES grid covers the contribution of the sub-grid pulsation to the surface tension. The formula is:

[0037]

[0038] in, is the volume average surface density, σ is the surface tension coefficient of the liquid;

[0039] Step 10: Establish the Farve average surface area per unit mass The transport equation, closed to the surface density, is:

[0040]

[0041] in, is the source term, which is modeled by the following formula:

[0042]

[0043] in, is a dimensionless number, τ V is the turbulence time scale;

[0044] Step 11. Σ * is the equilibrium value of local surface density generation and dissipation in the turbulent time scale, which is determined by the local minimum surface density Σ min Peace MeasureΣ equ Composition, formula:

[0045] ∑ * =∑ min +∑ equ

[0046]

[0047] in, is the turbulent kinetic energy, We * The critical Weber number is obtained from the subgrid pulsation velocity:

[0048]

[0049] Where F is defined as the local grid size l V The ratio of the characteristic scale to the volume of droplets contained in the local grid when the phase fraction is 5%;

[0050] Step 12: Since We in the surface density control equation * is a priori dimensionless parameter obtained when the average phase fraction is 5%. When the grid is relatively coarse, the phase fraction is too small and the source term cannot correctly analyze the change of surface density in the grid. The F factor can transform the scale to the priori We * The corresponding scale is:

[0051] 0.05 l 5a 3=αl V 3

[0052]

[0053] Step 13: The Sotel mean diameter (SMD) D can be calculated from the surface density 32 distributed:

[0054]

[0055] The present invention has the following advantages: the present invention adds a grid-related artificial diffusion term, the diffusion coefficient, to the phase volume fraction control equation, and by thickening the interface thickness, enables the fine droplet interface to be captured under a coarser grid, thereby correctly analyzing the interface dynamics of the two-phase flow field.

[0056] The present invention obtains the interface momentum diffusivity through artificial diffusion coefficient in the LES momentum control equation to model the sub-grid stress generated by the interface inhomogeneous phase.

[0057] The present invention models the surface tension term through the gas-liquid interface density, constructs an additional surface density transport equation to close it, and models the sub-grid scale surface tension pulsation through the surface density. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 This is a diagram of the present invention with an axial distance of 10 mm from the head;

[0059] Figure 2 This is a diagram of the present invention with an axial distance of 20 mm from the head;

[0060] Figure 3 This is the diagram of the present invention with an axial distance of 30 mm from the head;

[0061] Figure 4 This is a diagram of the present invention with an axial distance of 40 mm from the head;

[0062] Figure 5 This is a diagram of the present invention with an axial distance of 50 mm from the head;

[0063] Figure 6 It is a schematic diagram of a three-head model of the present invention;

[0064] Figure 7 It is a schematic cross-sectional view of the three heads of the present invention in the XZ plane, Y=0;

[0065] Figure 8 YZ plane of the present invention, X=0 Sauter mean diameter (D32, unit m) cloud diagram;

[0066] Fig. 9It is the XZ plane of the present invention, Y=0 Sauter mean diameter (D32, unit: m) cloud diagram. DETAILED DESCRIPTION

[0067] The present invention provides an explicit volume Σ-Y atomization large eddy simulation method for swirl cup atomization, which is used for the calculation of primary atomization and secondary atomization under a coarser grid with a ratio of average grid size to minimum droplet diameter greater than 10. The primary atomization and secondary atomization contain multi-scale droplet flows, defined as:

[0068]

[0069] Furthermore, the specific steps of the method are as follows:

[0070] Step 1: Perform LES filtering on the VOF phase transport equation and introduce Farve averaging to obtain:

[0071]

[0072] in, is the volume average density, is the Favre average phase fraction, is the average speed of Favre, is the subgrid average phase fraction flux, α″ is the Favre average phase fraction pulsation, u″ is the Favre average velocity pulsation;

[0073] Step 2: Volume averaging will produce unenclosed subgrid α fluxes, by introducing an artificial diffusion coefficient D V Model it:

[0074]

[0075] Step 3: On the one hand, it measures the impact of sub-grid scale phase fraction fluctuations on large-scale flow, and on the other hand, it is equivalent to artificially adding diffusion terms. Appropriate artificial diffusion terms can suppress the impact of numerical dissipation caused by discreteness on the calculation results. V The calculation formula is as follows:

[0076]

[0077] Among them, l V is the grid-related characteristic length, taking the cube root of the grid volume, C αu is a constant, |Q / E| is the high vorticity characteristic factor;

[0078] Step 4: Limit D V It can only appear in unstable areas such as strong turbulence. Before the interface instability occurs, no interface curvature will be generated, and |Q / E| approaches 0, so it approaches 0. The formula is expressed as:

[0079]

[0080] in, is the strain rate tensor, is the vorticity tensor;

[0081] Step 5: Using the conservation of liquid mass, the phase fraction of the Farve average can be obtained Convert to

[0082]

[0083] Step 6: Apply LES filtering to the VOF momentum control equation and introduce Farve averaging to obtain: is the volume mean stress, is the sub-grid stress, is the volume average surface stress:

[0084]

[0085] Step 7: Stress of unclosed subgrid Use the following model:

[0086]

[0087] Among them, ν V is the diffusion coefficient related to the inhomogeneity caused by the existence of the interface, and the calculation formula is:

[0088] ν V =D V Sc V

[0089] Step 8: is the subgrid-scale stress mainly caused by isotropic turbulence, and its calculation formula is:

[0090]

[0091] Among them, the turbulent viscosity coefficient ν t Using the conventional Smagorinsky model, C s is a dimensionless empirical coefficient. The turbulence intensity at the interface gradually decays to 0. The interface function Φ(α) is used to calculate v at the interface. t impose restrictions;

[0092] Step 9: When near the interface, the diffusion coefficient is mainly v V , and far away from the interface D V = 0, the diffusion coefficient degenerates to v t, based on the continuous surface tension model, the surface density is introduced to close the surface tension term, in which the sub-grid pulsation of the interface curvature is ignored. The surface density of each LES grid covers the contribution of the sub-grid pulsation to the surface tension. The formula is:

[0093]

[0094] in, is the volume average surface density, σ is the surface tension coefficient of the liquid;

[0095] Step 10: Establish the Farve average surface area per unit mass The transport equation, closed to the surface density, is:

[0096] in, is the source term, which is modeled by the following formula:

[0097]

[0098] in, is a dimensionless number, τ V is the turbulence time scale;

[0099] Step 11. Σ * is the equilibrium value of local surface density generation and dissipation in the turbulent time scale, which is determined by the local minimum surface density ∑ min Peace Measurement∑ equ Composition, formula:

[0100] ∑ * =∑ min +∑ equ

[0101]

[0102] in, is the turbulent kinetic energy, We * The critical Weber number is obtained from the subgrid pulsation velocity:

[0103]

[0104] Where F is defined as the local grid size l V The ratio of the characteristic scale to the volume of droplets contained in the local grid when the phase fraction is 5%;

[0105] Step 12: Since We in the surface density control equation * is a priori dimensionless parameter obtained when the average phase fraction is 5%. When the grid is relatively coarse, the phase fraction is too small and the source term cannot correctly analyze the change of surface density in the grid. The F factor can transform the scale to the priori We *The corresponding scale is:

[0106] 0.05 l 5a 3 =α l V 3

[0107]

[0108] Step 13: The Sotel mean diameter (SMD) D can be calculated from the surface density 32 distributed:

[0109]

[0110] In this embodiment, as shown in Table 1, there is a comparison between the conventional LES+VOF model and the multi-scale LES+VOF model proposed in this application, where Sc V is the interface-related Schmidt number, ∑ is the interface surface density, which is defined as the ratio of the gas-liquid interface surface area in a subgrid to the volume of the grid in which it is located.

[0111]

[0112]

[0113] In this embodiment, LES+VOF can be used to capture large eddy structures above the filter scale under coarser and finer grids, and the influence of sub-grid scale pulsating surface tension on the mainstream can also be calculated.

[0114] In this embodiment, the verification experiment uses the three-head aircraft engine combustion chamber air atomization model of the University of Connecticut, and the numerical calculation platform uses the open source computational fluid dynamics simulation platform OpenFOAM-5.0. Figure 1 , Figure 2 , Figure 3 , Figure 4 , Figure 5 As shown, it is a comparison between the calculation results based on the LES+ improved VOF model and the PDPA experimental measurement results, showing the distribution of the flow field SMD at different axial positions.

[0115] In this embodiment, Figure 6 , Figure 7 As shown, the example grid is 900W, the average grid size is about 0.3mm, the phase fraction in the core area of ​​the swirl zone is below 0.05%, the simulation results are at the same order of magnitude as the experimental results, and can capture the axial and radial transformation trends of the swirl droplet removal SMD. However, due to the influence of the artificial diffusion term, the absolute value of the radial velocity gradient is large, and some areas are over-smeared, while the axial change velocity is small, and the prediction accuracy is higher.

[0116] Although the specific embodiments of the present invention are described in detail in conjunction with the accompanying drawings, it should not be understood as limiting the scope of protection of the present invention. Within the scope described in the claims, various modifications and variations that can be made by those skilled in the art without creative work still belong to the scope of protection of the present invention.

Claims

1. An explicit volume Σ-Y atomization large eddy simulation method for swirl cup atomization, characterized by: A coarser grid with a ratio of average grid size to minimum droplet diameter greater than 10 is used for calculations of primary and secondary atomization containing multi-scale droplet flows, defined as: 。 2. The method for simulating explicit volume Σ-Y atomization large eddy simulation for swirl cup atomization according to claim 1, characterized in that: The specific steps of the method are as follows: Step 1: Perform LES filtering on the VOF phase transport equation and introduce Farve averaging to obtain: in, is the volume average density, is the Favre average phase fraction, is the average speed of Favre, is the subgrid average phase fraction flux, α″ is the Favre average phase fraction pulsation, u″ is the Favre average velocity pulsation; Step 2: Volume averaging will produce unenclosed subgrid α fluxes, by introducing an artificial diffusion coefficient D V Model it: Step 3: On the one hand, it measures the impact of sub-grid scale phase fraction fluctuations on large-scale flow, and on the other hand, it is equivalent to artificially adding diffusion terms. Appropriate artificial diffusion terms can suppress the impact of numerical dissipation caused by discreteness on the calculation results. V The calculation formula is as follows: Among them, l V is the grid-related characteristic length, taking the cube root of the grid volume, C αu is a constant, |Q / E| is the high vorticity characteristic factor; Step 4: Limit D V It can only appear in unstable areas such as strong turbulence. Before the interface instability occurs, no interface curvature will be generated, and |Q / E| approaches 0, so it approaches 0. The formula is expressed as: in, is the strain rate tensor, is the vorticity tensor; Step 5: Using the conservation of liquid mass, the phase fraction of the Farve average can be obtained Convert to Step 6: Apply LES filtering to the VOF momentum control equation and introduce Farve averaging to obtain: is the volume mean stress, is the sub-grid stress, is the volume average surface stress: Step 7: Stress of unclosed subgrid Use the following model: Among them, v V is the diffusion coefficient related to the inhomogeneity caused by the existence of the interface, and the calculation formula is: v V =D V Sc V Step 8: is the subgrid-scale stress mainly caused by isotropic turbulence, and its calculation formula is: Among them, the turbulent viscosity coefficient ν t Using the conventional Smagorinsky model, C s is a dimensionless empirical coefficient. The turbulence intensity at the interface gradually decays to 0. The interface function Φ(α) is used to calculate v at the interface. t impose restrictions; Step 9: When near the interface, the diffusion coefficient is mainly v V , and far away from the interface D V = 0, the diffusion coefficient degenerates to v t , based on the continuous surface tension model, the surface density is introduced to close the surface tension term, in which the sub-grid pulsation of the interface curvature is ignored. The surface density of each LES grid covers the contribution of the sub-grid pulsation to the surface tension. The formula is: in, is the volume average surface density, σ is the surface tension coefficient of the liquid; Step 10: Establish the Farve average surface area per unit mass The transport equation, closed to the surface density, is: in, is the source term, which is modeled by the following formula: in, is a dimensionless number, τ V is the turbulence time scale; Step 11, ∑ * is the equilibrium value of local surface density generation and dissipation in the turbulent time scale, which is determined by the local minimum surface density ∑ min Peace Measurement∑ equ Composition, formula: S * =S min +S equ in, is the turbulent kinetic energy, We * The critical Weber number is obtained from the subgrid pulsation velocity: Where F is defined as the local grid size l V The ratio of the characteristic scale to the volume of droplets contained in the local grid when the phase fraction is 5%; Step 12: Since We in the surface density control equation * is a priori dimensionless parameter obtained when the average phase fraction is 5%. When the grid is relatively coarse, the phase fraction is too small and the source term cannot correctly analyze the change of surface density in the grid. The F factor can transform the scale to the priori We * The corresponding scale is: 0.05 l 5a 3 ＝αl V 3 Step 13: The Sotel mean diameter (SMD) D can be calculated from the surface density 32 distributed:

Images