An unsteady numerical simulation method based on improved delayed detached eddy simulation

Abstract

The embodiments of this specification provide an unsteady numerical simulation method, device, and storage medium based on an improved delayed detached eddy simulation, wherein the method includes: constructing an eddy tilt measurement function; obtaining turbulent kinetic energy and turbulent specific dissipation rate in real time; obtaining a turbulent length scale based on the turbulent kinetic energy and turbulent specific dissipation rate based on the RANS method; obtaining the maximum grid scale of the local vortex direction; obtaining a shear layer adaptive grid filter scale based on the eddy tilt measurement function and the maximum grid scale based on the LES method; obtaining a new shielding function based on the standard shielding function, the additional shielding function, and the suppressed shielding function; and obtaining flow field data based on the shear layer adaptive grid filter scale, the new shielding function, and the turbulent length scale based on the improved DDES method. The improved DDES method can improve the resolution of unsteady large-separation turbulent flow structures in a wide speed domain and improve the computational accuracy of unsteady flow simulation with large separation in a wide speed domain.

Description

Technical Field

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

[0002] In the aerospace field, aircraft encounter a variety of complex unsteady flow phenomena when flying across a wide range of speeds. Numerical simulations can predict and evaluate the impact of these phenomena on aircraft performance, thereby guiding optimal design and improving aircraft stability and safety.

[0003] The existing technology mainly uses the DES / DDES hybrid method 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 the 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 in turn induce 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 a shielding function f d Redefining the length scale avoids the MSD problem to a great extent. The DDES method has achieved considerable success. However, in-depth research shows that when the flow grid size is less than 0.3 times the local boundary layer thickness δ, the screening 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 the grid 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 shielding function can be optimized by adjusting the coefficients to improve boundary layer protection. However, studies have also found that inappropriate coefficients can lead to overprotection, delaying the development of mixing layer instabilities and intensifying the second type of gray zone-induced problems. Furthermore, different coefficient values ​​may need to be calibrated for different flow phenomena and mesh densities, significantly increasing the difficulty of applications targeting complex geometries.

[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 faster in the freestream shear layer simulation, thus minimizing the delay in the construction of freestream shear layer instabilities. However, this also makes it more difficult for the aforementioned shielding function to protect the attached flow boundary layer, making the first type of gray zone induced problem 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] Construct 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] Obtaining a new shielding function according to the standard shielding function, the additional shielding function and the suppressed shielding function; 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, symmetric reconstruction variables are used;

[0018] If the flow field is discontinuous, 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 monotonic reconstruction variable;

[0024] ψ represents the dissipative control function; and

[0025] represents the symmetric reconstruction variable.

[0026] Furthermore, 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] Obtaining an instability coefficient function according to the vortex tilt measurement function; 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 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 shielding 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 an improved delayed detached eddy simulation, comprising:

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

[0057] Acquisition module, used to obtain turbulent kinetic energy and turbulent specific dissipation rate in real time;

[0058] A RANS scale module is used to obtain 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 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;

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

[0062] An improved calculation module is used to obtain flow field data based on an 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 existing technology, 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 the simulation of unsteady flows with 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 following briefly introduces the drawings required for use in the embodiments or the description of the prior art. 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 any creative work.

[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 Schematic diagram of subsonic back-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 in one or more embodiments of this specification;

[0071] Figure 5 A schematic diagram comparing time-averaged flow velocity profiles at different x-sections provided in one or more embodiments of this specification;

[0072] Figure 6 A schematic diagram comparing the average shear stress u'u' at different x-sections provided in one or more embodiments of this specification;

[0073] Figure 7 A schematic diagram comparing the average shear stress u'v' at different x-sections provided in 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] Figure 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] Figure 10 A schematic diagram of the velocity profile for the first 1 mm of the cylinder bottom provided for one or more embodiments of this specification;

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

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

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

[0080] Figure 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] Figure 15 Schematic diagram of time-averaged radial velocity of different flow sections under dense grid provided for one or more embodiments of this specification

[0082] Figure 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] Figure 17 A schematic diagram of turbulent shear stress in different flow direction sections under a dense grid provided in one or more embodiments of this specification;

[0084] Figure 18 A schematic structural diagram 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 technical solutions in one or more embodiments of this specification will be clearly and completely described below in conjunction with the drawings 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 making creative efforts should fall within the scope of protection of this document.

[0086] Currently, existing turbulence simulation technologies primarily use a RANS-LES hybrid method for calculations. The RANS-LES hybrid method combines the advantages of both RANS and LES. The RANS method is used for simulation near the wall to reduce computational costs, while the LES method is used for simulation far from the wall to capture the detailed characteristics of turbulence. However, during the conversion from RANS (Reynolds-Averaged Navier-Stokes) to LES (Large Eddy Simulation), existing technologies produce gray areas in the attached flow boundary layer, resulting in a flow field solution that is neither a pure RANS solution nor a LES solution, leading to low flow field calculation accuracy.

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

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

[0089] In an embodiment of the present application, the velocity trace, vortex vector, vortex size, kinematic viscosity coefficient and turbulent kinematic viscosity coefficient are obtained; and a vortex inclination measurement function is constructed based on the velocity trace, the vortex vector, the vortex 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 indicating velocity;

[0096] ν represents the kinematic viscosity coefficient;

[0097] ν t represents the turbulent viscosity coefficient.

[0098] Step S2: acquiring 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 a turbulence model (such as a k-omega model, LES, etc.).

[0100] For example, in RANS methods, 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 solves two differential equations to close the time-averaged governing equations for turbulence, thereby obtaining the distribution of the turbulent kinetic energy k.

[0101] Step S3: Based on the RANS method, the turbulent length scale is obtained according to the turbulent kinetic energy and the turbulent 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 of 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, and the flow characteristics cannot be accurately described.

[0110] Step S5: 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.

[0111] In an embodiment of the present application, a vortex tilt metric function is calculated within the entire computational domain in a large eddy simulation (LES) to obtain a vortex tilt metric value. This vortex tilt metric value is then used to capture key turbulence characteristics in the shear layer. Based on the size and distribution of the vortex tilt metric value, an adaptive mesh refinement criterion is developed that better resolves the turbulence structure in the shear layer. Based on the adaptive mesh refinement criterion and the maximum mesh size, an adaptive mesh filter scale is then calculated for each mesh cell.

[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 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 metric function;

[0122]

[0123] ω represents the vortex vector;

[0124] Ω represents the vorticity size;

[0125] A trace indicating velocity;

[0126] ν represents the kinematic viscosity coefficient;

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

[0128] The subscripts i and j indicate the numbers of 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 a boundary layer outer layer detector; an inhibition 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 inhibition 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; an inhibition 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; and a standard shielding function is constructed by detecting the inner layer data of the boundary layer through the detector. d Based on the coupling, additional shielding function f p2 and the suppression function f R , obtain the 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 shielding 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 Indicates 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 composed of α;

[0161] α represents the Intermediate variable combined with c4.

[0162] In this 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 an improved DDES (Delayed DES) method, 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 of 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, a symmetrical reconstruction variable is used; if the flow field is a discontinuous area, a monotonic reconstruction variable is used; and based on the symmetrical reconstruction variable and the monotonic reconstruction variable, a mixed reconstruction variable is obtained.

[0173] The smooth area of ​​the flow field requires low dissipation and requires symmetric reconstruction variables To maintain computational stability in discontinuous areas 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 the 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 monotonic 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 the symmetric reconstruction variable;

[0186]

[0187] β represents the coefficient of the symmetric reconstruction format, which is usually 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 monotonic 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 between smooth and discontinuous regions within the LES. In smooth regions, the flow is relatively stable, allowing for lower resolution and less computational resources. In contrast, in discontinuous regions, the flow varies dramatically, requiring higher resolution and more computational resources. The dissipative detector dynamically adjusts computational resources based on regional characteristics, improving computational efficiency and resolving unsteady, widely separated 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. Simulation of subsonic back-step separation flow

[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, as shown in Figure 2 As shown; where H is the step height; M ∞ is the incoming flow Mach number; W is the step width, which here 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 refined. The total number of grid elements is about 7 million. Figure 4 The x-section shown extracts the corresponding flow field information, including the average 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. As can be seen from the figure, 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 in the rear of the step, the improved DDES method has the best results, and is in good agreement with the experiment at all x-sections. Figure 5 It can also be seen that the results of the standard DDES method are slightly worse at x / H=1 and x / H=6; 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 time-averaged shear stress distributions u'u' and u'v' at different x-section locations are given. The figure again shows that the improved DDES method significantly outperforms the standard DDES method. The improved DDES method provides better agreement with the experimental results for both the peak and shape of shear stress at each section. The standard DDES method significantly underestimates the shear stress at x / H = 1 and overestimates the shear stress downstream at x / H = 6 due to the hysteresis of shear layer instability.

[0203] Finally, through the unsteady simulation of subsonic back-step separation flow, it is shown that the improved DDES method using a new shielding function and shear layer adaptive grid filtering can provide statistically more consistent results with the experiment in terms 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 pulsation quantity.

[0204] 2. Simulation of supersonic bottom separation flow in 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 calculations used two grids, a sparse grid and a dense grid, to analyze the grid's influence. The dense grid had approximately 9.8 million cells, while the sparse grid had approximately 4.7 million cells. The flow distributions of the two grids were essentially identical, with the sparse grid having approximately half the number of circumferential points as the dense grid. Figure 8 The distribution of the symmetry surface mesh is given; Figure 9 The spatial distribution of two grids at the bottom section of the cylinder is presented: a dense grid and a sparse grid. The RANS method uses the SST turbulence model, while the RANS-LES hybrid method uses the SST-DDES method based on the SST turbulence model. For the unsteady DDES simulation, a low-dissipation numerical scheme developed in this project was used, with a physical time step of 1.0e-6s and a statistical time of 0.04s for the time-averaged flow field analysis. Figure 10 The comparison between the experimental results and the calculation of the velocity profile in the first 1mm of the bottom is given. As can be seen from the figure, the inlet velocity profiles obtained with different grids and different methods are in good agreement with the experiment, ensuring that the initial conditions of the separated flow at the bottom are consistent with the experiment. Figure 11 (standard DDES) and Figure 12 (Improved DDES) It can be seen that the improved DDES method can obtain a richer flow field structure than the standard DDES method, mainly due to the improvement of the grid filtering method. The improved DDES method obviously promotes shear layer instability more quickly. The time-averaged field and the pulsating field are further analyzed based on the calculation results. Figure 13 Four flow section locations are given for analysis. Figure 14 and Figure 15 The time-averaged streamwise and radial velocity distributions at different flow section locations on the dense grid are given. As can be seen from the figures, the improved DDES method simulates the velocity distribution within the shear layer quite well, significantly outperforming the standard DDES method. On the other hand, the SST turbulence model also simulates the average velocity field well.

[0206] from Figure 14 (a) It can be seen that the streamwise velocity predicted by the improved DDES method is slightly too large 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 size may not meet the simulation requirements for the initial development of the shear layer. Figure 14 It can be seen again that the improved DDES method predicts a larger recirculation velocity at the centerline. 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, Figure 16 and Figure 17 The distributions of turbulent kinetic energy k and turbulent shear stress at different flow-direction cross-sectional locations under a dense grid are given. As can be seen from the figure, the standard DDES method does not accurately predict the pulsation magnitude at any cross-sectional location; the magnitudes are too small, and the location of the maximum pulsation deviates significantly from the experimental results. The improved DDES method significantly improves the calculation of unsteady pulsations. For the section x / R = 0.1575, the previous analysis indicated that the mesh required for simulating the flow characteristics may not be sufficient, and the pulsation calculated by the improved DDES method is slightly smaller. For the section x / R = 0.9449, the pulsation in the central region is also greater due to the larger predicted backflow velocity.

[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 more accurately obtains the position of the stagnation point after the flow; (2) Whether it is time-averaged or pulsating quantity, 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 than the standard DDES method for simulating flows with complex shapes.

[0209] The embodiment of the present application provides an unsteady numerical simulation device based on improved delayed detached eddy simulation, such as Figure 18 Shown, including:

[0210] Metric module 101, used to construct a vortex tilt metric 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 is configured to obtain 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 configured 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 for obtaining 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 embodiment of the storage medium in this specification and the embodiment 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 implementation of the corresponding blockchain-based service provision method mentioned above, and the repeated parts will not be repeated.

[0219] The foregoing description of this specification describes specific embodiments. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims can be performed in an order different from that described in the embodiments and still achieve the desired results. Furthermore, the processes depicted in the accompanying drawings do not necessarily require the specific order shown or the sequential order to achieve the desired results. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0220] In the 1930s, technological improvements could be clearly distinguished as either hardware improvements (for example, improvements to circuit structures like diodes, transistors, and switches) or software improvements (improvements to process flows). However, with the advancement of technology, many process flow improvements today can now be considered direct improvements to hardware circuit structures. Designers almost always create the corresponding hardware circuit structure by programming the improved process flow into the hardware circuit. Therefore, it cannot be said that a process flow improvement cannot be implemented using hardware modules. For example, a programmable logic device (PLD), such as a field programmable gate array (FPGA), is an integrated circuit whose logical function is determined by user programming. Designers can "integrate" a digital system on a PLD by programming it themselves, without having to hire a chip manufacturer to design and manufacture a dedicated integrated circuit chip. Moreover, nowadays, instead of manually fabricating integrated circuit chips, this programming is mostly done using "logic compiler" software. This is similar to the software compiler used when developing programs. Before compilation, the original code must also be written in a specific programming language, called a hardware description language (HDL). There is not just one HDL, but many, 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 will also understand that by simply programming the method flow in one of these hardware description languages ​​and then programming it into an integrated circuit, a hardware circuit that implements the logic method flow can be easily obtained.

[0221] The controller can be implemented in any suitable manner. For example, the controller can take the form of a microprocessor or processor and a computer-readable medium storing computer-readable program code (e.g., software or firmware) executable by the (micro)processor, logic gates, switches, application-specific integrated circuits (ASICs), programmable logic controllers, and embedded microcontrollers. 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 will also know that in addition to implementing the controller in a purely computer-readable program code format, the controller can be implemented in the form of logic gates, switches, application-specific integrated circuits, programmable logic controllers, and embedded microcontrollers by logically programming the method steps. Therefore, such a controller can be considered a hardware component, and the devices included therein for implementing various functions can also be considered as structures within the hardware component. Or even, the devices for implementing various functions can be considered as both software modules that implement 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 having 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 smartphone, 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 as being divided into various units according to their functions. 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] Those skilled in the art will appreciate that one or more embodiments of this specification may be provided as a method, system, or computer program product. Thus, one or more embodiments of this specification may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, 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 magnetic 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 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 that can direct a computer or other programmable data processing device to work 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 The 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 operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. 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 (CPUs), input / output interfaces, network interfaces, and memory.

[0229] 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. Memory is an example of a computer-readable medium.

[0230] Computer-readable media includes permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. The 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 disc read-only memory (CD-ROM), digital versatile disc (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 transitory computer-readable media (transitory media), such as modulated data signals and carrier waves.

[0231] It should also be noted that the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, commodity, or apparatus that includes a series of elements includes not only those elements but also other elements not explicitly listed, or includes elements inherent to such process, method, commodity, or apparatus. In the absence of further limitations, an element defined by the phrase "comprises a ..." does not exclude the presence of other identical elements in the process, method, commodity, or apparatus that includes the element.

[0232] One or more embodiments of this 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, and the like that perform specific tasks or implement specific abstract data types. One or more embodiments of this specification may also be practiced in distributed computing environments where tasks are performed by remote processing devices connected via a communications network. In a distributed computing environment, program modules may be located in local and remote computer storage media, including storage devices.

[0233] The various embodiments in this specification are described in a progressive manner. Similar parts between the various embodiments can be referred to in conjunction with each other. Each embodiment focuses on the differences between the other embodiments. In particular, the system embodiments are generally similar to the method embodiments, so the description is relatively simple. For relevant parts, refer to the description of the method embodiments.

[0234] The foregoing description is merely an example of the present invention and is not intended to limit the present invention. Persons skilled in the art will readily appreciate that various modifications and variations of the present invention are possible. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be included within the scope of the claims herein.

Claims

1. An unsteady numerical simulation method based on improved delayed detached eddy simulation, characterized in that include: Construct 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 is obtained according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale; The additional shielding function is constructed by the turbulent viscosity gradient detector, and the suppression function is constructed by the vortex normal gradient detector; the inner layer data of the boundary layer is detected by the detector to construct the standard shielding function; in the standard shielding function f d Based on the coupling, additional shielding function f p2 and the suppression function f R , obtain the new shielding function f p ; 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.

2. The method according to claim 1, characterized in that Also includes: Using the dissipation detector, the flow field is determined to be a smooth area or a discontinuous area; If the flow field is a smooth region, symmetric reconstruction variables are used; If the flow field is discontinuous, 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 calculation formula of the hybrid reconstruction variable is as follows: in, represents a mixed reconstruction variable; q L,R represents the traditional monotonic reconstruction variable; ψ represents the dissipative control function; and represents the symmetric reconstruction variable.

4. The method according to claim 1, wherein The shear layer adaptive grid filter scale is obtained based on the vortex tilt metric function and the maximum grid scale, including: 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; Obtaining an instability coefficient function according to the vortex tilt measurement function; 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 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 new 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 shielding function; f P2 represents an additional masking function; and f R represents the suppression function.

8. An unsteady numerical simulation device based on improved delayed detached eddy simulation, characterized in that include: Metrics module, used to construct vortex tilt metric function; Acquisition module, used to obtain turbulent kinetic energy and turbulent specific dissipation rate in real time; A RANS scale module is used to obtain 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 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; 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 An improved calculation module is used to obtain flow field data based on an improved DDES method according to the shear layer adaptive grid filter scale, the new shielding function and the turbulence length scale; The additional shielding function is constructed by the turbulent viscosity gradient detector, and the suppression function is constructed by the vorticity normal gradient detector; the inner layer data of the boundary layer is detected by the detector to construct the standard shielding function; in the standard shielding function f d Based on the coupling, additional shielding function f p2 and the suppression function f R , obtain the new shielding function f p ; 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. 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 are implemented according to any one of claims 1 to 7.

Images