Low Mach number incompressible flow noise prediction method, storage medium and computer equipment

Abstract

The invention provides a low-Mach-number incompressible flow noise prediction method, a storage medium and computer equipment. The noise prediction method comprises the following steps: after a flow field to be evaluated reaches dynamic stability, performing transient numerical simulation on the flow field by using a low-speed incompressible flow solver to obtain a transient fluid pressure fluctuation value; using an average field solver to calculate an average pressure field and an average velocity field of at least two periods of the flow field after dynamic stabilization; solving an acoustic control equation after the time step flow calculation is finished, wherein the fluctuation pressure is obtained by the difference between the transient fluid pressure fluctuation value and the average pressure field; and then a boundary acoustic variable is updated by using a non-reflection boundary formula. The technical scheme of the invention can be used in the field of numerical prediction of complex engineering flow noise.

Description

Technical Field

[0001] The invention relates to the technical field of noise prediction, in particular to the technical field of noise prediction for low-speed incompressible flow. Background Art

[0002] When flowing gas passes through an obstacle, it will generate noise, which is called aerodynamic noise. Many devices, such as fans, cars, aircraft, etc., generally generate aerodynamic noise during operation. Aerodynamic noise can lead to a series of problems such as noise pollution, so it is necessary to evaluate the aerodynamic noise generated by the equipment.

[0003] People expect to predict the aerodynamic noise generated by the equipment before the equipment is produced, so as to avoid actual testing of aerodynamic noise with the manufactured equipment and save costs. Therefore, a variety of noise prediction methods have been developed. Computational Fluid Dynamics (CFD) is a type of software program that simulates and analyzes physical phenomena such as fluid flow, heat transfer and mass transfer through numerical calculation methods. It can simulate various complex flow conditions and can be used to predict the aerodynamic noise of equipment.

[0004] OpenFOAM is an open source CFD tool with high flexibility and scalability. The main calculation process of this tool is:

[0005] 1. Divide the calculated area into discrete "grids";

[0006] 2. Given physical information on the boundary grid, simulate the actual computing environment;

[0007] 3. Solve the governing equations (partial differential equations) that satisfy the physical phenomena on each small grid;

[0008] 4. Obtain the calculation results and visualize them.

[0009] In step 2 of the above calculation process, the perfectly matched layer (PML) and buffer zone solutions are often used. PML and buffer zone are methods used to simulate infinite extension areas. Their main function is to effectively absorb outward propagating waves, reduce reflection and scattering, and thus improve the accuracy of the simulation. However, in the implementation process of PML and buffer zone, additional computational grids need to be added to apply damping to achieve the effect of eliminating reflection. The additional grids will reduce the calculation speed and consume more computing resources. Summary of the invention

[0010] In order to solve the problem of consuming computing resources due to the existence of PML and Buffer zone in the prior art, the present invention provides a low Mach number incompressible flow noise prediction method, storage medium and computer equipment.

[0011] The technical solution of the present invention is as follows:

[0012] The low Mach number incompressible flow noise prediction method comprises the following steps:

[0013] A. When the flow field to be evaluated reaches dynamic stability, a low-speed incompressible flow solver is used to perform transient numerical simulation on the flow field to obtain transient fluid pressure fluctuation values;

[0014] B. calculating the average pressure field and the average velocity field of at least two cycles of the flow field after dynamic stability is achieved using a mean field solver;

[0015] C. Solving the acoustic control equation after the time step flow calculation is completed, and the fluctuating pressure is obtained by the difference between the transient fluid pressure fluctuation value and the average pressure field;

[0016] D. After completing step C, update the boundary acoustic variables using the non-reflecting boundary formula, where:

[0017]

[0018] Among them, p a To calculate the sound pressure, t is the time, u a To calculate the speed of sound, u a ,v a ,w a for u a The three components of λ4 are the characteristic propagation velocity, ξ is the calculation coordinate, c0 is the speed of sound, ρ0 is the ambient density, ξ x To calculate the partial derivatives of the coordinates with respect to the physical coordinates, v0 is the background flow velocity in the y direction, and w0 is the background flow velocity in the z direction.

[0019] Optionally, the low-speed incompressible flow solver in step A uses a turbulence model of IDDES based on the SST k-omega two-equation model, and the control equation of the turbulence model is:

[0020]

[0021] Among them, ρ is the fluid density, k is the turbulent kinetic energy, ω is the turbulent specific dissipation rate, u is the fluid flow velocity, μ is the fluid kinematic viscosity, and μ t is the turbulent viscosity, P k is the turbulent kinetic energy generation term, F1 and F2 are SST mixing functions, S is the strain rate tensor, σ k , σ ω, σ ω2 , α, β and a1 are model constants, and the dissipation term in the first equation of the turbulence model is corrected by introducing the mixing length scale:

[0022]

[0023] l LES =C DES Δ,

[0024]

[0025] Among them, l IDDES is the length scale of IDDES, Δ is the length scale of LES, l RANS and l LES is a hybrid function that takes wall modeling into account; is the empirical mixing function, f e is the lifting function, C DES and C μ is a model constant.

[0026] Optionally, the acoustic control equation in step C is:

[0027]

[0028] in, is the calculated average velocity field, p ic The pressure is fluctuating.

[0029] Optionally, after the damping term is introduced into the acoustic control equation in step C:

[0030]

[0031] Among them, the expression of damping coefficient is:

[0032] ν=(R v -1)c0dx / R Δ +0.01,

[0033] R v is the maximum volume growth rate between the current unit and the adjacent units, R Δ is 3.14, dx is the unit characteristic length.

[0034] Optionally, the time advancement numerical format K in the process of solving the acoustic control equation in step C is i :

[0035] K i =Δt F(R n +β i K i-1 )i=1,2,...m,

[0036] R n+1 =R n +K m ,

[0037] Where m = 4, R n ,R n+1 are the original variable vectors before and after solving the above acoustic control equations, Δt is the time step, F() represents the spatial discrete term, β i The expression is:

[0038]

[0039] A storage medium stores a computer program, which, when executed by a processor, implements the aforementioned low Mach number incompressible flow noise prediction method.

[0040] A computer device comprises: a memory, a processor and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the aforementioned low Mach number incompressible flow noise prediction method.

[0041] The technical effects of the present invention are as follows:

[0042] The low Mach number incompressible flow noise prediction method of the present invention adopts a non-reflecting boundary formula to update the boundary acoustic variables, and the effect is similar to that of the existing PML and Buffer zone methods. The non-reflecting boundary adopted is directly applied to the boundary of the calculation domain, and there is no need to use an additional calculation area like PML and Buffer zone, which reduces the number of grids required for calculation, improves the calculation efficiency, and achieves the purpose of the present invention.

[0043] The further effects of the above optional manner will be described below in conjunction with the specific implementation manner. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1 The figure is a flow chart of an embodiment of the noise prediction method of the present invention.

[0045] Figure 2 This is a diagram of the numerical clutter phenomenon.

[0046] Figure 3 This is the effect diagram after adopting numerical clutter suppression measures.

[0047] Figure 4 Schematic diagram of simulation experiment calculation.

[0048] Figure 5 Schematic diagram of the time variation of an oblique plane wave in a 2D pipe.

[0049] Figure 6 This is a comparison chart of the reflection intensity of the monitoring points. DETAILED DESCRIPTION

[0050] The technical solution of the present invention is described in detail below in conjunction with the accompanying drawings.

[0051] The low Mach number incompressible flow noise prediction method of the present invention is performed based on the OpenFOAM platform. Before executing the method of the present invention, a conventional meshing step is required. This step belongs to the existing technology and will not be repeated here.

[0052] like Figure 1 As shown, the low Mach number incompressible flow noise prediction method of the present invention includes the following steps.

[0053] Perform transient numerical simulation

[0054] In this step, when the flow field to be evaluated reaches dynamic stability, a low-speed incompressible flow solver is used to perform transient numerical simulation on the flow field to obtain transient fluid pressure fluctuation values. The low-speed incompressible flow solver is OpenFOAMpisoFoam, which belongs to the existing technology and will not be described in detail here.

[0055] In this step, the turbulence model used is based on the IDDES (improved delayed detached-eddy simulation) turbulence model of the SST (Shear Stress Transport) k-omega two-equation model. The control equation of the turbulence model is:

[0056]

[0057] Among them, ρ is the fluid density, k is the turbulent kinetic energy, ω is the turbulent specific dissipation rate, u is the fluid flow velocity, μ is the fluid kinematic viscosity, and μ t is the turbulent viscosity, P k is the turbulent kinetic energy generation term, F1 and F2 are SST mixing functions, S is the strain rate tensor, σ k , σ ω , σ ω2 , α, β and a1 are model constants, and the dissipation term in the first equation of the turbulence model is corrected by introducing the mixing length scale:

[0058]

[0059] l LES =C DES Δ,

[0060]

[0061] Among them, l IDDES is the length scale of IDDES, Δ is the length scale of LES (Large Eddy Simulation), l RANS and l LES is a mixed function that takes wall modeling into account, is the empirical mixing function, f e is the lifting function, C DES and C μ is a model constant.

[0062] Calculate the mean flow field

[0063] In this step, a mean field solver is used to calculate the mean pressure field and the mean velocity field of at least two cycles of the flow field after dynamic stability is achieved.

[0064] The specific method for implementing this step is illustrated by an example: define a mean field variable pMean and initialize it to 0, superimpose the transient fluid pressure fluctuation value obtained in each time step to pMean, and store the total number of time steps; when the specified time is reached to start solving the acoustic control equation, divide pMean by the number of time steps used for averaging to obtain the average value within a period of time.

[0065] Solving the governing equations for acoustics

[0066] In this step, the acoustic control equation is solved after the time step flow calculation is completed, and the fluctuating pressure is obtained by the difference between the transient fluid pressure fluctuation value and the average pressure field. The acoustic control equation is:

[0067]

[0068] where p a To calculate the sound pressure, u a To calculate the speed of sound, is the average velocity field calculated in step 2, p ic is the fluctuation pressure, ρ0 is the far-field environment density, and c0 is the far-field sound speed.

[0069] In this step, when solving the acoustic equations, the time-marching scheme uses the low-dispersion and low-dissipation Runge–Kutta format, which can reduce the loss during sound wave propagation and ensure the stability of the solution when Courant numbers C0 = c0Δt / Δx < 1. Numerical format K i as follows:

[0070] K i =Δt F(R n +β i K i-1 )i=1,2,...m,

[0071] R n+1 =R n +K m ,

[0072] Where m = 4, R n ,R n+1 are the original variable vectors before and after solving the acoustic control equations, Δt is the time step, F(·) represents the spatial discrete term, β i is the optimization parameter and can be expressed as:

[0073]

[0074] Since the numerical results of acoustics are easily affected by numerical errors, high-frequency clutter will be generated in areas with strong discontinuities. Therefore, numerical clutter will be generated in areas with large grid growth rates and at boundaries. In this step, artificial damping is added to eliminate numerical clutter. The damping term introduced into the acoustic control equation is transformed into the following formula:

[0075]

[0076] The expression of the damping coefficient is:

[0077] ν=(R v -1)c0dx / R Δ +0.01,

[0078] Where R v is the maximum volume growth rate between the current unit and the adjacent units, R Δ is a constant coefficient of 3.14, d x is the unit characteristic length.

[0079] Figure 2 The figure shows the phenomenon of numerical clutter when no artificial damping is added. As the Gaussian wave propagates outward, non-physical fluctuations occur in the local area with a large unit volume growth rate, and clutter appears before the sound wave reaches the boundary. Figure 3 As shown in the figure, the sound wave is still very stable when it reaches the boundary, without any clutter. After the sound wave propagates to the boundary of the computational domain, no non-physical phenomenon occurs, which effectively suppresses the generation of numerical clutter.

[0080] Update boundary acoustic variables

[0081] In this step, after solving the acoustic control equation, the boundary acoustic variables are updated using the non-reflecting boundary formula. The following is the derivation process of the non-reflecting boundary formula.

[0082] The acoustic control equations are solved using the step-wise convective wave equation (PCWE), which is written in a conservative form:

[0083]

[0084] Where U is the conservation state vector, F x ,F y ,F z is the flux vector in three directions. The specific form is as follows:

[0085] U={p a ,ρ0u a ,ρ0v a ,ρ0w a}, (2)

[0086]

[0087] where u a ,v a ,w a Represents u a The three components of represent The three components.

[0088] Define a mapping from physical coordinates to calculated coordinates:

[0089] ξ=ξ(x,y,z),η=η(x,y,z),ζ=ζ(x,y,z), (4)

[0090] Where (ξ,η,ζ) is the computational space coordinate, (x,y,z) is the physical space coordinate. To simplify the writing, the following symbols are defined:

[0091]

[0092] Through linear approximation, the conservation control equation in the physical coordinate system can be converted into a non-conservation equation in the computational coordinate system:

[0093]

[0094] The vector V is the original variable vector V = {p a ,u a ,v a ,w a}, the characteristic matrix B in the x direction ξ for:

[0095]

[0096] The characteristic matrices in the other two directions can also be obtained in the same way. ξDecomposing it, we can get:

[0097]

[0098] in yes The right eigenvector matrix of , the columns of which are right eigenvectors, is included Diagonal matrix of eigenvalues:

[0099]

[0100] The above formula defines the modified convection velocity Corrected Speed ​​of Sound The diagonal elements of are the propagation velocities of the characteristic waves, which are λ2=λ3=0, The physical meanings of the diagonal matrix are the sound waves propagating along the positive and negative ξ directions and the vortex waves generated by the hydrodynamic force. The specific expressions of the corrected convection velocity and the corrected sound velocity are:

[0101]

[0102] Considering the assumption of local one-dimensional inviscidity, ignoring the transverse spatial derivatives in the other two directions in formula (6), and using formula (8) to decompose the matrix in formula (6), we can obtain:

[0103]

[0104] Define two new variables:

[0105]

[0106] Formula (11) can be transformed into the following form:

[0107]

[0108] Among them, the time change rate of the characteristic amplitude value L ξ and the local spatial derivative term d ξ The component forms are:

[0109]

[0110] Substituting formula (14) into formula (13), we can obtain the complete local one-dimensional inviscid relation (LODI):

[0111]

[0112] In actual use, formula (15) will cause the calculation of acoustic variables at the boundary to diverge. The reason is that the tangential acoustic velocity at the boundary does not change with time and only depends on the initial value. However, it is impossible to directly give a suitable initial value for each simulation. During the program running, the boundary value is unstable, resulting in calculation crash. To overcome this problem, the last two equations in the formula are modified as follows:

[0113]

[0114] The physical meaning of formula (16) is that the transverse sound velocity is the sum of the sound velocity and the background flow velocity outflow boundary. This assumption can ensure that the transverse velocity does not cause adverse disturbances to the boundary and maintain numerical stability.

[0115] The final control equation at the boundary is:

[0116]

[0117] The set of equations (17) is the characteristic non-reflecting boundary (LMCBC) formulation applicable to the low Mach number flow noise governing equations (PCWE).

[0118] The following comparative experiments are conducted to compare the method of the present invention with the method using PML and Buffer zone type boundaries to verify the technical effect of the method of the present invention.

[0119] Case description:

[0120] The simulation domain is a 1m*3m two-dimensional pipeline, and the grid spacing is a constant 0.00625m, ensuring that the shortest wavelength studied can also contain 50 grid points. The number of grids fully meets the calculation requirements of the second-order finite volume. It can be considered that numerical dissipation and numerical dispersion will not have a significant impact on the research results. The calculation diagram is shown in Figure 4 The computational grid consists of three parts: region 1 is used to generate plane waves, region 2 is the propagation region, and region 3 is used to set the PML and Buffer zone (region 3 is removed when testing the LMCBC method of the present invention). The upper and lower surfaces of the computational domain are set as periodic boundaries. The performance of the three boundaries at different frequencies, wave angles, and Mach numbers is tested. The monitoring point at (2.3m, 0.5m) is used to evaluate the reflection intensity of the three boundaries.

[0121] Figure 5 A schematic diagram showing the time-varying propagation of an inclined plane wave in a 2D pipe shows that the acoustic calculation conforms to the laws of physics and that the sound wave propagates along the x-axis at the speed of sound. Figure 5 c It can be seen that the sound wave flows out smoothly after reaching the outlet boundary without obvious reflection.

[0122] Figure 6The reflection intensity monitored at the aforementioned monitoring points is shown, and the reflection intensity of three different boundaries at different Mach numbers Ma and different incident angles θ is shown. The reflection intensity of the LMCBC boundary proposed in the present invention is lower than that of the Buffer Zone at all tested Mach numbers and incident angles, achieving a complete improvement in the effect. Compared with PML, the reflection of LMCBC is slightly higher than that of PML when the incident angle is greater than 45°, and is better than PML at low incident angles, and the overall effect is similar. However, LMCBC does not require the exit area 3 during the calculation process, which reduces the number of grids and effectively improves the calculation efficiency.

[0123] The computer program stored on the storage medium of the embodiment of the present invention and the computer program executed by the computer device of the embodiment of the present invention can be executed as follows: Figure 1 An embodiment of the low Mach number incompressible flow noise prediction method of the present invention is shown.

[0124] The technical solution of the present invention can take into account the diffraction, reflection and absorption effects of rigid walls and structural materials on sound waves, realize the numerical prediction of complex engineering flow noise, and provide a numerical means for noise prediction and noise reduction design of related products.

[0125] It is worth noting that the above description is only a preferred embodiment of the present invention, and does not limit the scope of patent protection of the present invention. The present invention can also be replaced by equivalent technologies. Therefore, any equivalent changes made by using the description and illustrations of the present invention, or directly or indirectly applied to other related technical fields, are included in the scope of the present invention.

Claims

1. A method for predicting low Mach number incompressible flow noise, characterized by: The steps include: A. When the flow field to be evaluated reaches dynamic stability, a low-speed incompressible flow solver is used to perform transient numerical simulation on the flow field to obtain transient fluid pressure fluctuation values; B. calculating the average pressure field and the average velocity field of at least two cycles of the flow field after dynamic stability is achieved using a mean field solver; C. Solving the acoustic control equation after the time step flow calculation is completed, and the fluctuating pressure is obtained by the difference between the transient fluid pressure fluctuation value and the average pressure field; D. After completing step C, update the boundary acoustic variables using the non-reflecting boundary formula, where: Among them, p a To calculate the sound pressure, t is the time, u a To calculate the speed of sound, u a ,v a ,w a for u a The three components of λ4 are the characteristic propagation velocity, ξ is the calculation coordinate, c0 is the speed of sound, ρ0 is the ambient density, ξ x To calculate the partial derivatives of the coordinates with respect to the physical coordinates, v0 is the background flow velocity in the y direction, and w0 is the background flow velocity in the z direction.

2. The low Mach number incompressible flow noise prediction method according to claim 1, characterized in that: The low-speed incompressible flow solver in step A uses the turbulence model of IDDES based on the SST k-omega two-equation model. The governing equation of the turbulence model is: Among them, ρ is the fluid density, k is the turbulent kinetic energy, ω is the turbulent specific dissipation rate, u is the fluid flow velocity, μ is the fluid kinematic viscosity, and μ t is the turbulent viscosity, P k is the turbulent kinetic energy generation term, F1 and F2 are SST mixing functions, S is the strain rate tensor, σ k , σ ω , σ ω2 , α, β and a1 are model constants, and the dissipation term in the first equation of the turbulence model is corrected by introducing the mixing length scale: l LES =C DES D, Among them, l IDDES is the length scale of IDDES, Δ is the length scale of LES, l RANS and l LES is a hybrid function that takes wall modeling into account; is the empirical mixing function, f e is the lifting function, C DES and C μ is a model constant.

3. The low Mach number incompressible flow noise prediction method according to claim 2, characterized in that: The acoustic control equation described in step C is: in, is the calculated average velocity field, p ic The pressure is fluctuating.

4. The method for predicting low Mach number incompressible flow noise according to claim 3, characterized in that: After the damping term is introduced into the acoustic control equation described in step C: Among them, the expression of damping coefficient is: ν=(R v -1)c0dx / R Δ +0.01, R v is the maximum volume growth rate between the current unit and the adjacent units, R Δ is 3.14, dx is the unit characteristic length.

5. The low Mach number incompressible flow noise prediction method according to claim 4, characterized in that: The time-marching numerical format K in the process of solving the acoustic control equation in step C is i : K i =ΔtF(R n +β i K i-1 ) i=1,2,...m, R n+1 =R n +K m , Where m = 4, R n ,R n+1 are the original variable vectors before and after solving the above acoustic control equations, Δt is the time step, F( ) represents the spatial discrete term, β i The expression is:

6. A storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the low Mach number incompressible flow noise prediction method according to any one of claims 1 to 5 is implemented.

7. A computer device, characterized in that: include: A memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the low Mach number incompressible flow noise prediction method according to any one of claims 1 to 5.

Images