A weakly compressible aeroacoustic numerical solution method based on Roe's Riemann solver

Abstract

This invention provides a numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculations, comprising the following steps: establishing a sound field computational domain; obtaining the ground-state flow field of the weakly compressible flow; mapping the ground-state flow field to the sound field computational domain; interpolating the center values ​​of the cells in the sound field computational domain to the interface between the cell and adjacent cells; calculating the mean value at the interface using Roe averaging; calculating the flux at the interface; and solving the acoustic control equations established using linearized Euler equations. This invention can be widely applied in the field of aerodynamic noise calculation technology.

Description

Technical Field

[0001] This invention relates to the field of acoustic computing technology, and in particular to the field of aerodynamic noise calculation technology for turbomachinery with complex fluid flow. Background Technology

[0002] In engineering fields such as aviation, shipbuilding, and energy equipment, bladed rotating machinery (such as axial flow fans, compressors, propellers, turbines, and propellers) is widely used in pneumatic conveying and power conversion. With the continuous increase in equipment power density and rotational speed, aerodynamic noise has become one of the key bottlenecks affecting overall machine performance, environmental compliance, and comfort. Especially in fan and compressor systems, the rotating fundamental frequency noise (BPF) and its harmonics generated when blades pass through unsteady flow fields will undergo multiple reflections and couplings within complex pipe, casing, and hull structures, forming a broadband, strong noise field. Aerodynamic noise not only affects the equipment's operating environment but also leads to structural fatigue and reduced efficiency.

[0003] Traditional aerodynamic noise prediction methods mainly rely on semi-empirical formulas or acoustic analogy models (such as the Ffowcs Williams–Hawkings equations) under the incompressible Navier–Stokes approximation. While these methods can estimate far-field sound pressure levels under low Mach number conditions, they struggle to accurately describe the refraction, reflection, and attenuation of sound waves in complex background flows, and they are highly dependent on the unsteadiness of rotating machinery flows and the treatment of acoustic boundaries. On the other hand, while fully compressible solvers can capture complete wave characteristics, their computational cost is extremely high, making them difficult to apply on an engineering scale. Summary of the Invention

[0004] To address the issues of low accuracy and high computational cost in existing aerodynamic noise prediction methods, this invention provides a numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculations.

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

[0006] A numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculations includes the following steps:

[0007] S1. Establish the sound field calculation domain;

[0008] S2. Obtain the original term data of the weakly compressible flow field;

[0009] S3. Map the original data of the weakly compressible flow field to the grid of the acoustic field calculation domain;

[0010] Its characteristic is that it also includes the following steps:

[0011] S4. Interpolate the center value of the cell of the mesh of the sound field calculation domain to the interface between the cell and the adjacent cell;

[0012] S5. Calculate the flux on the interface using Roe average;

[0013] S6. Calculate the flux on the interface;

[0014] S7. Solve the acoustic control equations established using the linearized Euler equations.

[0015] Optionally, the interpolation is performed using a BJ limiter in step S4.

[0016] Optionally, the interpolation expression in step S4 includes:

[0017]

[0018]

[0019] in, For the total variation decreasing slope limiter, The coordinates of the center of the interface are... The coordinates of the center of the unit; and These are the disturbance amounts on the two surfaces of the interface, respectively; and These are the disturbance amounts at the centers of the units on both sides of the interface; and These are the perturbation gradients along the normal direction at the center of the units on both sides of the interface.

[0020] Optionally, the expression for calculating the mean on the interface in step S5 includes:

[0021]

[0022] in, These are the densities of the units on both sides of the interface; The average density of Roe; These are the velocity components in the x-direction of the units on both sides of the interface; These are the velocity components in the y-direction of the units on both sides of the interface; These are the velocity components in the z-direction of the units on both sides of the interface; The velocity component in the x-direction after Roe's average; The velocity component in the y-direction after Roe's average; The z-direction velocity component after Roe's average; This is the velocity vector after Roe's average. These are the total enthalpy of the units on both sides of the interface; This is the total enthalpy after Roe's average. Specific heat ratio; The pressure is derived from the average state of Roe; The speed of sound corresponding to the average state of Roe; Let be the unit normal vector of the interface; The normal velocity component of the Roe average velocity on the interface.

[0023] Optionally, the formula for calculating the flux at the interface in step S6 is as follows:

[0024]

[0025] in, These are the fluxes on both sides of the interface. These are the state perturbation vectors on both sides of the interface.

[0026]

[0027] The variable marked with an upper right apostrophe indicates that the variable is an unsteady disturbance. , These represent the pressure disturbance amounts on both sides of the interface.

[0028] Optionally, the acoustic control equations in step S7 are:

[0029]

[0030] in: This refers to the gradient term in the average flow term; ; , For the velocity perturbation in the x-direction, For the velocity perturbation in the y-direction, The velocity perturbation is in the z-direction. This represents the pressure disturbance in the z-direction.

[0031] Optionally, the following steps may be included between step S5 and step S6:

[0032] S51. Determine whether the cells of the grid in the current sound field calculation domain are sound field boundaries;

[0033] If step S51 determines that the element of the current sound field calculation domain's grid is not a sound field boundary, then step S6 is executed; if step S51 determines that the element of the current sound field calculation domain's grid is a sound field boundary, then the flux is calculated using the NRBC condition and boundary damping is applied.

[0034] Optionally, the formula for calculating the flux under the NRBC condition is:

[0035] When the sound field boundary is a non-reflective boundary, the flux formula is:

[0036]

[0037] When the sound field boundary is a total reflection boundary, the flux formula is:

[0038]

[0039] in, The numerical flux in the normal direction of the interface; To average the flow state in the background Below, by the disturbance state quantity The calculated normal physical flux; and These are the left and right disturbance states obtained by interface reconstruction of the left and right units on the interface, respectively. The background average flow state at the interface is used for linearization and flux assessment.

[0040] Optionally, the numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculations further includes the following steps after step S7:

[0041] S8. Determine whether the maximum Coulomb number of the unit is greater than a preset value. The expression:

[0042]

[0043] Among them, the acoustic Courant number The expression:

[0044]

[0045] in, For time step, The volume of the unit. The average velocity at the center of the unit. Let be the surface normal vector of the interface. Let be the area of ​​the interface. For reference speed of sound.

[0046] Optionally, if the maximum Courant number of the unit is less than or equal to the preset value, then the following steps are further performed:

[0047] S9. Determine whether the termination time has been reached;

[0048] If the result of step S9 is negative, then steps S5 to S9 are repeated.

[0049] The technical effects of this invention are as follows:

[0050] The technical solution of this invention, based on the linearized Euler equation (LEE) and combined with the Roe-averaged Riemann flux, achieves high-precision, low-dissipation acoustic simulation of the sound field propagation process. Common LEE solutions often employ a central scheme with artificial dissipation or a simple upwind scheme, which is prone to excessive dissipation or dispersion errors in weak amplitude sound fields. This method directly uses the "Roe average + characteristic matrix" for LEE flux calculation and provides an achievable surface flux formula and amplitude limiting control on an unstructured mesh, explicitly aligning with the "characteristic direction of weakly compressible sound waves," resulting in a better balance between accuracy and stability. In summary, the technical solution of this invention achieves the objective of this invention.

[0051] The further effects of the above-mentioned alternative methods will be explained in detail below with reference to specific implementation methods. Attached Figure Description

[0052] Figure 1 This is a flowchart of one embodiment of the method of the present invention.

[0053] Figure 2 This is a comparison chart showing the results of different methods used to process one-dimensional convection equations.

[0054] Figure 3 This is a sound pressure cloud map of a point sound source obtained by processing a shear flow background using the method of the present invention.

[0055] Figure 4 This is a comparison chart showing the accuracy of numerical simulation results of a point sound source under shear flow conditions with the results processed by the method of this invention.

[0056] Figure 5 This is a schematic diagram of the computational domain for a two-dimensional Gaussian pulse propagation example (open boundary).

[0057] Figure 6 The pressure pulsation diagram and result comparison diagram for a two-dimensional Gaussian pulse propagation example (open boundary) at t=0.0001 seconds are shown.

[0058] Figure 7 The pressure pulsation diagram and result comparison diagram for a two-dimensional Gaussian pulse propagation example (open boundary) at t=0.0003 seconds are shown.

[0059] Figure 8 The pressure pulsation diagram and result comparison diagram are shown for a two-dimensional Gaussian pulse propagation example (open boundary) at t=0.0006 seconds.

[0060] Figure 9The diagram shows the pressure pulsation and results comparison for a two-dimensional Gaussian pulse propagation example (open boundary) at t=0.001 seconds.

[0061] Figure 10 A schematic diagram of the computational domain for a two-dimensional Gaussian pulse propagation example (with boundaries set).

[0062] Figure 11 Pressure pulsation diagram and result comparison diagram for a two-dimensional Gaussian pulse propagation example (with boundaries set) at t=0.0004 seconds.

[0063] Figure 12 Pressure pulsation diagram and result comparison diagram for a two-dimensional Gaussian pulse propagation example (with boundaries set) at t=0.0009 seconds.

[0064] Figure 13 Pressure pulsation diagram and result comparison diagram for a two-dimensional Gaussian pulse propagation example (with boundaries set) at t=0.0015 seconds.

[0065] Figure 14 This is a schematic diagram showing the setup of the inclined plate model.

[0066] Figure 15 The flow field calculation contour plot is for the inclined plate model.

[0067] Figure 16 The results are the acoustic calculations for the inclined plate model. Detailed Implementation

[0068] The technical solution of the present invention will be described in detail below.

[0069] The linearized Euler equations can correctly explain convection effects, refraction in inhomogeneous media and inhomogeneous flows. Its main assumption is that instantaneous variables (density ρ, velocity u, and pressure p) can be decomposed into time-averaged quantities (variables with a right-hand subscript of 0) and small unsteady disturbances (variables with a right-hand superscript of ').

[0070] (1)

[0071] To obtain the system with wave variables, substitute equation (1) into the Euler equation (the effect of viscosity on wave propagation can be ignored), and subtract the steady-state mean flow equation. To linearize the system, the second-order terms (multiplication of first-order waves) can be ignored, as they are much smaller than first-order waves. Thus, the Euler equation is obtained in the following form:

[0072]

[0073] Where γ is the specific heat ratio of the fluid. In the case of weakly compressible flow, only isentropic flow is considered, i.e.:

[0074] (5)

[0075] Under the premise of equation (5), the density equation can be ignored, there is no large density jump, the pressure can be normalized relative to the density, equation (3) can be simplified to equation (6), and equation (4) can be simplified to equation (7):

[0076] (6)

[0077] (7)

[0078] The left-hand side of equations (6) and (7) is a linear form that includes all propagation effects. Generally, the system is excited by the nonlinear source terms on the right-hand side (i.e., monopole sources in the energy equation and dipole and quadrupole sources in the momentum equation). For data from a weakly compressible flow field, the following momentum source term exists:

[0079] (8)

[0080] Equations (6) and (7) are defined as follows:

[0081] (9)

[0082] In equation (9), the state vector is defined as , The gradient term in the average flow term is continuously updated. For the sound source term, the combination of A, B, and C forms the direction matrix of the Jacobian matrix A, defined as:

[0083] (10)

[0084] Let x be the normal vector of the plane in the x-direction. Let be the normal vector of the plane in the y-direction. Let be the normal vector of the z-direction plane. For the finite volume method, the governing equations are integrated over an arbitrary volume Ω, where the constant coefficient matrix is ​​extracted in the derivative using the divergence theorem, while the pressure and velocity in the variable coefficient matrix are piecewise constants within each cell. The solution to the Riemann problem explains the discontinuity of the average variable, thus yielding the following integral relationship:

[0085] (11)

[0086] Source Item Defined as:

[0087] (12)

[0088] Riemannian solutions using the complete Jacobian matrix on unstructured grids for:

[0089] (13)

[0090] Specifically, the background state variables are constructed using the Roe average (represented by an overline above the variable): It is given by averaging the square roots of the left and right states. For weighted velocity average, The weighted average of total enthalpy; based on this, the average pressure is defined. With the square of the average speed of sound Furthermore, the normal unit vector n indicates the interface orientation. This is the velocity vector after Roe's average. This represents the component of the average velocity in the normal direction, expressed by the following fractions:

[0091]

[0092] The eigenvalues ​​of the Jacobian matrix can be derived from... Find:

[0093] (twenty three)

[0094] Take two tangential vectors that are orthogonal to n and unitarily orthogonal: ,satisfy The left and right characteristic matrices of matrix A are obtained as follows:

[0095] (twenty four)

[0096] (25)

[0097] According to Obtain the matrix :

[0098] (26)

[0099] in, .

[0100] Since the Jacobian matrix is ​​the only input to the flux, the total flux through the cell boundaries (i.e., interfaces with adjacent cells) of the mesh in the acoustic field computation domain is given by the integral over all interfaces:

[0101] (27)

[0102] The surface integral is approximated using the midpoint, and the source term and time derivative are approximated using the cell center value. The time progression of each cell adopts a four-stage explicit Runge-Kutta scheme.

[0103] The flux at the interface between two adjacent cells is solved using the Riemann problem. For the variable coefficient matrix, the Roe method is employed.

[0104] (28)

[0105] in, For the flux in the left state, For the flux in the right state, respectively Let be the left-state perturbation vector. Let be the right-state perturbation vector, and its expressions are as follows:

[0106]

[0107] To minimize the diffusion of the solver implemented according to this method, higher-order extensions are required. This is accomplished by extrapolating the state variables from the center of each element to the interface in a finite linear manner:

[0108] (32)

[0109] (33)

[0110] Ψ represents the total variation decreasing slope limiter. When the gradient changes smoothly, Ψ = 1, and the method is second-order; when there is a significant gradient in the fluid, 0 < Ψ < 1, and the method is first-order. In one embodiment of the invention, a Barth-Jespersen type slope limiter (BJ limiter) is used.

[0111] To avoid non-physical reflection of sound waves by the outer boundary of the finite computation domain, this invention employs a non-reflection boundary condition based on eigenvalue decomposition, building upon the linearized Euler equations. The LEE is then placed in the boundary normal direction... Direction written in one-dimensional form:

[0112] (34)

[0113] in Let be the normal flux Jacobian matrix. For Perform eigenvalue decomposition:

[0114] (35)

[0115] The characteristic variables propagating along the normal direction are obtained:

[0116] (36)

[0117] in It contains two acoustic features And two convection characteristics .

[0118] For open boundaries, the feature quantities pointing to the interior of the computational domain will be... Set to zero, retaining only the characteristics of outward propagation. This approximates reflection-free absorption. Numerically, the flux is decomposed along the surface normal using the Roe eigenma matrix, retaining only the outward characteristic contribution. Therefore, the following approach is adopted at the boundary:

[0119] (37)

[0120] in Indicates all that satisfy The "incident" characteristic quantity, These are reference values ​​obtained by extrapolation from the internal solution; Let be the boundary reflection coefficient, and take . This corresponds to an ideal non-reflective boundary. This corresponds to a hard-walled boundary with total internal reflection.

[0121] No reflection boundary, i.e. This yields the commonly used Robin form, which can be considered a first-order approximation of the characteristic non-reflective boundary condition of this invention in a simple homogeneous medium, denoted as... ,but:

[0122] (38)

[0123] Based on this, a boundary virtual state is constructed. Then, with the inner surface value Substitute them together into the one-sided conserved Roe flux:

[0124] (39)

[0125] For a hard-walled boundary with total reflection, apply This corresponds to total sound pressure reflection and a zero normal velocity. In terms of flux realization, we take... but

[0126] (40)

[0127] To measure the stability conditions in acoustic numerical simulations, the acoustic Courant number is introduced. For mesh element i in the computational domain, the acoustic Courant number is... The definition of is:

[0128] (41)

[0129] in, For time step, For unit volume, The average velocity at the center of the unit. and These are the surface normal vector and the area, respectively. The reference speed of sound. Equation (41) takes into account the propagation effect of the convection term and the speed of sound term in the normal direction.

[0130] The global acoustic Courant number can be characterized by the maximum and average acoustic Courant numbers of the individual cells across the entire computational domain.

[0131] (42)

[0132] In equation (42), N represents the total number of grid cells in the computational domain. Typically, the stability requirements of numerical simulations... It does not exceed a certain critical value.

[0133] Figure 1 This paper illustrates a specific flowchart of an example of the numerical solution method for weakly compressible aeroacoustics based on the Roe Riemannian algorithm of this invention. The following section details... Figure 1 The specific process shown will be explained.

[0134] S1. Establish the sound field computation domain

[0135] In this step, a computational domain model of the sound field is built using SolidWorks, and imported into Workbench for mesh generation and boundary naming. The binary mesh data is then imported into OpenFoam to complete initialization. This step is implemented using existing technologies and will not be elaborated further here.

[0136] S2. Obtain the original terms data of the weakly compressible flow field.

[0137] In this step, a computational domain model of the flow field is established using SolidWorks and imported into Workbench for mesh generation and boundary naming. The binary mesh data is then imported into OpenFoam, and the rhoPimplefoam (OpenFoam's built-in weakly compressible flow field solver) is used to solve the flow field, thereby obtaining the original terms data of the weakly compressible flow field. The implementation of this step falls within the existing technical field and will not be elaborated further here.

[0138] S3, Mapping

[0139] In this step, the ground-state flow field is mapped to the acoustic field computation domain, which is achieved using existing technologies.

[0140] S4, Interpolation

[0141] In this step, the center values ​​of the cells in the sound field computational domain are interpolated to the interface between the cell and its adjacent cells. The expression used for interpolation is:

[0142]

[0143]

[0144] in, For the total variation decreasing slope limiter, The coordinates of the center of the interface are... The coordinates of the center of the unit; and These are the disturbance amounts on the two surfaces of the interface, respectively; and These are the disturbance amounts at the centers of the units on both sides of the interface; and These are the perturbation gradients along the normal direction at the center of the units on both sides of the interface.

[0145] S5. Calculate the mean value at the interface.

[0146] In this step, the mean value on the interface is calculated using the Roe average. The expression for the mean value includes:

[0147]

[0148] in, These are the densities of the units on both sides of the interface; The average density of Roe; These are the velocity components in the x-direction of the units on both sides of the interface; These are the velocity components in the y-direction of the units on both sides of the interface; These are the velocity components in the z-direction of the units on both sides of the interface; The velocity component in the x-direction after Roe's average; The velocity component in the y-direction after Roe's average; The z-direction velocity component after Roe's average; This is the velocity vector after Roe's average. These are the total enthalpy of the units on both sides of the interface; This is the total enthalpy after Roe's average. Specific heat ratio; The pressure is derived from the average state of Roe; The speed of sound corresponding to the average state of Roe; Let be the unit normal vector of the interface; The normal velocity component of the Roe average velocity on the interface.

[0149] S51, Is it a boundary?

[0150] In this step, it is determined whether the currently processed unit is the boundary of the far field of the sound field. If it is not the boundary of the far field of the sound field, then proceed to step S6. If it is determined that the currently processed unit is the boundary, then the "boundary processing" step is executed.

[0151] Boundary processing

[0152] In this step, the flux on the interface is calculated based on the following determination of the boundary type.

[0153] If the boundary is an open boundary, then the flux formula on the interface is:

[0154]

[0155] If the boundary is a hard-walled boundary with total internal reflection, then the flux formula on the interface is:

[0156]

[0157] in, The numerical flux in the normal direction of the interface; To average the flow state in the background Below, by the disturbance state quantity The calculated normal physical flux; and These are the left and right disturbance states obtained by interface reconstruction of the left and right units on the interface, respectively. The background average flow state at the interface is used for linearization and flux assessment.

[0158] S6. Calculate the flux at the interface.

[0159] In this step, the flux at the interface is calculated. The formula for calculating the flux at the interface is:

[0160]

[0161] in, These represent the flux in the left and right states, respectively. These are the left and right state perturbation vectors, respectively.

[0162]

[0163] S7. Solve the acoustic control equations

[0164] In this step, the following acoustic control equations are solved:

[0165]

[0166] in: This refers to the gradient term in the average flow term; ; , For the velocity perturbation in the x-direction, For the velocity perturbation in the y-direction, The velocity perturbation is in the z-direction. This represents the pressure disturbance in the z-direction.

[0167]

[0168]

[0169] The variable marked with an upper right apostrophe indicates that the variable is an unsteady disturbance. , These represent the pressure disturbance amounts on both sides of the interface.

[0170] S8. Does it satisfy the Courant number?

[0171] In this step, the following judgment is made on the current unit:

[0172] If the maximum Coulomb number of the unit is less than or equal to the preset value, then step S9 is executed; if the maximum Coulomb number of the unit is greater than the preset value, then an alarm message indicating an error is output (output error).

[0173] The maximum Courant number of the unit The expression:

[0174]

[0175] Among them, the acoustic Courant number The expression:

[0176]

[0177] in, For time step, The volume of the unit. The average velocity at the center of the unit. Let be the surface normal vector of the interface. Let be the area of ​​the interface. For reference speed of sound.

[0178] S9. Has the end time been reached?

[0179] In this step, if the termination time is determined to have been reached, proceed to the next step "output result"; if the termination time is determined not to have been reached, repeat steps S4 to S9.

[0180] Output

[0181] Utilize existing software to output results. For example, use Openfoam to output result files.

[0182] The technical effectiveness of the present invention is verified using multiple verification examples. The verification examples are described below.

[0183] Example 1: One-dimensional convection equation

[0184] One-dimensional convection equations are a common verification example for high-precision discrete schemes. Figure 2 A comparison is presented between the results obtained using the method of this invention (Roe–LEE) and the numerical results of the one-dimensional convection equation. The numerical results of the one-dimensional convection equation include the central difference baseline and the analytical solution. Figure 2 The vertical axis represents the disturbance variable. .

[0185] like Figure 2 As shown, the Roe–LEE curve and the analytical solution almost coincide at the wavefront and main peak, with no discernible peak shift. This indicates that the phase velocity error of the interfacial flux is effectively suppressed under the support of the face-centered Roe average and uniform eigendecomposition. In contrast, the central difference exhibits significant phase lag and peak misalignment upstream of the main peak, indicating that this scheme displays typical dispersion error in the absence of an upwind stabilizing term.

[0186] like Figure 2 As shown, the Roe-LEE solution exhibits minimal amplitude attenuation of the main peak, with no overshoot or ringing on either side of the peak, and the numerical solution maintains the same monotonicity as the analytical solution. In contrast, the central difference solution generates a series of oscillations upstream of the main peak, while the main peak is significantly weakened, indicating that its high-frequency pseudo-modes are not effectively suppressed. This comparison clearly demonstrates that the Roe-LEE solver can suppress cross-plane phase jumps and spurious attenuation, and reduce the accumulation of numerical dispersion during propagation, without introducing excessive dissipation.

[0187] This example demonstrates that the method of the present invention has both small phase and amplitude errors under the same grid and time step conditions, and its numerical stability and fidelity are better than those of the central difference method. It provides a reliable reference for the preservation of phase, pointing map and spectral shape in subsequent two-dimensional / three-dimensional non-uniform medium propagation and engineering geometry examples.

[0188] Example 2: Point Sound Source in Shear Flow Background

[0189] This example uses a point sound source within a shear flow. Shear flow is a typical form of non-uniform flow, where the velocity distribution is typically hyperbolic tangential along a certain direction, significantly affecting the propagation and refraction characteristics of sound waves. This example not only verifies the adaptability of the invention to complex background flows but also evaluates the stability of Roe averaging and eigendecomposition under non-uniform flow field conditions.

[0190] like Figure 3 and Figure 4 As shown, shear flow significantly affects the propagation of sound waves from a point source. Sound waves propagate faster on the high-velocity side and slower on the low-velocity side, resulting in wavefront deflection. Compared to the example with a point source without a flow field, the wavefront shape gradually shifts from concentric circles to an ellipse. The corresponding one-dimensional tangent profile results show that the numerical solution and the reference solution are highly consistent in terms of peak position and phase change, especially with minimal amplitude error in the main propagation direction, exhibiting only a slight attenuation deviation in the high-gradient region of the far field. This result further verifies the high-resolution characteristics of Roe Riemann averaging and eigenvalue decomposition flux calculations, as well as the stability and low dissipation of fourth-order time-progression in non-uniform flow fields.

[0191] Example 3: Two-dimensional Gaussian pulse propagation (open boundary)

[0192] Two-dimensional Gaussian pulse propagation is a common verification example in computational aeroacoustics. For example... Figure 5 As shown, there is an artificial sound source, namely the Gaussian initial pulse, at the center (0,0,0) of the square computational domain, and the Gaussian pulse propagates in a uniform flow field with a velocity of 100 m / s. Figure 5 In the middle, the top, bottom, left and right edges are all far field. In order to avoid the reflected waves generated after the sound waves reach the far field boundary and "pollute" the sound field by returning to the computational domain, these four edges all adopt the characteristic NRBC non-reflection boundary condition, and a damping sponge layer is activated at a distance of 0.1m from the boundary.

[0193] like Figure 6 , Figure 7 As shown, the pulse loop propagates outward in a concentric circle and shifts to the right along with the average flow. The numerical solution completely coincides with the analytical solution in terms of wavefront position and peak amplitude. This indicates that the Roe-LEE method effectively suppresses the interface phase velocity error and does not exhibit the upstream ringing and phase lag common in center difference classification methods. It can obtain the amplitude and phase of a high-fidelity two-dimensional propagation field.

[0194] like Figure 8 , Figure 9As shown, the outward-expanding pulse loop begins to enter the damping sponge layer region near the boundary. Subsequently, a systematic deviation from the analytical solution appears in the cross-sectional comparison, with a decrease in peak amplitude. This difference is not due to algorithmic dispersion, but rather the expected effect of the combined absorption layer-boundary processing. Specifically, the sponge layer dissipates incident energy within the body, and the boundary's characteristic conditions suppress reflection, thus preventing reflected waves from the outer edge of the sound field from returning and contaminating the computational domain. In real engineering calculations, the sound field boundary is usually placed further away, ensuring that the central computational domain does not interfere with the absorption layer during the simulation period. In this case, the absorption layer does not affect the solution within the domain; its only function is to block the return flow of wave energy. Within the computation time window of this example, no secondary wave systems caused by boundary reflection were observed entering the computational domain, indicating that the method combining the NRBC boundary and the damping sponge layer effectively suppresses wave reflection at the boundary.

[0195] Example 4: Two-dimensional Gaussian pulse propagation (with boundary settings)

[0196] This example is based on example 3, such as Figure 10 As shown, the four boundaries of the square computational domain are set as walls, and the total reflection hard wall boundary condition is adopted. The damping sponge layer is turned off and the artificial sound source is also a Gaussian pulse, which is placed at the position (0,0.2,0) in the computational domain. The background flow field is a uniform flow field with an incoming flow velocity of 100m / s, and the line segment AB is the set of sampling points.

[0197] like Figures 11 to 13 As shown, the numerical results and analytical solutions completely coincide in the three stages. At 0.0004s, the boundary is touched, and the numerical profile almost perfectly matches the analytical solution in terms of peak position and amplitude, indicating that the Roe-LEE method can accurately capture the phase velocity and peak value of the incident wave. Secondly, at 0.0009s, the pulse reaches another wall and is reflected, and the curve shows a parallel structure of the incident and reflected peaks. The phase relationship between the peaks matches the analytical solution, and the pressure amplitude near the wall is enhanced. Finally, at 0.0015s, the reflected wave reaches the third wall, and the two sets of peaks and valleys are consistent with the analytical solution. No pseudo-ringing or non-physical energy gain was observed, indicating that the Roe-LEE method can maintain energy balance and spectral consistency even in rigid wall reflection interaction scenarios.

[0198] Example 5: Aeroacoustic Radiation Based on the Wake of an Inclined Plate

[0199] This example selects a typical inclined plate wake radiation case for verification, which can verify the applicability of the method of the present invention in complex unsteady flow and wall interference acoustic fields. The inclined plate is set as follows: Figure 14 As shown, the angle of attack is 25° and the chord length is 0.18 m. Figure 15 The obtained flow field cloud map. Using the method of this invention... Figure 15The result represents the sound source term input into the acoustic calculation domain. The acoustic calculation results obtained using the method of this invention are as follows: Figure 16 As shown, the main noise source of the inclined plate wake is located in the shear layer vortex shedding region. The sound waves radiate downstream and outward from the trailing edge, exhibiting typical dipole radiation characteristics. The sound field forms a local high sound pressure region near the trailing edge of the inclined plate and spreads outward to form a ring-shaped wavefront structure. The calculation results are consistent with the acoustic distribution trends reported in existing experiments and literature, indicating that the method of the present invention can accurately capture broadband aerodynamic noise induced by unsteady flow.

[0200] This example verifies the stability and applicability of the invention under complex geometry, non-uniform flow, and near-wall interference sound field conditions, proving that it can not only be used for ideal sound source verification, but also directly applied to engineering aerodynamic noise problems such as wind turbine blades, compressor channels, and cabin walls, providing a reliable numerical means for noise prediction and noise reduction design of rotating machinery.

[0201] 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, all equivalent changes made based on the description and figures of the present invention, or direct or indirect applications to other related technical fields, are included within the scope of the present invention.

Claims

1. A numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculations, comprising the following steps: S1. Establish the sound field calculation domain; S2. Obtain the original term data of the weakly compressible flow field; S3. Map the original data of the weakly compressible flow field to the grid of the acoustic field calculation domain; Its characteristic is that it also includes the following steps: S4. Interpolate the center value of the cell of the mesh of the sound field calculation domain to the interface between the cell and the adjacent cell; S5. Calculate the flux on the interface using Roe average; S6. Calculate the flux on the interface; S7. Solve the acoustic control equations established using the linearized Euler equations.

2. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 1, characterized in that: In step S4, the interpolation is performed using a BJ limiter.

3. The numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculation as described in claim 2, characterized in that: The interpolation expression in step S4 includes: in, For the total variation decreasing slope limiter, The coordinates of the center of the interface are... The coordinates of the center of the unit; and These are the disturbance amounts on the two surfaces of the interface, respectively; and These are the disturbance amounts at the centers of the units on both sides of the interface; and These are the perturbation gradients along the normal direction at the center of the units on both sides of the interface.

4. The numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculation as described in claim 3, characterized in that: The expression for calculating the mean on the interface in step S5 includes: in, These are the densities of the units on both sides of the interface; The average density of Roe; These are the velocity components in the x-direction of the units on both sides of the interface; These are the velocity components in the y-direction of the units on both sides of the interface; These are the velocity components in the z-direction of the units on both sides of the interface; The velocity component in the x-direction after Roe's average; The velocity component in the y-direction after Roe's average; The z-direction velocity component after Roe's average; This is the velocity vector after Roe's average. These are the total enthalpy of the units on both sides of the interface; This is the total enthalpy after Roe's average. Specific heat ratio; The pressure is derived from the average state of Roe; The speed of sound corresponding to the average state of Roe; Let be the unit normal vector of the interface; The normal velocity component of the Roe average velocity on the interface.

5. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 4, characterized in that: The formula for calculating the flux at the interface in step S6 is as follows: in, These are the fluxes on both sides of the interface. These are the state perturbation vectors on both sides of the interface. The variable marked with an upper right apostrophe indicates that the variable is an unsteady disturbance. , These represent the pressure disturbance amounts on both sides of the interface.

6. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 5, characterized in that: The acoustic control equations described in step S7 are as follows: in: This refers to the gradient term in the average flow term; ; , For the velocity perturbation in the x-direction, For the velocity perturbation in the y-direction, The velocity perturbation is in the z-direction. This represents the pressure disturbance in the z-direction.

7. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 1, characterized in that: The following steps are also included between step S5 and step S6: S51. Determine whether the cells of the grid in the current sound field calculation domain are sound field boundaries; If step S51 determines that the element of the current sound field calculation domain's grid is not a sound field boundary, then step S6 is executed; if step S51 determines that the element of the current sound field calculation domain's grid is a sound field boundary, then the flux is calculated using the NRBC condition and boundary damping is applied.

8. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 7, characterized in that: The formula for calculating the flux under the NRBC condition is as follows: When the sound field boundary is a non-reflective boundary, the flux formula is: When the sound field boundary is a total reflection boundary, the flux formula is: in, The numerical flux in the normal direction of the interface; To average the flow state in the background Below, by the disturbance state quantity The calculated normal physical flux; and These are the left and right disturbance states obtained by interface reconstruction of the left and right units on the interface, respectively. The background average flow state at the interface is used for linearization and flux assessment.

9. The numerical solution method for weakly compressible aeroacoustics based on Roe Riemannian calculation as described in claim 1, characterized in that: The steps following step S7 are also included: S8. Determine whether the maximum Coulomb number of the unit is greater than a preset value. The expression: Among them, the acoustic Courant number The expression: in, For time step, The volume of the unit. The average velocity at the center of the unit. Let be the surface normal vector of the interface. Let be the area of ​​the interface. For reference speed of sound.

10. The numerical solution method for weakly compressible aeroacoustic based on Roe Riemannian calculation as described in claim 9, characterized in that: If the maximum Courant number of the unit is less than or equal to the preset value, then the following steps are further performed: S9. Determine whether the termination time has been reached; If the result of step S9 is negative, then steps S5 to S9 are repeated.

Images