A method for predicting rotating turbulent flows using the BSTMQ turbulence model.

Abstract

A method for predicting rotating turbulent flows using a BSTMQ turbulence model is provided, which can improve the effectiveness of predicting the internal flow field of fluid machinery. The method includes the steps of performing solid modeling of a rotating turbulent flow prediction target by 3D modeling software to obtain a simulation model, meshing the simulation model, improving the eddy viscosity damping function of an existing STRUCT turbulence model to construct a BSTMQ turbulence model, and performing a numerical simulation calculation of the rotating turbulence using the BSTMQ turbulence model by a calculation program based on OpenFoam to obtain a prediction result of the rotating turbulence of the prediction target. The BSTMQ model can better simulate the rotating turbulent flow, improve the prediction effect of the internal flow field of fluid machinery, and provide a new idea for efficiently and accurately calculating the rotating turbulent flow in fluid machinery.

Description

[Technical field]

[0001] The present invention relates to turbulence simulation in engineering fluid mechanics, and more particularly to a method for predicting rotating turbulence using the BSTMQ turbulence model. [Background technology]

[0002] In recent years, turbulent mixed flow models that take into account both computational accuracy and computational cost have shown promising future applications. In addition to various turbulent mixed flow models based on grid spatial scales, the STRUCT model proposed by Professor Emilion based on time scales can achieve monotonic convergence of the grid and has excellent computational accuracy and robustness.

[0003] However, in fluid machinery, the impeller rotates at high speed, the blades twist, and the wall surface is complex, so the internal flow field is subject to strong rotational and curvature effects. The existing STRUCT model cannot be directly applied and has the following problems: (1) The eddy viscosity damping of the STRUCT model is close to 1 in the near-wall region, so the model reverts to the traditional Reynolds-averaged Navier-Stokes (RANS) model, which is not suitable for analyzing the near-wall region of rotating machinery. (2) In the existing STRUCT model, the analysis time scale calculated based on the traditional vortex detection criterion does not reflect the rotational and curvature effects. The above deficiencies may cause the model to be unable to accurately predict the development process of rotating turbulence. Summary of the Invention [Problem to be solved by the invention]

[0004] Object of the Invention: In response to the above problems, the present invention provides a method for predicting rotating turbulent flows using a BSTMQ turbulence model, which can improve the prediction effect of the internal flow field of fluid machinery. [Means for solving the problem]

[0005] Technical solution: In order to solve the above problems, the present invention provides a rotating turbulent flow prediction method using the BSTMQ turbulence model, comprising: (1) performing solid modeling of a target for prediction of rotating turbulent flow using 3D modeling software to obtain a simulation model; (2) meshing the simulation model; A step of constructing a BSTMQ turbulence model by performing rotational curvature correction based on a bifurcation method on an existing STRUCT turbulence model. 、B STMQ turbulence model equation is as follows:

number

number

number

number

number

[0006] Furthermore, in the step (3), the step of constructing the BSTMQ turbulence model specifically includes the steps of: Step (31): Introduce the strain rate tensor and the rotation rate tensor in the rotating coordinate system; Step (32): Spalart -Shur tensor is used to unify the rotation and curvature effects to obtain a corrected rotation rate tensor, Step (33): Create a velocity gradient invariant from the strain rate tensor introduced in step (31) and the corrected rotation rate tensor obtained in step (32); Step (34): Determine a corrected analysis time scale from the velocity gradient invariant obtained in step (33); Step (35): Determine the modeling time scale for rotating turbulence from the specific dissipation rate; Step (36): According to steps (34) and (35), calculate the eddy viscosity damping function D, which reflects the rotation and curvature effects. f-b Build

number

number

number

number

number

[0007] Furthermore, in the inertial coordinate system, the strain rate tensor S ij and the rotational velocity tensor Ω ij The calculation formula is as follows:

number

number

number

[0008] moreover, Spalart By unifying the rotation effect and the curvature effect based on the -Shur tensor, the corrected rotation rate tensor Ω is obtained as follows: mod ij The formula for calculation is obtained,

number

[0009] In addition, the correction variable W A ij The decision was For two-dimensional flow, W A ij =Ω SS ij And, TIFF0007621020000017.tif11170

number

number

number

[0010] moreover,

number

number

[0011] Furthermore, the corrected analysis time scale is:

number

number

number

number

[0012] Furthermore, the conventional eddy viscosity damping function D f The eddy viscosity damping function D reflects the rotation and curvature effects based on f-b and construct the conventional eddy viscosity damping function D f teeth,

number

[0013] The present invention also provides a computer apparatus including a memory, a processor, and a computer program stored in the memory and operable on the processor, the computer apparatus comprising: The processor employs a computing device that, when executing the computer program, implements the steps of the above-described method.

[0014] The present invention also provides a computer-readable storage medium having a computer program stored thereon, The computer program employs a computer readable storage medium, characterized in that, when executed by a processor, it implements the steps of the above-described method. Effect of the Invention

[0015] The beneficial effects are as follows: Compared with the prior art, the key advantages of the present invention are that the BSTMQ model can better simulate rotating turbulent flows, improve the prediction effect of the internal flow field of fluid machinery, and provide a new idea for efficiently and accurately calculating rotating turbulent flows. The STRUCT model with rotational curvature correction based on the bifurcation method uses the BSkO model as the base model and improves the analysis ability of the near-wall region. Spalart By unifying the rotation and curvature effects using the -Shur tensor and constructing a new analytical time scale, the constructed BSTMQ model can better simulate rotating turbulence. [Brief description of the drawings]

[0016] [Figure 1] 1 is a schematic flowchart of a rotating turbulent flow prediction method in the present invention. [Diagram 2] 1 is a schematic flow chart of a turbulence model compilation process according to the present invention. [Diagram 3] FIG. 2 is a schematic diagram of the computational domain for an instance of rotating channel flow in accordance with the present invention; [Figure 4] FIG. 2 is a schematic diagram of the computational domain of an instance of a rotating-flow sudden expansion in accordance with the present invention. [Diagram 5] FIG. 11 is a diagram comparing average velocity distributions of instances of rotating channel flow in the present invention. [Figure 6] FIG. 10 is a graph comparing dimensionless normal pulsation velocities of instances of rotating channel flow in the present invention. [Figure 7] FIG. 11 is a diagram comparing the resolution of instances of a rotating channel flow in the normal direction according to the present invention. [Figure 8] FIG. 11 is a comparison of the average axial velocity pulsation root mean square at Z / D=0.25 for instances of a rotating flow sudden expansion pipe according to the present invention. [Figure 9] FIG. 11 is a graph comparing the root mean square of the circumferential average velocity pulsation at the Z / D=4 position of instances of a rotational flow sudden expansion pipe in the present invention. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0017] As shown in FIG. 1, in this embodiment, the rotating turbulent flow simulation method using the BSTMQ model includes the following steps (1) to (4).

[0018] Step (1): Using 3D modeling software, solid modeling is performed on the object of rotational turbulence prediction (classical basic instances such as rotating channel flow and rotating flow sudden expansion pipes, and rotating machines such as centrifugal pumps and axial pumps) to obtain a simulation model.

[0019] Step (2): OpenFOAM (Registered Trademark) Mesh the simulation model using proprietary or third-party software.

[0020] Step (3): A rotational curvature correction is performed based on the bifurcation method for the existing STRUCT turbulence model to construct a BSTMQ turbulence model.

[0021] Step (4): OpenFOAM (Registered Trademark) Based on the program of the rotating turbulent flow numerical simulation method using the flatbed BSTMQ turbulence model, the numerical simulation calculation was carried out, and finally, the OpenFOAM (Registered Trademark) Post-processing using proprietary or third-party software provides a more accurate prediction of the rotating turbulent flow being predicted.

[0022] The BSTMQ turbulence model is a STRUCT model corrected based on the rotational curvature of the bifurcation method, and the steps are specifically as follows:

[0023] Step (31): Introduce the strain rate tensor and the rotation rate tensor in the rotating coordinate system.

[0024] In this step, we define the strain rate tensor S ij and the rotational velocity tensor Ω ij is obtained from the following formula:

number

[0025] Here, i and j are two free indices in the tensor representation, which take on the values ​​1, 2, and 3, respectively, and represent the need to traverse all values ​​of the index. i and U j is the velocity vector, and the velocity components in the three directions are U 1 =U,U 2 =V, U 3 = W, and x i and x j are the coordinate axes, and x 1 =x, x 2 =y,x 3 = y,

number

[0026] Angular velocity vector Ω F k For a coordinate system rotating with , the rotation rate tensor is defined as

number

[0027] TIFF0007621020000032.tif7170

[0028] Step (32): Spalart -Unifying rotation and curvature effects using the Shur tensor.

[0029] In this step, we compute the Spalart-Shur tensor (Ω SS ij ) is obtained by the following formula:

number

[0030] Where:

number

[0031] Unifying the rotation and curvature effects gives the corrected rotation rate tensor.

number

[0032] Here, C r is a model constant, and W A ij teeth Spalart - is a correction variable based on the Shur tensor, For two-dimensional flow,

number

number

number

number

number

[0033] Tensor X ij are intermediate variables, where δ ij is the Kronecker symbol tensor, and S ik and S kj physically it also represents the strain rate tensor, and the indices are changed to represent tensor operations. S ik S kj has two free indexes and one summation index, so it must first sum and then traverse all components.

[0034] II S and III S are obtained by the following formulas, respectively.

number

[0035] Here, S 11 , S 22 and S 33 are the three diagonal components of the strain rate tensor.

[0036] Step (33): Create a dimensionless velocity gradient invariant from the strain rate tensor introduced in step (31) and the corrected rotation rate tensor in step (32), thereby ensuring that the subsequent model equations are Galilean invariant and can reflect the effects of rotation and curvature.

[0037] In this step, the velocity gradient invariant is given by:

number

[0038] where k is the turbulent kinetic energy, ε is the dissipation rate, S is the mode of the strain rate tensor, and Ω mod are the modes of the corrected rotation rate tensor,

number

number

number

[0039] Step (34): From the velocity gradient invariant of step (33), we obtain the analytical time scale t r Determine.

[0040] The Q vortex detection criterion considering the rotation effect and curvature effect is given by the following equation:

number

[0041] where k is the turbulent kinetic energy, ε is the dissipation rate, S and Ωmod are the modes of the strain rate tensor and the corrected rotation rate tensor, respectively.

[0042] The corrected analysis time scale is:

number

[0043] Step (35): Modeling time scale t of rotating turbulence from specific dissipation rate m Determine the modeling time scale t m is obtained by the following formula:

number

[0044] where α is an empirical constant for the modeled time scale, ω is the specific dissipation rate,

number

[0045] Step (36): According to steps (34) and (35), calculate the eddy viscosity damping function D, which reflects the rotation and curvature effects. f-b Build.

[0046] Conventional eddy viscosity damping function D f is as follows.

number

[0047] Finally, the eddy viscosity damping function D, which reflects the rotation and curvature effects, f-b Build.

number

[0048] Step (37): Construct the BSTMQ turbulence model according to the eddy viscosity decay function in step (36).

[0049] The turbulent viscosity is corrected using the BSkO model as the base model. The equation for the eddy viscosity coefficient, constructed using the invariant of the velocity gradient tensor, is given as follows:

number

[0050] There are two basic constraints that the expression for eddy viscosity must satisfy: (1) the linear constitutive relationship must be maintained, and (2) for homogeneous shear flow in an inertial reference frame, the model simplifies to its original form, which is a constant.

[0051] Adjusting the model constants in the equation through the bifurcation diagram gives:

number

[0052] Finally, the equations of the BSTMQ turbulence model based on the BSkO model are as follows:

number

[0053] OpenFOAM (Registered Trademark) The calculation programs based on OpenFOAM are as follows: (Registered Trademark) OpenFOAM is an open source CFD computing platform that uses the finite volume method to discretize the equations and an unstructured grid to discretize the computational domain. (Registered Trademark) OpenFOAM already comes equipped with a relatively complete set of solvers, turbulence models, boundary conditions and post-processing procedures, and can calculate the most complex flows. (Registered Trademark)is divided into steady-state and unsteady calculations, where the steady-state calculation does not include time discretization, and the unsteady calculation time discretization is divided into a series of time steps according to the user's instructions. The specific time steps are usually estimated based on the Courant-Friedrichs-Lewy (CFL) number. The spatial discretization of the computational domain is a meshing, and OpenFOAM (Registered Trademark) You can use your own meshing tools or import computational domain grids from third-party software. The resulting discrete equations can be meshed using OpenFOAM. (Registered Trademark) The pressure can be calculated using three pressure correction methods provided by: SIMPLE (Semi-Implicit Method for Pressure Linked Equations), PISO (Pressure Implicit with Splitting of Operator), and PIMPLE (Pressure Implicit Method for Pressure Linked Equations). The SIMPLE algorithm is used for steady calculations, and the PISO or PIMPLE algorithm is usually selected for unsteady calculations of incompressible turbulent flows. When solving algebraic equations, the pressure equation is generally an antisymmetric matrix, while the velocity, turbulent kinetic energy, etc. are symmetric matrices. Since the asymmetric matrix is ​​not a diagonally dominant matrix and is prone to divergence during calculation, a multigrid preconditioner is used during calculation to convert it to a diagonally dominant matrix, then a diagonally incomplete Gauss-Seidel smooth solver is used, and finally, a conditional conjugate gradient solver is selected to calculate the pressure equation. When calculating symmetric matrices such as velocity, a smooth solver is used to calculate the velocity equation.

[0054] As shown in Figure 2, the turbulence model is compiled using OpenFOAM. (Registered Trademark)All turbulence models in share the base class turbulenceModel, and three different subclasses are defined: Laminar, RASModel, and LESModel. These three subclasses represent laminar, RANS, and LES models, respectively. Since the BSTMQ model only corrects for the turbulent viscosity of the RANS model, new model programs are compiled under the RASModel subclass.

[0055] As shown in Figure 2, the dotted box in the middle represents the implementation process of steps (31) and (32), the right side represents the implementation process of steps (33) to (36), and the right side represents the implementation process of the foundation BSkO model in step (37).

[0056] Spalart When using the -Shur tensor to unify the rotation and curvature effects and construct an invariant for the velocity gradient tensor, specifically, we first define a symmetric tensor field whose data is stored at the center of the body, that is, the strain rate tensor. Then, II S =S ij S ji and III S =S ij S jk S ki The calculation is based on the formula that constructs the intermediate variable X ij is a very complicated symmetric tensor. First, X ij We need to create an empty symmetric tensor field according to the dimension of W, then calculate the result for each component of this tensor separately and replace the empty symmetric tensor field. ij is also defined in this way, and the corrected rotational velocity tensor Ω mod and recalculate it to obtain the corrected rotational velocity tensor invariant

number

[0057] When constructing a BSkO-based model that takes into account the rotation and curvature effects based on the invariant of the velocity gradient tensor constructed in step (33), the specific procedure is as follows: First, a private member function: tmp <volscalarfield>Declare f()const, this function is declared as a template class of scalar fields, and f() is to return the result of the calculation of the eddy viscosity equation. This private function variable can only be called by member functions.

[0058] To improve the eddy viscosity damping function of the existing STRUCT model, the following specific BSTMQ model is proposed. As mentioned above, the private member function: tmp <volscalarfield>Declare Df()const, this function is declared as a template class of scalar field, Df() is to return the calculation result of eddy viscosity damping function formula. Next, multiply the returned calculation result of eddy viscosity coefficient formula f() and eddy viscosity damping function formula Df() by the original eddy viscosity formula. Recalculate the turbulent kinetic energy equation and specific dissipation rate equation using this eddy viscosity, and re-update the private functions f() and Df() based on the results obtained at each step.

[0059] After building the model, OpenFOAM (Registered Trademark) Then, the eddy and dynamic viscosities are combined into an effective viscosity, the turbulent stresses are calculated, and the solver calls the results of the turbulence model using the turbulence→divDevReff(U) statement.

[0060] OpenFOAM (Registered Trademark) Based on the unique functions and the custom compiled turbulence model, a program for numerical simulation of rotating turbulent flows using the BSTMQ model can be implemented.

[0061] To better demonstrate the advantages of the embodiment of the present invention, the SSTk-ω model, which is commonly used in current rotational turbulent flow calculations in fluid engineering, is compared with the conventional STRUCT model (STSST). As shown in Figures 3 and 4, a rotating channel flow and a rotating flow sudden expansion pipe are taken as examples, respectively, and the finite volume method is used to discretize the control equations, and a high-quality hexahedral grid is used to spatially discretize the calculation domain. The total number of grids in the instance of the rotating channel flow is about 390,000, and the total number of meshes in the instance of the rotating flow sudden expansion pipe is about 1,660,000. After introducing BSTMQ in the embodiment of the present invention, a transient calculation is performed to ensure that the maximum coolant number is less than 1 by the time step. In order to quantitatively compare the advantages of the BSTMQ model, the attached Figures 5 to 9 are the results of different variables in the two instances.

[0062] Figure 5 shows the average velocity distribution along the normal direction of the rotating channel calculated by various models. The average velocity calculated by the BSTMQ model is closest to the result of direct numerical simulation (DNS), which can predict the velocity distribution of the asymmetric distribution. The results of the traditional SSTk-ω model and STSST model are significantly different from the DNS results. Figure 6 shows the root-mean-square distribution of the average pulsating velocity along the flow direction of the rotating channel flow calculated by various models, and near the suction surface (y / h=2) where the rotation effect becomes more significant, the calculation results of the BSTMQ model are much better than those of the SSTk-ω model and the STSST model. Figure 7 shows the ratio of the turbulent viscosity nut to the molecular viscosity nu of the rotating channel flow calculated by various models, which can be used to show the degree of analysis of the flow field by the models. In the conventional SSTk-ω model, the nut / nu value is over 10 in most areas of the flow passage, and the analysis level is close to that of the RANS model. In the STSST model, the nut / nu values ​​are all over 10 near the suction surface, and the analysis level is close to that of the RANS model. In other locations, the nut / nu value is less than 10, and the analysis level is close to that of the large eddy simulation. In the BSTMQ model, the Nut / nu value of the flow passage is less than 1, and the analysis level is close to that of the DNS. Figure 8 compares the average axial velocity at Z / D=0.25 for the instance of the rotating sudden expansion. Figure 9 compares the root mean square of the circumferential average velocity pulsation at Z / D=4 for the instance of the rotating sudden expansion. In the case of the rotating sudden expansion, the calculation results of the BSTMQ model were also superior to those of the SSTk-ω model and the STSST model.< / volscalarfield> < / volscalarfield>

Claims

1. A method for predicting rotating turbulence using a BSTMQ turbulence model, comprising the steps of: (1) performing solid modeling of a prediction target of rotating turbulent flow using three-dimensional modeling software to obtain a simulation model; (2) meshing the simulation model; A step of performing rotation curvature correction on the existing STRUCT turbulence model based on the bifurcation method to construct a BSTMQ turbulence model, the equation of the BSTMQ turbulence model is as follows: [0051] where k is the turbulent kinetic energy and ω is the specific dissipation rate. [0052] is the partial derivative of the turbulent kinetic energy with respect to time, [0,53] is the partial derivative of the turbulent kinetic energy in the coordinate direction, [0,54] is the partial derivative of the specific dissipation rate with respect to time, [0.55] is the partial derivative of the specific dissipation rate in the coordinate direction, and u ｊ is the velocity, and C ＊ μ is the equation for the eddy viscosity coefficient, and C μ is the eddy viscosity, S is the mode of the strain rate tensor, v is the dynamic viscosity, and v Ｔ is the eddy viscosity coefficient, P ｋ is the turbulent kinetic energy production term, and D ｆ－ｂ is the eddy viscosity damping function, and б ｋ , B ω , α, β, F １ , B ω２ is a constant coefficient, step (3); (4) performing a numerical simulation calculation of the rotating turbulent flow using the BSTMQ turbulence model by a calculation program based on the fluid analysis software to obtain a prediction result of the rotating turbulent flow to be predicted; In the step (3), the step of constructing the BSTMQ turbulence model specifically includes the steps of: Step (31): Introduce the strain rate tensor and the rotation rate tensor in the rotating coordinate system; Step (32): Unify the rotation and curvature effects using the Spaart-Shur tensor to obtain a corrected rotation rate tensor; Step (33): Create a velocity gradient invariant from the strain rate tensor introduced in step (31) and the corrected rotation rate tensor obtained in step (32); Step (34): determining a corrected analysis time scale from the velocity gradient invariant obtained in step (33); Step (35): determining a modeling time scale for rotating turbulence from the specific dissipation rate; Step (36): Construct an eddy viscosity damping function D f−b that reflects the rotation and curvature effects according to steps (34) and (35); [0,56] Where: [0.57] and [0,58] is the dimensionless velocity gradient invariant III, [0,59] is the dimensionless velocity gradient invariant I constructed based on the strain rate tensor, [0.60] is the dimensionless velocity gradient invariant II constructed based on the corrected rotational velocity tensor, α is an empirical constant for the modeling time scale, Step (37): A method for predicting rotating turbulence, comprising: using the BSkO model as a base model, correcting the turbulent viscosity, and constructing a STRUCT turbulence model, i.e., a BSTMQ turbulence model, which is corrected based on the rotational curvature of the bifurcation method according to the eddy viscosity damping function D f-b .

2. In the inertial coordinate system, the strain rate tensor S ｉｊ and the rotational velocity tensor Ω ｉｊ The calculation formula is as follows: [0061] Here, i and j are two free indices in the tensor representation, both of which can take on values ​​1, 2, and 3, which means that we need to traverse all values ​​of the index, U ｉ and U ｊ is a velocity vector, and the three velocity components in the inertial coordinate system are U １ = U, U ２ = V, U ３ = W, and x ｉ and x ｊ are the coordinate axes in the inertial coordinate system, and x １ = x, x ２ = y, x ３ = z, [0062] are the partial derivatives of the three velocity components in the three directions, Angular velocity vector Ω Ｆ ｋ In a coordinate system rotating with , the rotational velocity tensor Ω Ａ ｉｊ The calculation formula is as follows: [0063] 【number】

3. When the rotation effect and the curvature effect are unified based on the Spalart-Shur tensor, the following corrected rotation rate tensor Ω ｍｏｄ ｉｊ The formula for calculation is obtained, [0,64] Here, C ｒ is a model constant, and W Ａ ｉｊ 3. The method for predicting a rotating turbulent flow according to claim 2, wherein: is a correction variable based on the Spalart-Shur tensor.

4. Correction variable W Ａ ｉｊ The decision In the case of two-dimensional flow, W Ａ ｉｊ =Ω ＳＳ ｉｊ And, 【number】 [0.65] 【number】 [0.66] Here, δ ｉｊ is the Kronecker symbol tensor, and S ｉｋ , and S ｋｊ , S ｉｊ are strain rate tensors with different indices, [0.67] 4. The method for predicting a rotating turbulent flow according to claim 3, wherein: 【Request 5】 【Number 68】 The calculation formula is as follows: [0,69] where k is the turbulent kinetic energy, ε is the dissipation rate, S is the mode of the strain rate tensor, and Ω ｍｏｄ 2. The method for predicting rotating turbulence according to claim 1, wherein: is the mode of the corrected rotational velocity tensor.

6. The corrected analysis time scale is: [Number 70] Where: [Equation 71] and Modeling rotating turbulence from specific dissipation rate time scale t ｍ Determine, [Equation 72] where α is an empirical constant for the modeled time scale, ω is the specific dissipation rate, [73] 6. The method of claim 5, wherein: represents performing a geometric mean operation in the flow region of interest.

7. Conventional eddy viscosity damping function D ｆ Based on this, the eddy viscosity damping function D reflecting the rotation effect and the curvature effect is ｆ－ｂ and the conventional eddy viscosity damping function D ｆ teeth, [74] 7. The method for predicting a rotating turbulent flow according to claim 6, wherein:

8. A computing device including a memory, a processor, and a computer program stored in the memory and operable on the processor, A computer device characterized in that said processor, when executing said computer program, implements the steps of the method according to any one of claims 1 to 7.

9. A computer-readable storage medium having a computer program stored thereon, the computer program being characterized in that, when executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 7.

Images