Unsteady numerical simulation method based on improved delay separation vortex simulation

Abstract

The embodiment of the invention provides an unsteady numerical simulation method and device based on improved delay separation vortex simulation and a storage medium, and the method comprises the steps: constructing a vortex inclination metric value function; turbulent kinetic energy and turbulence ratio dissipation rate are obtained in real time; on the basis of an RANS method, according to the turbulence energy and the turbulence ratio dissipation rate, the turbulence length scale is obtained; obtaining the maximum grid scale in the local vortex direction; on the basis of an LES method, according to the vortex inclination metric value function and the maximum grid scale, a shear layer self-adaptive grid filtering scale is obtained; obtaining a novel shielding function according to the standard shielding function, the additional shielding function and the suppression shielding function; and based on an improved DDES method, obtaining flow field data according to the shear layer adaptive grid filtering scale, the novel shielding function and the turbulence length scale. Based on the improved DDES method, the resolving power of the unsteady large-separation turbulence structure in the wide speed domain can be improved, and the calculation precision of wide-speed-domain large-separation unsteady flow simulation is improved.

Description

Technical Field

[0001] This document relates to the field of fluid mechanics technology, and in particular to a non-steady numerical simulation method, device and storage medium based on improved delayed detached eddy simulation. Background Art

[0002] In the field of aerospace, aircraft will encounter various complex unsteady flow phenomena when flying in a wide speed range. Through numerical simulation, the impact of these phenomena on aircraft performance can be predicted and evaluated, thereby guiding the optimal design of the aircraft and improving the stability and safety of the aircraft.

[0003] The existing technology mainly uses DES / DDES hybrid methods as the most mainstream RANS-LES hybrid method for numerical simulation. The RANS method is used to simulate the near-wall area where high-frequency small-scale motion dominates, and LES is used to calculate the unsteady separated flow area under low-frequency large-scale motion.

[0004] However, in the process of converting RANS to LES, the existing technology may encounter the problem that the solution of the flow field is neither RANS nor LES. This area is called the "gray zone". The induction of the gray zone will cause two types of problems: (1) modeled stress dissipation in the attached flow area; (2) the large eddy viscosity generated by the front RANS delays the development of mixing layer instability, resulting in a large deviation between the flow field solution and the actual situation.

[0005] Among them, the first type of gray zone induced problem is mainly caused by inappropriate grid scale in the boundary layer. When the grid is refined, the conversion from RANS to LES will produce gray zones in the attached flow boundary layer. At this time, if the grid density is not enough to support the resolution required for LES calculation, it will produce too small turbulent viscosity, which will cause MSD. The direct consequence of MSD is the forward movement of the separation line, leading to grid-induced separation (GIS). In order to avoid MSD and GIS as much as possible, the classic DES places high demands on the grid topology and distribution. The DDES method introduces the shielding function f d Redefining the length scale avoids the MSD problem to a great extent. The DDES method has achieved considerable success, but in-depth research shows that when the streamwise grid size is less than 0.3 times the local boundary layer thickness δ, the shielding function f d The shielding function gradually fails and cannot properly protect the attached flow boundary layer, and the friction resistance drops sharply. This characteristic of the shielding function makes mesh convergence very difficult. The study also found that when the flow has a strong adverse pressure gradient, the effect of the shielding function also decreases significantly. Many studies have adjusted the shielding function f dThe protection of the shielding function to the boundary layer can be improved by using the coefficients. However, the study also found that inappropriate coefficients can lead to over-protection, thereby delaying the development of mixing layer instability, which in turn strengthens the second type of problem induced by the gray zone. On the other hand, different coefficient values ​​may need to be calibrated for different flow phenomena and grids of different densities, which greatly increases the difficulty of applications for complex shapes.

[0006] The second type of gray zone induced problem is mainly solved by modifying the subgrid length scale Δ of the LES part in the hybrid method. max Adjusted to Δ based on the local vorticity weight ω , which makes the transition from RANS to LES in the free stream shear layer simulation faster, thus minimizing the delay in the construction of the free stream shear layer instability. However, this also makes it more difficult for the aforementioned shielding function to protect the attached flow boundary layer, making the first type of problem induced by the gray zone more prominent. Summary of the invention

[0007] In view of the above scheme, the present application aims to propose an unsteady numerical simulation method, device and storage medium based on improved delayed detached eddy simulation to solve at least one of the above technical problems.

[0008] In a first aspect, one or more embodiments of this specification provide an unsteady numerical simulation method based on improved delayed detached eddy simulation, comprising:

[0009] Constructing vortex tilt measurement function;

[0010] Obtain turbulent kinetic energy and turbulent specific dissipation rate in real time;

[0011] Based on the RANS method, a turbulence length scale is obtained according to the turbulent kinetic energy and the turbulence specific dissipation rate;

[0012] The maximum grid size for obtaining the local vortex direction;

[0013] Based on the LES method, a shear layer adaptive grid filter scale is obtained according to the vortex tilt metric function and the maximum grid scale;

[0014] According to the standard shielding function, the additional shielding function and the suppressed shielding function, a new shielding function is obtained; and

[0015] Based on the improved DDES method, flow field data are obtained according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

[0016] Furthermore, using a dissipation detector, the flow field is determined to be a smooth region or a discontinuous region;

[0017] If the flow field is a smooth region, the symmetric reconstruction variables are used;

[0018] If the flow field is a discontinuous region, use monotonic reconstruction variables; and

[0019] A mixed reconstruction variable is obtained according to the symmetric reconstruction variable and the monotonic reconstruction variable.

[0020] Furthermore, the hybrid reconstruction variable calculation formula is as follows:

[0021]

[0022] in, represents a mixed reconstruction variable;

[0023] q L,R represents the traditional monotone reconstruction variable;

[0024] ψ represents the dissipative control function; and

[0025] represents a symmetric reconstruction variable.

[0026] Further, the velocity trace, vortex vector, vorticity magnitude, kinematic viscosity coefficient and turbulent kinematic viscosity coefficient are obtained;

[0027] Constructing a vortex inclination measurement function according to the velocity trace, the vortex vector, the vorticity magnitude, the kinematic viscosity coefficient and the turbulent kinematic viscosity coefficient;

[0028] According to the vortex inclination measurement function, an instability coefficient function is obtained; and

[0029] The shear layer adaptive grid filter scale is determined according to the instability coefficient function and the maximum grid scale.

[0030] Furthermore, the shear layer adaptive grid filter scale calculation formula is as follows:

[0031]

[0032] Among them, Δ SLA represents the scale of the shear layer adaptive grid filter;

[0033] represents the maximum grid size based on the local vortex direction;

[0034] F KH represents the Kelvin-Helmholtz instability coefficient function;

[0035] i and j represent the numbers of units i and j;

[0036] nb(i) represents the neighbor unit number of unit i; and

[0037] I ij Represents the length vector of the line connecting the centers of cells i and j.

[0038] Furthermore, the data measured by the outer boundary layer detector are used to construct an additional shielding function;

[0039] Using the data measured by the shear layer detector, a suppression function is constructed;

[0040] Using data measured by the basic detector, a standard shielding function is constructed; and

[0041] A new shielding function is obtained according to the additional shielding function, the suppression function and the standard shielding function.

[0042] Furthermore, the calculation formula of the new shielding function is as follows:

[0043] f P =f d ·(1-(1-f P2 )·f R )

[0044] Among them, f P Represents a new type of shielding function;

[0045] f d Represents the standard masking function;

[0046] f P2 represents an additional masking function; and

[0047] f R Represents the suppression function.

[0048] Furthermore, the calculation formula of the improved DDES method is as follows:

[0049] l DDES = l RANS -f p max{0,l RANS -C DES Δ SLA}

[0050] Among them, l DDES represents the turbulence length scale of the DDES method;

[0051] l RANS represents the turbulence length scale of the RANS method;

[0052] f p Represents a new type of shielding function;

[0053] CDES is the calibration factor for the DDES method; and

[0054] Δ SLA Represents the scale of the shear layer adaptive grid filter.

[0055] In a second aspect, an embodiment of the present application provides an unsteady numerical simulation device based on improved delayed detached eddy simulation, comprising:

[0056] Metric module, used to construct vortex tilt metric function;

[0057] An acquisition module is used to obtain turbulent kinetic energy and turbulent specific dissipation rate in real time;

[0058] A RANS scale module, for obtaining a turbulence length scale based on a RANS method according to the turbulent kinetic energy and the turbulence specific dissipation rate;

[0059] The grid scale module is used to obtain the maximum grid scale of the local vortex direction;

[0060] A filter scale module, used for obtaining a shear layer adaptive grid filter scale based on the LES method according to the vortex tilt metric function and the maximum grid scale;

[0061] A shielding module, used for obtaining a new shielding function according to a standard shielding function, an additional shielding function and a suppressed shielding function; and

[0062] The improved calculation module is used to obtain flow field data based on the improved DDES method according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

[0063] In a third aspect, an embodiment of the present application provides a storage medium for storing computer executable instructions, characterized in that when the computer executable instructions are executed, the steps of the unsteady numerical simulation method based on improved delayed detached eddy simulation described in any one of the first aspects are implemented.

[0064] Compared with the prior art, this application can at least achieve the following technical effects:

[0065] The present application can utilize the vortex inclination measurement function to promote the development of Kelvin-Helmholtz instability in the early stage of shear layer development; then strengthen the protection of the outer layer of the boundary layer by adding a shielding function; and avoid the activation of the additional shielding function in the shear layer and the generation of unnecessary shielding effects by suppressing the shielding function; thereby enabling the improved DDES method to improve the two types of contradictory gray zone problems faced by existing DDES-type hybrid methods, while strengthening the protection of the attached flow boundary layer, promoting the instability of the shear layer, thereby effectively improving the computational accuracy of unsteady flow simulation with a wide speed domain and large separation. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] In order to more clearly illustrate one or more embodiments of this specification or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this specification. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0067] Figure 1 A flow chart of an unsteady numerical simulation method based on improved delayed detached eddy simulation provided for one or more embodiments of this specification;

[0068] Figure 2 A schematic diagram of a subsonic rear step flow structure and simulation conditions provided for one or more embodiments of this specification;

[0069] Figure 3 A schematic diagram of a two-dimensional cross section of a reference grid provided for one or more embodiments of this specification;

[0070] Figure 4 A schematic cross-sectional diagram of flow field information extraction provided by one or more embodiments of this specification;

[0071] Figure 5 A comparative schematic diagram of average flow velocity profiles at different x-sections provided for one or more embodiments of the present specification;

[0072] Figure 6 A schematic diagram showing a comparison of average shear stress u'u' at different x-sections provided by one or more embodiments of this specification;

[0073] Figure 7 A schematic diagram showing a comparison of average shear stress u'v' at different x-sections provided by one or more embodiments of this specification;

[0074] Figure 8 A schematic diagram of the symmetric surface grid distribution provided for one or more embodiments of this specification;

[0075] Fig. 9 A schematic diagram of the grid distribution of the bottom cross section of a cylinder provided for one or more embodiments of this specification;

[0076] Fig.10 A schematic diagram of the velocity profile 1 mm before the bottom of a cylinder provided for one or more embodiments of this specification;

[0077] Fig.11 A schematic diagram of a standard DDES provided for one or more embodiments of this specification;

[0078] Fig.12A schematic diagram of an improved DDES provided for one or more embodiments of this specification;

[0079] Fig.13 A schematic diagram of the flow cross-sectional position provided for one or more embodiments of this specification;

[0080] Fig.14 A schematic diagram of the time-averaged flow velocity of different flow sections under a dense grid provided in one or more embodiments of this specification;

[0081] Fig.15 Schematic diagram of the time-averaged radial velocity of different flow direction sections under a dense grid provided for one or more embodiments of this specification

[0082] Fig.16 A schematic diagram of turbulent kinetic energy of different flow direction cross sections under a dense grid provided for one or more embodiments of this specification;

[0083] Fig.17 A schematic diagram of turbulent shear stress of different flow direction cross sections under a dense grid provided for one or more embodiments of this specification;

[0084] Fig.18 A schematic diagram of the structure of an unsteady numerical simulation device based on improved delayed detached eddy simulation provided in one or more embodiments of this specification. DETAILED DESCRIPTION

[0085] In order to enable those skilled in the art to better understand the technical solutions in one or more embodiments of this specification, the following will be combined with the drawings in one or more embodiments of this specification to clearly and completely describe the technical solutions in one or more embodiments of this specification. Obviously, the described embodiments are only part of the embodiments of this specification, not all of the embodiments. Based on one or more embodiments of this specification, all other embodiments obtained by ordinary technicians in this field without creative work should fall within the scope of protection of this document.

[0086] At present, the existing turbulence simulation technology mainly uses the RANS-LES hybrid method for calculation. The RANS-LES hybrid method combines the advantages of both RANS and LES. The RANS method is used for simulation in the near-wall area to reduce the calculation cost; the LES method is used for simulation in the area far from the wall to capture the detailed characteristics of turbulence. However, the existing technology will produce gray areas in the attached flow boundary layer during the conversion process from RANS (Reynolds-Averaged Navier-Stokes) to LES (Large Eddy Simulation), making the solution of the flow field neither a pure RANS solution nor a LES solution, resulting in low accuracy in flow field calculation.

[0087] In view of the above technical problems, this application proposes an unsteady numerical simulation method based on delayed separation simulation, such as Figure 1 As shown, the specific steps are as follows:

[0088] Step S1, constructing a vortex inclination measurement function.

[0089] In an embodiment of the present application, the trace of velocity, the vortex vector, the vorticity size, the kinematic viscosity coefficient and the turbulent kinematic viscosity coefficient are obtained; and a vortex inclination measurement function is constructed based on the trace of velocity, the vortex vector, the vorticity size, the kinematic viscosity coefficient and the turbulent kinematic viscosity coefficient.

[0090] Specifically, the calculation formula of the vortex tilt metric function is as follows:

[0091]

[0092] Wherein, VTM represents the vortex tilt measurement function;

[0093] ω represents the vortex vector;

[0094] Ω represents the vorticity size;

[0095] A trace that indicates speed;

[0096] ν represents the kinematic viscosity coefficient;

[0097] ν t represents the turbulent viscosity coefficient.

[0098] Step S2, obtaining turbulent kinetic energy and turbulent specific dissipation rate in real time.

[0099] In the embodiment of the present application, when performing numerical simulation calculations, the turbulent kinetic energy and turbulent specific dissipation rate can be obtained in real time by solving turbulence models (such as k-omega model, LES, etc.).

[0100] For example, in the RANS method, the k-omega model (where k represents the turbulent kinetic energy and omega represents the dissipation rate of the turbulent kinetic energy) can be used to calculate the turbulent kinetic energy. This model closes the time-averaged governing equations of turbulence by solving two differential equations to obtain the distribution of the turbulent kinetic energy k.

[0101] Step S3, based on the RANS method, obtaining the turbulence length scale according to the turbulent kinetic energy and the turbulence specific dissipation rate.

[0102] In the embodiment of the present application, the distribution of turbulent kinetic energy and the turbulent specific dissipation rate are determined by the turbulent kinetic energy to determine the length scale of the turbulence.

[0103] The turbulence length scale calculation formula of the RANS method is:

[0104]

[0105] Where k is the turbulent kinetic energy;

[0106] omega is the turbulence specific dissipation rate;

[0107] Coefficient β * =0.09.

[0108] Step S4, obtaining the maximum grid size in the local vortex direction.

[0109] In an embodiment of the present application, the large-scale vortex structure is captured to obtain the maximum grid scale in the local vortex direction. The grid scale must be smaller than the maximum vortex scale in the flow. If the grid scale is too large, the large-scale vortex structure cannot be captured, thereby failing to accurately describe the flow characteristics.

[0110] Step S5, based on the LES method, according to the vortex tilt metric function and the maximum grid scale, obtain the shear layer adaptive grid filter scale.

[0111] In an embodiment of the present application, in a large eddy simulation (LES), a vortex tilt metric function is calculated in the entire computational domain to obtain a vortex tilt metric value, and the key turbulence characteristics in the shear layer are captured by the vortex tilt metric function. According to the size and distribution of the vortex tilt metric value, a criterion for adaptive mesh encryption is formulated, which can better resolve the turbulence structure in the shear layer. Based on the adaptive mesh encryption criterion and the maximum mesh scale, the adaptive mesh filter scale of each mesh unit is calculated.

[0112] Shear layer adaptive grid filter scale Δ SLA The calculation formula is as follows:

[0113]

[0114] in, represents the maximum grid size based on the local vortex direction;

[0115] The subscripts i and j indicate the numbers of the units i and j;

[0116] F KH represents the Kelvin-Helmholtz instability coefficient function;

[0117]

[0118] α2 represents the upper threshold of the Kelvin-Helmholtz instability coefficient function;

[0119] α1 represents the lower threshold of the Kelvin-Helmholtz instability coefficient function;

[0120] Represents a constant coefficient, used to limit F KH The scope of the function;

[0121] VTM represents the vortex tilt measurement function;

[0122]

[0123] ω represents the vortex vector;

[0124] Ω represents the vorticity size;

[0125] A trace that indicates speed;

[0126] ν represents the kinematic viscosity coefficient;

[0127] ν t represents the turbulent motion viscosity coefficient;

[0128] The subscripts i and j indicate the numbers of the units i and j;

[0129] nb(i) represents the neighbor unit number of unit i;

[0130] I ij represents the length vector of the line connecting the centers of cells i and j;

[0131] I ij =n ω ×(r i -r j )

[0132] n ω represents the dimensionless vortex direction vector;

[0133] r i represents the position vector of unit i;

[0134] r j represents the position vector of cell j.

[0135] In this application, the rapid destabilization of the shear layer in the separation region can be effectively promoted by considering the grid filter scale affected by the flow vortex structure.

[0136] Step S6, obtaining a new shielding function according to the standard shielding function, the additional shielding function and the suppressed shielding function.

[0137] In an embodiment of the present application, an additional shielding function is constructed using data measured by an outer boundary layer detector; a suppression function is constructed using data measured by a shear layer detector; a standard shielding function is constructed using data measured by a basic detector; and a new shielding function is obtained based on the additional shielding function, the suppression function and the standard shielding function.

[0138] Specifically, an additional shielding function is constructed through the turbulent viscosity gradient detector to strengthen the protection of the outer layer of the boundary layer; a suppression function is constructed through the vorticity normal gradient detector to avoid the additional shielding function being activated in the shear layer and producing unnecessary shielding effects; the detector detects the data of the inner layer of the boundary layer and constructs a standard shielding function. In the standard shielding function f d Based on the coupling, an additional shielding function f p2 and the suppression function f R , and obtain a new shielding function f p .

[0139] The calculation formula of the new shielding function is as follows:

[0140] f P =f d ·(1-(1-f P2 )·f R )

[0141] Among them, f p Represents a new type of shielding function;

[0142] f d Represents the standard masking function;

[0143]

[0144] represents a constant coefficient, with a value of 8;

[0145] r d Indicates a defined detection function;

[0146]

[0147] v represents the molecular viscosity coefficient;

[0148] v t represents the dynamic eddy viscosity coefficient;

[0149] u i,j represents the velocity gradient;

[0150] k represents the Karman constant;

[0151] d represents the distance to the wall;

[0152] f p2 Represents an additional shielding function;

[0153] represents the outer boundary layer probe;

[0154] c3 represents the calibration factor of the outer boundary layer detector;

[0155] represents the gradient of turbulent viscosity in the normal direction of the wall;

[0156] f R represents the suppression function;

[0157]

[0158] represents the shear layer detector;

[0159] c4 represents the shear layer detector threshold;

[0160] η represents the intermediate function formed by α;

[0161] α represents Intermediate variable combined with c4.

[0162] In the present application, the additional shielding function is used to effectively protect the outer layer of the boundary layer, while the suppression function is used to prevent the additional shielding function from being opened in the free shear flow.

[0163] Step S7, based on the improved DDES (Delayed DES) method, the flow field data is obtained according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

[0164] In the embodiment of the present application, DDES is used as a hybrid method, which automatically adjusts between the RANS method and the LES method and converts to the corresponding method as needed; when the new shielding function f p When =0, the RANS method is selected and the length scale becomes the RANS length scale; otherwise, the LES method is selected and the length scale becomes the LES length scale, thereby obtaining the required flow field data.

[0165] The calculation formula of the improved DDES method is as follows:

[0166] l DDES = l RANS -f p max{0,l RANS -C DES Δ SLA}

[0167] Among them, l DDES Represents the turbulence length scale of the DDES method.

[0168] l RANS Represents the turbulence length scale for the RANS method.

[0169] f p Represents a new type of shielding function;

[0170] C DES is the calibration coefficient of the DDES method;

[0171] Δ SLA Represents the scale of the shear layer adaptive grid filter.

[0172] Preferably, in LES, a dissipation detector is used to determine whether the flow field is a smooth area or a discontinuous area; if the flow field is a smooth area, symmetric reconstruction variables are used; if the flow field is a discontinuous area, monotonic reconstruction variables are used; and based on the symmetric reconstruction variables and the monotonic reconstruction variables, a mixed reconstruction variable is obtained.

[0173] The smooth area of ​​the flow field requires low dissipation and requires the use of symmetric reconstruction variables To maintain computational stability in the discontinuous region of the flow field, it is necessary to use the traditional monotonic reconstruction variable q L,R ; The above symmetric reconstruction variables and monotonic reconstruction variables are mixed through a dissipation detector Ψ (the value of Ψ is between 0 and 1).

[0174] The formula for calculating the hybrid reconstruction variable is as follows:

[0175]

[0176] in, represents a mixed reconstruction variable;

[0177] q L,R represents the traditional monotone reconstruction variable;

[0178]

[0179] q L represents the left side value of the face obtained by monotonic reconstruction;

[0180] represents the limiter of the left unit;

[0181] q R represents the right side value of the face obtained by monotonic reconstruction;

[0182] represents the limiter of the right unit;

[0183] q i The variable q representing unit i;

[0184] q k The variable q representing unit k;

[0185] represents a symmetric reconstruction variable;

[0186]

[0187] β represents the coefficient of the symmetric reconstruction format, usually taking the value of 2 / 3;

[0188] Represents the length vector from the face center to the left unit center;

[0189] represents the gradient of a variable;

[0190] Represents the length vector from the face center to the right unit center;

[0191] ψ represents the dissipative control function;

[0192]

[0193] represents the Hamiltonian operator;

[0194] u represents the velocity vector;

[0195] ε represents a small positive value set to avoid the denominator being 0;

[0196] In this application, when Ψ=1, it indicates that the region is a discontinuous region, and the mixed variable Degenerates into a monotone reconstruction variable q L,R , thus ensuring the stability of calculation.

[0197] When Ψ=0, it indicates that the region is a smooth region, and the mixed variable Degenerates into symmetric reconstruction variables This ensures low dissipation.

[0198] In this application, a dissipative detector is used to distinguish the smooth area and the discontinuous area of ​​the LES part. In the smooth area, the flow is relatively stable, and a lower resolution and less computing resources can be used. In the discontinuous area, the flow changes violently, requiring a higher resolution and more computing resources. The dissipative detector can dynamically adjust the computing resources according to the regional characteristics, thereby improving the computing efficiency and improving the resolution of unsteady large separation turbulent structures.

[0199] In the embodiments of the present application, the improved DDES numerical simulation method is verified by the following calculation model. Compared with the existing DDES method, the calculation results are more accurate.

[0200] 1. Subsonic back-step separation flow simulation

[0201] At H = 1.27 cm, M ∞ =0.128, W=12H, δ=1.9cm, U ∞ =44.2m / s, Re H = 37000, the flow structure of the subsonic back step research example is simulated, such as Figure 2 As shown; where H is the step height; M ∞ is the incoming flow Mach number; W is the step width, which indicates that the step width is 12 times the step height; δ is the incoming flow boundary layer thickness; U ∞ is the incoming flow velocity; Re H is the Reynolds number based on the step height; the physical time step is taken as 6.25e-6s to capture the unsteady effects of the flow. Given a two-dimensional cross section of the base grid, such as Figure 3 As shown in the figure, the span width z / H=4, periodic boundary conditions are used, and 100 points are evenly arranged so that the grid spacing Δx≈Δz is used to capture the three-dimensional unsteady flow. The grids in the shear layer and flow separation zone are also encrypted. The total number of grid elements is about 7 million. Figure 4 The x-section shown extracts the corresponding flow field information, including the mean velocity pattern and shear stress. Figure 5 The distribution of the flow velocity u at different x-section positions in the time-averaged flow field is given, and the results of the RANS model are also given for comparison. It can be seen from the figure that at x / H = -4 in front of the step, the velocity pattern of the incoming flow is in good agreement with the experiment, ensuring the correctness of the inlet conditions. Among the flow velocity patterns at different x-sections at the rear of the step, the improved DDES method has the best result, and it is in good agreement with the experiment at each x-section. Figure 5 It can also be seen that the results of the standard DDES method at x / H=1 and x / H=6 are slightly worse; while the underlying velocity obtained by the pure RANS method is smaller, which is particularly obvious at x / H=6 and x / H=10.

[0202] pass Figure 6and Figure 7 The shear stress u'u' and u'v' distributions at different x-section positions in the time-averaged flow field are given. It can be seen from the figure again that the results of the improved DDES method are significantly better than those of the standard DDES method. The shear stress peaks and shapes of each section of the improved DDES method are in better agreement with the experiment; while the standard DDES method significantly underestimates the shear stress at x / H=1 due to the hysteresis of shear layer instability, and overestimates the shear stress at the downstream x / H=6.

[0203] Finally, the unsteady simulation of subsonic back-step separation flow shows that the improved DDES method using a new shielding function and shear layer adaptive grid filtering can give results that are more consistent with the experiment in terms of the statistics of time-averaged and pulsating quantities, and is superior to the standard DDES method using a standard shielding function and maximum grid size filtering. Compared with the technical indicators, the test results are as follows: (1) The test Mach number is within the indicator range; (2) Both the improved DDES method and the standard DDES method accurately obtain the flow reattachment position; (3) The improved DDES method is superior to the standard DDES method in obtaining flow field information such as average velocity type and turbulent pulsating quantity.

[0204] 2. Simulation of supersonic bottom separation flow of a cylinder

[0205] The flow simulation condition is the incoming flow Mach number M ∞ =2.46, corresponding to U ∞ =593.8m / s, the incoming flow temperature and pressure are T ∞ =145K, P ∞ =31415Pa. The bottom radius R of the cylinder is 31.75mm. The Reynolds number Re based on the bottom radius R is R =1.42875e+6. The calculation used two sets of grids, sparse and dense, to analyze the influence of the grid. The total number of cells in the dense grid is about 9.8 million; the total number of cells in the sparse grid is about 4.7 million. The flow distribution of the two sets of grids is basically the same, and the number of circumferential points of the sparse grid is about half of that of the dense grid. Figure 8 The distribution of the symmetry surface mesh is given; Fig. 9 The spatial distribution of two sets of grids at the bottom section of the cylinder is given, which are dense grids and sparse grids. The RANS method uses the SST turbulence model, and the RANS-LES hybrid method uses the SST-DDES method based on the SST turbulence model. For the unsteady DDES simulation, the low dissipation numerical format developed in this project is used, the physical time step is 1.0e-6s, and the statistical time for time-averaged flow field analysis is 0.04s. Fig.10 The comparison between the experimental results and the calculation of the velocity type in the first 1mm of the bottom is given. As can be seen from the figure, the inlet velocity type obtained by different grids and different methods is in good agreement with the experiment, ensuring that the initial conditions of the separated flow at the bottom are consistent with the experiment. Fig.11 (Standard DDES) and Fig.12 (Improved DDES) It can be seen that the improved DDES method can obtain richer flow field structures than the standard DDES method, which is mainly due to the improvement of the grid filtering method. The improved DDES method obviously promotes shear layer instability faster. The time-averaged field and pulsating field are further analyzed based on the calculation results. Fig.13 Four flow section locations are given for analysis. Fig.14 and Fig.15 The time-averaged streamwise velocity and radial velocity distribution at different streamwise cross-sectional positions of the dense grid are given. As can be seen from the figure, the improved DDES method simulates the velocity distribution in the shear layer quite well, significantly better than the standard DDES method; on the other hand, the SST turbulence model can also simulate the average velocity field well.

[0206] from Fig.14 (a) It can be seen that the streamwise velocity predicted by the improved DDES method is slightly larger when the flow just begins to separate. This is because the local angular vortex generated by the improved DDES method has a larger range. At this time, the flow has just left the object surface, and the current grid scale may not meet the simulation requirements for the initial development of the shear layer. Fig.14 It can be seen again that the recirculation velocity predicted by the improved DDES method at the centerline is too large. Some literature points out that the simulation of the recirculation zone at the bottom of the supersonic flow may require a higher-order format or a denser grid.

[0207] For turbulent pulsation, Fig.16 and Fig.17 The distribution of turbulent kinetic energy k and turbulent shear stress at different flow section positions under dense grids is given. As can be seen from the figure, the standard DDES method cannot predict the pulsation amount well at each section position, the value is small, and the maximum pulsation position is also far from the experiment. The improved DDES method significantly improves the calculation of unsteady pulsation; for the x / R=0.1575 section, the previous analysis pointed out that the grid required for the flow characteristic simulation may not be enough, and the pulsation calculated by the improved DDES method is slightly smaller; for the x / R=0.9449 section, because the predicted backflow velocity is too large, the pulsation in the central area is also larger.

[0208] Finally, the unsteady simulation of cylindrical supersonic bottom separation flow shows that: (1) Compared with the standard DDES method, the improved DDES method significantly improves the simulation accuracy of cylindrical supersonic bottom separation flow by better simulating the shear layer instability of the flow, and obtains the position of the stagnation point after the flow more accurately; (2) Whether it is time-averaged or pulsating quantities, the results of the improved DDES method are better than those of the standard DDES method and are in better agreement with the experiment; (3) The improved DDES method is less affected by the circumferential grid size and is more advantageous for simulating flows with complex shapes than the standard DDES method.

[0209] The present application embodiment provides an unsteady numerical simulation device based on improved delayed detached eddy simulation, such as Fig.18 As shown, including:

[0210] A measurement module 101 is used to construct a vortex inclination measurement value function;

[0211] An acquisition module 102 is used to acquire turbulent kinetic energy and turbulent specific dissipation rate in real time;

[0212] A RANS scale module 103, used for obtaining a turbulence length scale based on a RANS method according to the turbulent kinetic energy and the turbulence specific dissipation rate;

[0213] A grid scale module 104 is used to obtain the maximum grid scale of the local vortex direction;

[0214] A filter scale module 105 is used to obtain a shear layer adaptive grid filter scale based on the LES method according to the vortex tilt metric function and the maximum grid scale;

[0215] A shielding module 106, configured to obtain a new shielding function according to the standard shielding function, the additional shielding function and the suppressed shielding function; and

[0216] The improved calculation module 107 is used to obtain flow field data based on the improved DDES method according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

[0217] An embodiment of the present application provides a storage medium for storing computer executable instructions, characterized in that when the computer executable instructions are executed, the steps of the unsteady numerical simulation method based on improved delayed detached eddy simulation described in any one of the above embodiments are implemented.

[0218] It should be noted that the embodiments of the storage medium in this specification and the embodiments of the blockchain-based service provision method in this specification are based on the same inventive concept. Therefore, the specific implementation of this embodiment can refer to the aforementioned corresponding implementation of the blockchain-based service provision method, and the repeated parts will not be repeated.

[0219] The above is a description of a specific embodiment of the specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in an order different from that in the embodiments and still achieve the desired results. In addition, the processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0220] In the 1930s, improvements to a technology could be clearly distinguished as hardware improvements (for example, improvements to the circuit structure of diodes, transistors, switches, etc.) or software improvements (improvements to the method flow). However, with the development of technology, many improvements to the method flow today can be regarded as direct improvements to the hardware circuit structure. Designers almost always obtain the corresponding hardware circuit structure by programming the improved method flow into the hardware circuit. Therefore, it cannot be said that an improvement in a method flow cannot be implemented using a hardware entity module. For example, a programmable logic device (PLD) (such as a field programmable gate array (FPGA)) is such an integrated circuit whose logical function is determined by the user's programming of the device. Designers can "integrate" a digital system on a PLD by programming it themselves, without having to ask a chip manufacturer to design and produce a dedicated integrated circuit chip. Moreover, nowadays, instead of manually making integrated circuit chips, this kind of programming is mostly implemented by "logic compiler" software, which is similar to the software compiler used when developing and writing programs, and the original code before compilation must also be written in a specific programming language, which is called hardware description language (HDL). There is not only one HDL, but many kinds, such as ABEL (Advanced Boolean Expression Language), AHDL (Altera Hardware Description Language), Confluence, CUPL (Cornell University Programming Language), HDCal, JHDL (Java Hardware Description Language), Lava, Lola, MyHDL, PALASM, RHDL (Ruby Hardware Description Language), etc. The most commonly used ones are VHDL (Very-High-Speed ​​Integrated Circuit Hardware Description Language) and Verilog. Those skilled in the art should also know that it is only necessary to program the method flow slightly in the above-mentioned hardware description languages ​​and program it into the integrated circuit, and then it is easy to obtain the hardware circuit that implements the logic method flow.

[0221] The controller can be implemented in any appropriate manner, for example, the controller can take the form of a microprocessor or processor and a computer-readable medium storing a computer-readable program code (such as software or firmware) that can be executed by the (micro)processor, a logic gate, a switch, an application-specific integrated circuit (ASIC), a programmable logic controller, and an embedded microcontroller. Examples of controllers include, but are not limited to, the following microcontrollers: ARC 625D, Atmel AT91SAM, Microchip PIC18F26K20, and Silicone Labs C8051F320. The memory controller can also be implemented as part of the control logic of the memory. Those skilled in the art also know that in addition to implementing the controller in a purely computer-readable program code manner, the controller can be implemented in the form of a logic gate, a switch, an application-specific integrated circuit, a programmable logic controller, and an embedded microcontroller by logically programming the method steps. Therefore, this controller can be considered as a hardware component, and the devices included therein for implementing various functions can also be regarded as structures within the hardware component. Or even, the devices for implementing various functions can be regarded as both software modules for implementing the method and structures within the hardware component.

[0222] The systems, devices, modules or units described in the above embodiments may be implemented by computer chips or entities, or by products with certain functions. A typical implementation device is a computer. Specifically, the computer may be, for example, a personal computer, a laptop computer, a cellular phone, a camera phone, a smart phone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.

[0223] For the convenience of description, the above devices are described in terms of functions and are divided into various units. Of course, when implementing the embodiments of this specification, the functions of each unit can be implemented in the same or multiple software and / or hardware.

[0224] It should be understood by those skilled in the art that one or more embodiments of this specification may be provided as a method, system or computer program product. Therefore, one or more embodiments of this specification may take the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware. Moreover, this specification may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0225] This specification is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of this specification. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0226] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0227] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0228] In a typical configuration, a computing device includes one or more processors (CPU), input / output interfaces, network interfaces, and memory.

[0229] The memory may include non-permanent storage in a computer-readable medium, random access memory (RAM) and / or non-volatile memory in the form of read-only memory (ROM) or flash RAM. The memory is an example of a computer-readable medium.

[0230] Computer readable media include permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. Information can be computer readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, compact disk read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cassettes, magnetic tape magnetic disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device. As defined herein, computer readable media does not include temporary computer readable media (transitory media), such as modulated data signals and carrier waves.

[0231] It should also be noted that the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, commodity or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, commodity or device. In the absence of more restrictions, the elements defined by the sentence "comprises a ..." do not exclude the existence of other identical elements in the process, method, commodity or device including the elements.

[0232] One or more embodiments of the present specification may be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform specific tasks or implement specific abstract data types. One or more embodiments of the present specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules may be located in local and remote computer storage media, including storage devices.

[0233] Each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. In particular, for the system embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and the relevant parts can be referred to the partial description of the method embodiment.

[0234] The above description is only an embodiment of this document and is not intended to limit this document. For those skilled in the art, this document may have various changes and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of this document should be included in the scope of the claims of this document.

Claims

1. An unsteady numerical simulation method based on improved delayed detached eddy simulation, characterized in that include: Constructing vortex tilt measurement function; Obtain turbulent kinetic energy and turbulent specific dissipation rate in real time; Based on the RANS method, a turbulence length scale is obtained according to the turbulent kinetic energy and the turbulence specific dissipation rate; The maximum grid size for obtaining the local vortex direction; Based on the LES method, a shear layer adaptive grid filter scale is obtained according to the vortex tilt metric function and the maximum grid scale; According to the standard shielding function, the additional shielding function and the suppressed shielding function, a new shielding function is obtained; as well as Based on the improved DDES method, flow field data are obtained according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

2. The method according to claim 1, characterized in that Also includes: Using the dissipation detector, determine whether the flow field is a smooth area or a discontinuous area; If the flow field is a smooth region, the symmetric reconstruction variables are used; If the flow field is a discontinuous region, use monotonic reconstruction variables; and A mixed reconstruction variable is obtained according to the symmetric reconstruction variable and the monotonic reconstruction variable.

3. The method according to claim 2, characterized in that The hybrid reconstruction variable calculation formula is as follows: in, represents a mixed reconstruction variable; q L,R represents the traditional monotone reconstruction variable; ψ represents the dissipative control function; and represents a symmetric reconstruction variable.

4. The method according to claim 1, characterized in that: The shear layer adaptive grid filter scale obtained based on the vortex tilt metric function and the maximum grid scale includes: Obtain velocity trace, vortex vector, vorticity magnitude, kinematic viscosity coefficient and turbulent kinematic viscosity coefficient; Constructing a vortex inclination measurement function according to the velocity trace, the vortex vector, the vorticity magnitude, the kinematic viscosity coefficient and the turbulent kinematic viscosity coefficient; According to the vortex inclination measurement function, an instability coefficient function is obtained; and The shear layer adaptive grid filter scale is determined according to the instability coefficient function and the maximum grid scale.

5. The method according to claim 4, characterized in that The calculation formula of the shear layer adaptive grid filter scale is as follows: Among them, Δ SLA represents the scale of the shear layer adaptive grid filter; represents the maximum grid size based on the local vortex direction; F KH represents the Kelvin-Helmholtz instability coefficient function; i and j represent the numbers of units i and j; nb(i) represents the neighbor unit number of unit i; and I ij Represents the length vector of the line connecting the centers of cells i and j.

6. The method according to claim 1, characterized in that The new shielding functions obtained based on the standard shielding function, additional shielding function and suppression shielding function include: Using the data measured by the outer boundary layer detector, an additional shielding function is constructed; Using the data measured by the shear layer detector, a suppression function is constructed; Using data measured by the basic detector, a standard shielding function is constructed; and A new shielding function is obtained according to the additional shielding function, the suppression function and the standard shielding function.

7. The method according to claim 6, characterized in that The calculation formula of the novel shielding function is as follows: f P =f d ·(1-(1-f P2 )·f R ) Among them, f P Represents a new type of shielding function; f d Represents the standard masking function; f P2 represents an additional masking function; and f R Represents the suppression function.

8. The method according to claim 1, characterized in that The calculation formula of the improved DDES method is as follows: L DDES l RANS -f p max{0,l RANS -C DES Δ SLA } Among them, l DDES represents the turbulence length scale of the DDES method; l RANS represents the turbulence length scale of the RANS method; f p Represents a new type of shielding function; C DES is the calibration factor for the DDES method; and Δ SLA Represents the scale of the shear layer adaptive grid filter.

9. An unsteady numerical simulation device based on improved delayed detached eddy simulation, characterized in that include: Metric module, used to construct vortex tilt metric function; An acquisition module is used to obtain turbulent kinetic energy and turbulent specific dissipation rate in real time; A RANS scale module, for obtaining a turbulence length scale based on a RANS method according to the turbulent kinetic energy and the turbulence specific dissipation rate; The grid scale module is used to obtain the maximum grid scale of the local vortex direction; A filter scale module, used for obtaining a shear layer adaptive grid filter scale based on the LES method according to the vortex tilt metric function and the maximum grid scale; A shielding module, used for obtaining a new shielding function according to a standard shielding function, an additional shielding function and a suppressed shielding function; as well as The improved calculation module is used to obtain flow field data based on the improved DDES method according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale.

10. A storage medium for storing computer executable instructions, characterized in that: When the computer executable instructions are executed, the steps of the unsteady numerical simulation method based on improved delayed detached eddy simulation according to any one of claims 1 to 8 are implemented.

Images