A modeling method of turbulent large eddy simulation considering combustion effect

Abstract

The application discloses a turbulent large eddy simulation modeling method considering combustion effect, which comprises the following steps: firstly, obtaining basic flow field velocity distribution data; calculating multi-point velocity gradient and considering turbulent anisotropy caused by flame interface, normalizing modeling of large eddy turbulent sub-grid stress characteristic direction tensor; then, considering local pressure expansion effect caused by combustion, modeling and solving sub-grid kinetic energy of large eddy sub-grid kinetic energy transport equation; finally, combining normalized sub-grid characteristic direction tensor and sub-grid kinetic energy, calculating sub-grid momentum source term and coupling into the momentum equation solving process, so as to improve the accuracy of large eddy turbulent simulation. The modeling method fully considers the effects of local anisotropy and pressure expansion caused by combustion, and is more suitable for simulation research on engine and gas turbine combustion chamber, and can significantly improve the simulation accuracy of reactive flow.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology of turbulent flow, specifically to a turbulent large eddy simulation modeling method that considers combustion effects, applicable to industrial simulation calculation scenarios where flame and turbulence are coupled in the combustion chamber of an engine or gas turbine. Background Technology

[0002] Traditional aero-engine and gas turbine combustor development systems rely primarily on testing to expose design flaws, resulting in long development cycles, high costs, and high risks. However, with the advancement of computational fluid dynamics (CFD) technology, numerical simulations of turbulent combustion processes within the combustor using CFD methods can effectively support forward engineering, enabling rapid design iteration and optimization across all stages of the entire lifecycle. This significantly reduces the amount of physical testing, shortens design and testing time, lowers development costs, and improves development efficiency.

[0003] In turbulent CFD simulations, commonly used simulation methods include three typical frameworks: Direct Numerical Simulation (DNS), Large Eddy Simulation (LES), and Reynolds-Averaged Method (RANS). The DNS method does not introduce any model assumptions, thus achieving high computational accuracy; however, it consumes significant computational resources, making it difficult to apply to real-world industrial scenarios. The LES method uses some model assumptions, resulting in lower computational costs than the DNS method while maintaining a high level of accuracy. The RANS method boasts the highest computational efficiency but the lowest simulation accuracy.

[0004] In current developments in turbulence simulation, turbulence models are largely explored based on the RANS framework. Patent No. 202110408435.X discloses a turbulence length-scale correction method for the k-epsion turbulence model based on the RANS framework, and Patent No. 202110408433.0 discloses a length-scale correction method for the SST turbulence model based on the RANS framework. In these methods, the inventors correct the magnitude of the source term in the transport equation using dimensionless velocity divergence and introduce control functions to control the region of influence of the source term, avoiding dependence on wall distance parameters and improving the prediction effect on the aerodynamic forces and aerothermal effects of hypersonic vehicles. Patent application number 201911309994.4 discloses a turbulence model for numerical simulation of vortex-induced vibration of a rigid cylinder. Under the RANS framework, and with reference to the constraint method of the SST k-ω model, the low Reynolds number k-ε turbulence model is modified, and a turbulence model applicable to the numerical simulation of vortex-induced vibration of a rigid cylinder in the subcritical Reynolds number range is established.

[0005] However, the RANS method has inherent limitations. Although it can be corrected through model adjustments, its computational accuracy remains low, and it can only provide steady-state flow field information, failing to obtain flow field information that evolves over time. Patent application number 202210480149.9 discloses an adaptive turbulence simulation method suitable for afterburner flow and combustion processes. This method achieves adaptive transition between unsteady RANS, LES, and DNS modes by calculating and adaptively adjusting the resolution control function, thus achieving adaptive turbulence simulation and improving the simulation accuracy of turbulent flow while maintaining a certain level of computational efficiency.

[0006] However, none of the above methods consider the combustion effect in the turbulent flow modeling process. Combustion leads to significant pressure expansion, which further causes anisotropy of the fluid along the flame surface normal within the flame space. Therefore, it is necessary to develop a large eddy simulation modeling method for turbulent flow that considers the combustion effect, improve the simulation accuracy of reactive turbulence in engine and gas turbine combustors, support the forward development of combustors, reduce development costs, and improve development efficiency. Summary of the Invention

[0007] The purpose of this invention is to address the shortcomings of existing technologies by providing a turbulent large eddy simulation modeling method that considers combustion effects, thereby solving the problem of large deviations in the simulation calculation results of reactive turbulence inside the combustion chamber in existing technologies.

[0008] The objective of this invention is achieved through the following technical solution:

[0009] A turbulent large eddy simulation modeling method considering combustion effects includes the following steps:

[0010] (1) Begin the calculation of the nth iteration step, and collect basic information data in the flow field (the calculation of the nth iteration step starts from the data information collection, and after collection, it enters the calculation process of each intermediate quantity); the basic information includes physical space coordinates x, y, z, velocity vector field U, pressure p, density ρ, scalar Sc isc＝1,n ;

[0011] (2) Based on the basic information data in step (1), the reaction kinetics are solved and calculated to obtain the thermodynamic state of the new iteration step, the scalar transport equation is solved and the flow field density information is updated.

[0012] (3) Normalize and model the subgrid feature direction tensor using the basic flow field information data from step (1); for each grid cell [x l y m , z nFirst, the velocity gradients on the local grid cell [p,q=0] and its surrounding grid cells [p,q=-1~1] are solved, and the overall combustion process variable C is calculated. Then, based on C, the flame association factor θ between the surrounding grid cells and the local grid cell in the flame space is determined. p,q Then, based on the structured modeling idea of ​​Taylor expansion, the velocity gradients on the local grid cell and surrounding grid cells are sampled at multiple points based on the flame correlation factor to obtain the subgrid feature direction tensor G of the local grid cell. ij Finally, calculate G. ij / G mm To G ij Perform normalized modeling;

[0013] (4) Based on the momentum equation, a subgrid kinetic energy transfer equation was derived, taking into account the expansion effect caused by combustion. Terms in the equation that could not be directly solved by large eddy simulation were modeled. Based on the basic flow field data information from step (1) and G... ij The transport equations were solved to obtain the subgrid kinetic energy k. sgs ;

[0014] (5) Based on the subgrid normalized feature tensor direction G in step (3) ij and the subgrid kinetic energy k in step (4) sgs Reconstructing the subgrid turbulent stress τ ij ;

[0015] (6) Calculate the subgrid filter source term of the momentum equation and substitute the source term into the iterative calculation of the kinetic energy equation to complete the turbulent momentum solution process considering the combustion effect, and realize the high-precision simulation study of the turbulent flame flow field in the combustion chamber of the engine and gas turbine.

[0016] (7) Judge the calculation result of sub-iteration step k. If it has not converged, repeat steps (2) to (6) to continue the calculation of the (k+1)th sub-iteration step. If it has converged, proceed to the calculation of the (n+1)th iteration step until the target time or the total number of iterations is reached.

[0017] Furthermore, the sub-lattice feature orientation tensor G described in step (3) ij The velocity gradient effect of the target mesh cell and its surrounding mesh is considered, and its expression is:

[0018]

[0019] Where i, j, k represent spatial geometric directions, u i u j Refers to the velocity components in the i and j directions, x k Refers to the geometric coordinate information in the k-direction.

[0020] In particular, α represents the spatial index of the local or adjacent grid cell, where p1p2 = p2p3 = p1p3 = q1q2 = q2q3 = q1q3 = 0, and α and β cannot simultaneously represent different adjacent grid cells. l, m, and n are the indices of the local grid cell in the geometric directions x1, x2, and x3, respectively. p and q refer to the grid cell information based on the local grid cell index, distributed in -1, 0, and 1. Only one of p1, p2, and p3 can be non-zero at the same time, and the same applies to q. For example, when p2 = 1, then p1 = p3 = 0, and α represents the adjacent grid cell with an index greater than 1 in the x2 direction; when p3 = -1, then p1 = p2 = 0, and α represents the adjacent grid cell with an index less than 1 in the x3 direction; and when p1 = p2 = p3 = 0, α represents the local grid cell. Δ α,k The length of the grid cell at position α in space along the k direction. C α,β θ represents the weighting coefficient of the velocity gradient between adjacent grid cells. α,β This represents the flame association factor between adjacent grid cells.

[0021] Furthermore, the adjacent grid cell velocity gradient weighting coefficient C α,β Defined as:

[0022]

[0023] Furthermore, the flame correlation factor θ of the adjacent grid cells α,β Defined as:

[0024]

[0025] In the formula θ represents the vector direction from the local grid cell to the adjacent grid cell (α or β), while θ represents the gradient vector of the flame process variables. and The angle between them.

[0026] Furthermore, the subgrid kinetic energy parameter equations are modeled and solved. The modeling method for the transport equations is as follows:

[0027]

[0028]

[0029] Furthermore, in the formula, the overline represents spatial filtering, and the cap represents density-weighted spatial filtering. Refers to the sub-lattice kinetic energy convection transport term. The term ε refers to the subgrid kinetic energy pressure expansion term, and ε refers to the subgrid kinetic energy dissipation term. Refers to the subscale flux term of subgrid kinetic energy. Refers to the subgrid kinetic stress transport term. The term refers to the three-dimensional unclosed term of the sublattice kinetic energy. In the formula... ε represents the Taylor expansion coefficients. The modeling of the subgrid kinetic energy dissipation term ε employs the scaling approximation approach, assuming a proportional relationship between subgrid kinetic energy and subgrid dissipation, with a correlation coefficient γ = 0.5. Where k... rs Large eddy scale kinetic energy, S ij Let be the stress tensor, expressed as σ ij =2μ(S) ij -1 / 3S kk δ ij μ is the kinematic viscosity, δ ij For Kroneck's symbol.

[0030] Furthermore, in steps (5) and (6), the normalized sub-lattice feature orientation tensor G is used. ij and subgrid kinetic energy k sgs Reconstructing the subgrid stress tensor τ ij The subgrid momentum source term SrcU is calculated in conjunction with the solution process of the precession momentum equation. The algebraic expression is as follows:

[0031]

[0032]

[0033] The beneficial effects of this invention are as follows: Addressing the shortcoming of current industrial combustion flow simulation programs that do not consider the significant impact of combustion on flow, this invention considers the thermal expansion and anisotropic effects of combustion during the modeling process, which helps improve the simulation accuracy of reactive turbulence. This invention can be applied to the numerical simulation of flow fields within industrial combustion chambers such as engine combustion chambers and gas turbine combustion chambers, supporting the forward development of combustion chambers, reducing development costs, and improving development efficiency. Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the model's operation steps;

[0035] Figure 2 A schematic diagram of the combustion chamber structure in the model, which includes a swirling premixed flame; Figure 2 In the middle, D refers to the diameter of the combustion chamber, and L0 refers to the length of the combustion chamber inlet section;

[0036] Figure 3 This is a verification diagram comparing the correlation coefficients between the subgrid feature orientation tensor and the target tensor; each subgroup in the diagram represents τ from left to right. 11 τ 22 τ 33 τ 12 τ 23τ 13 ;

[0037] Figure 4 A comparative verification diagram of subgrid kinetic energy and its dissipation;

[0038] Figure 5 A comparison diagram of the cross-sectional field of heat release rate obtained from the model calculation;

[0039] Figure 6 A comparison chart of the average and pulsating velocities calculated by the model. Detailed Implementation

[0040] The present invention will be further described and illustrated below with reference to the accompanying drawings and specific embodiments.

[0041] This invention primarily constructs a turbulent large eddy simulation modeling method considering combustion effects, and simulates and analyzes the turbulence in the combustion chamber environment coupled with combustion heat release effects. The specific steps are as follows:

[0042] Step (1) begins the calculation of the nth iteration step, collecting basic information data in the flow field; the basic information includes physical space coordinates x, y, z, velocity vector field U, pressure p, density ρ, and scalar Sc. isc＝1,n ;

[0043] Step (2) Based on the basic information data in step (1), the reaction kinetics are solved and calculated to obtain the thermodynamic state of the new iteration step, the scalar transport equation is solved and the flow field density information is updated;

[0044] Step (3) uses the basic flow field information data from step (1) to normalize and model the subgrid feature direction tensor; for each grid cell [x l y m , z n First, the velocity gradients on the local grid cell [p,q=0] and its surrounding grid cells [p,q=-1~1] are solved, and the overall combustion process variable C is calculated. Then, based on C, the flame association factor θ between the surrounding grid cells and the local grid cell in the flame space is determined. p,q Then, based on the structured modeling idea of ​​Taylor expansion, the velocity gradients on the local grid cell and surrounding grid cells are sampled at multiple points based on the flame correlation factor to obtain the subgrid feature direction tensor G of the local grid cell. ij Finally, calculate G. ij / G mm To G ij Perform normalized modeling;

[0045] Step (4) derives the subgrid kinetic energy transfer equation based on the momentum equation, considering the expansion effect caused by combustion, and models the terms in the equation that cannot be directly solved by large eddy simulation. Based on the basic flow field data information from step (1) and G... ij The transport equations were solved to obtain the subgrid kinetic energy k. sgs ;

[0046] Step (5) is based on the sub-lattice normalized feature tensor direction G in step (3). ij and the subgrid kinetic energy k in step (4) sgs Reconstructing the subgrid turbulent stress τ ij ;

[0047] Step (6) Calculate the subgrid filter source term of the momentum equation and substitute the source term into the iterative calculation of the kinetic energy equation to complete the turbulent momentum solution process considering the combustion effect, and realize the high-precision simulation study of the turbulent flame flow field in the combustion chamber of the engine and gas turbine.

[0048] Step (7) judges the calculation result of sub-iteration step k. If it has not converged, repeat steps (2) to (6) to continue the calculation of the (k+1)th sub-iteration step. If it has converged, proceed to the calculation of the (n+1)th iteration step until the target time or the total number of iterations is reached.

[0049] The sub-lattice feature orientation tensor G mentioned in step (3) ij The velocity gradient effect of the target mesh cell and its surrounding mesh is considered, and its expression is:

[0050]

[0051] Where i, j, k represent spatial geometric directions, u i u j Refers to the velocity components in the i and j directions, x k Refers to the geometric coordinate information in the k-direction.

[0052] In particular, α represents the spatial index of the local or adjacent grid cell, where p1p2 = p2p3 = p1p3 = q1q2 = q2q3 = q1q3 = 0, and α and β cannot simultaneously represent different adjacent grid cells. l, m, and n are the indices of the local grid cell in the geometric directions x1, x2, and x3, respectively. p and q refer to the grid cell information based on the local grid cell index, distributed in -1, 0, and 1. Only one of p1, p2, and p3 can be non-zero at the same time, and the same applies to q. For example, when p2 = 1, then p1 = p3 = 0, and α represents the adjacent grid cell with an index greater than 1 in the x2 direction; when p3 = -1, then p1 = p2 = 0, and α represents the adjacent grid cell with an index less than 1 in the x3 direction; and when p1 = p2 = p3 = 0, α represents the local grid cell. Δ α,k The length of the grid cell at position α in space along the k direction. C α,β θ represents the weighting coefficient of the velocity gradient between adjacent grid cells. α,β This represents the flame association factor between adjacent grid cells.

[0053] The adjacent grid cell velocity gradient weighting coefficient C α,β Defined as:

[0054]

[0055] The flame correlation factor θ of adjacent grid cells α,β Defined as:

[0056]

[0057] In the formula θ represents the vector direction from the local grid cell to the adjacent grid cell (α or β), while θ represents the gradient vector of the flame process variables. and The angle between them.

[0058] The subgrid kinetic energy parameter equations are modeled and solved. The modeling method for the transport equations is as follows:

[0059]

[0060] In the formula, the overline represents spatial filtering, and the cap represents density-weighted spatial filtering. Refers to the sub-lattice kinetic energy convection transport term. The term ε refers to the subgrid kinetic energy pressure expansion term, and ε refers to the subgrid kinetic energy dissipation term. Refers to the subscale flux term of subgrid kinetic energy. Refers to the subgrid kinetic stress transport term. The term refers to the three-dimensional unclosed term of the sublattice kinetic energy. In the formula... ε represents the Taylor expansion coefficients. The modeling of the subgrid kinetic energy dissipation term ε employs the scaling approximation approach, assuming a proportional relationship between subgrid kinetic energy and subgrid dissipation, with a correlation coefficient γ = 0.5. Where k... rs Large eddy scale kinetic energy, S ij Let be the stress tensor, expressed as σ ij =2μ(S) ij -1 / 3S kk δ ij μ is the kinematic viscosity, δ ij For Kroneck's symbol.

[0061] In steps (5) and (6), the normalized sub-lattice feature orientation tensor G is used. ij and subgrid kinetic energy k sgs Reconstructing the subgrid stress tensor τ ij The subgrid momentum source term SrcU is calculated in conjunction with the solution process of the precession momentum equation. The algebraic expression is as follows:

[0062]

[0063]

[0064] The specific implementation effect of the above method is demonstrated below with reference to the embodiments.

[0065] Example

[0066] In this embodiment, the direct numerical simulation results of the swirling premixed turbulent flame in the model combustion chamber are used as the target verification object. The model developed in this invention is used for simulation calculation within the framework of large eddy simulation, and the calculation results are compared with the direct numerical simulation results to illustrate the specific implementation effects and advantages of this invention.

[0067] In this embodiment, the model combustion chamber structure is as follows: Figure 2 As shown in the table below, the specific parameters of the swirling premixed turbulent flame are as follows:

[0068] Table 1. Swirl Premixed Flame Parameter Settings

[0069] pressure 20atm Inlet average axial velocity 40m / s Inlet temperature 760K swirl number 1.0 Equivalent ratio 0.6

[0070] The simulation calculations employ a low Mach number, weakly compressible finite difference program, with second-order precision in both time stepping and spatial discretization. The time stepping uses a semi-implicit Crank-Nicolson scheme, while the convection terms of the scalar transport equation and the subgrid kinetic energy transport equation are calculated using a third-order WENO scheme.

[0071] By comparing the calculated data from the model with the data from direct numerical simulation, Figure 3The correlation coefficients between the subgrid feature orientation tensor and the target tensor are compared and verified. The target tensor is taken from the spatially filtered data of direct numerical simulation. The results show that the large eddy simulation model (MCKGM) considering combustion effects described in this invention has a higher correlation coefficient than other traditional models (SSM, VGM, KGM) and provides more accurate predictions of the subgrid feature orientation tensor.

[0072] Figure 4 The paper presents a comparison between the subgrid kinetic energy calculated by the model and the subgrid kinetic energy obtained by direct numerical simulation with spatial filtering, and also compares the subgrid kinetic energy dissipation intensity. It can be found that the traditional model severely underestimates the magnitude of the subgrid kinetic energy, while the MCKGM model described in this invention is more accurate and also has high accuracy in simulating the subgrid kinetic energy dissipation intensity.

[0073] Figure 5 The results demonstrate the predictive effect of the model described in this invention on the intensity of combustion heat release, proving that the model described in this invention can accurately reproduce the results of direct numerical simulation.

[0074] Figure 6 The comparison and verification diagram of the average and fluctuating velocities calculated by the model is shown. The results prove that the model described in this invention has better performance than other traditional models for both the average velocity and fluctuating velocity of the flow field, especially in the prediction of fluctuating velocity, where the accuracy is extremely high.

[0075] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A method for modeling turbulent large eddy simulation considering combustion effects, characterized in that, Includes the following steps: 1) Begin the calculation for the nth iteration step, collecting basic information data in the flow field; the basic information includes physical space coordinates x, y, z, velocity vector field U, pressure p, density ρ, and scalar Sc. isc=1,n ; 2) Based on the basic information data in step 1), the reaction kinetics are solved to obtain the thermodynamic state of the new iteration step, the scalar transport equation is solved and the flow field density information is updated; 3) Use the basic flow field information data from step 1) to normalize and model the subgrid characteristic direction tensor; for each grid cell [x l y m , z n First, the velocity gradients on the local grid cell [p,q = 0] and its surrounding grid cells [p,q = -1~1] are solved, and the overall combustion process variable C is calculated. Then, based on C, the flame association factor between the surrounding grid cells and the local grid cell in the flame space is determined. Then, based on Taylor expansion-based structured modeling, the velocity gradients on the local grid cell and surrounding grid cells are sampled at multiple points using the flame correlation factor to obtain the subgrid feature orientation tensor G of the local grid cell. ij Finally, calculate G. ij / G mm To G ij Perform normalized modeling; where G mm =G 11 +G 22 +G 33 That is to say, G ij The sum of all terms in a tensor when i=j; 4) Based on the momentum equation, a subgrid kinetic energy transfer equation was derived, considering the expansion effect caused by combustion. Terms in the equation that could not be directly solved by large eddy simulation were modeled. Based on the basic information data from step 1) and the subgrid characteristic direction tensor G... ij The transport equations were solved to obtain the subgrid kinetic energy k. sgs ; The modeling method for the subgrid kinetic energy transfer equation is as follows: ； In the formula, the overline represents spatial filtering, and the cap represents density-weighted spatial filtering; Refers to the subgrid kinetic energy convection transport term. The term refers to the subgrid kinetic energy pressure expansion term. The term refers to the kinetic energy dissipation of the sublattice. Refers to the subscale flux term of subgrid kinetic energy. Refers to the subgrid kinetic stress transport term. The three unclosed terms of the sublattice kinetic energy are referred to; in the formula These are the Taylor expansion coefficients; Sub-grid kinetic energy dissipation term The modeling employs the scale approximation approach, assuming a proportional relationship between subgrid kinetic energy and subgrid dissipation, with a correlation coefficient γ = 0.5; where k rs For large eddy-scale kinetic energy, S ij Let be the stress tensor, denoted as , μ is the kinematic viscosity. The symbol for Kronecker. For average density, For average pressure, L mm The Leonard kinetic energy is approximately the same as the scale. Let x be the grid scale in the k direction. i , x j , x k Let i be the spatial coordinates of the i, j, and k directions; 5) Based on the sub-lattice normalized feature tensor direction G in step 3), ij and the subgrid kinetic energy k in step 4) sgs Reconstructing subgrid turbulent stress ; 6) Calculate the subgrid filter source term of the momentum equation and substitute the source term into the iterative calculation of the kinetic energy equation to complete the turbulent momentum solution process considering the combustion effect, and realize the high-precision simulation study of the turbulent flame flow field in the combustion chamber of engine and gas turbine. 7) Judge the calculation result of sub-iteration step k. If it has not converged, repeat steps 2) to 6) to continue the calculation of the (k+1)th sub-iteration step. If it has converged, proceed to the calculation of the (n+1)th iteration step until the target time or the total number of iterations is reached.

2. The turbulent large eddy simulation modeling method considering combustion effects as described in claim 1, characterized in that, The sublattice characteristic direction tensor G ij The velocity gradient effect of the target mesh cell and its surrounding mesh is considered, and its expression is: ； Where i, j, k represent spatial geometric directions, u i u j Refers to the velocity components in the i and j directions, x k Geometric coordinate information in the k-direction; , Indicates the spatial location index of the local or adjacent grid cell. Furthermore, α and β cannot simultaneously represent different adjacent grids; l, m, and n are the indices of the local grid cells in the geometric directions x1, x2, and x3; p and q refer to the grid cell information based on the local grid cell indices, distributed in -1, 0, and 1; only one of p1, p2, and p3 can be non-zero at the same time, and the same applies to q; when p2 = 1, then p1 = p3 = 0, and α represents the adjacent grid with an index greater than 1 in the x2 direction than the local grid cell; when p3 = -1, then p1 = p2 = 0, and α represents the adjacent grid with an index less than 1 in the x3 direction than the local grid cell; and when p1 = p2 = p3 = 0, α represents the local grid; This represents the length of the grid cell at position α in space along the k direction; This represents the weighting coefficient of the velocity gradient between adjacent grid cells. This represents the flame association factor between adjacent grid cells.

3. The turbulent large eddy simulation modeling method considering combustion effects as described in claim 2, characterized in that, The adjacent grid cell velocity gradient weighting coefficient Defined as: 。 4. The turbulent large eddy simulation modeling method considering combustion effects as described in claim 2, characterized in that, The adjacent grid cell flame correlation factor Defined as: ； In the formula θ represents the vector direction from the local grid cell to the adjacent grid cell, while θ represents the gradient vector of the flame process variables. and The angle between them.

5. The turbulent large eddy simulation modeling method considering combustion effects as described in claim 1, characterized in that, Based on the normalized sub-lattice feature orientation tensor G ij and subgrid kinetic energy k sgs Reconstructing the subgrid stress tensor The subgrid momentum source term SrcU is calculated using the coupled precession momentum equation solution process, and the expression is as follows: ； 。

Images