Unsteady flow field acceleration calculation method based on space-time hybrid matrix POD dimension reduction

Abstract

The invention discloses an unsteady flow field acceleration calculation method based on space-time hybrid matrix POD dimensionality reduction, which comprises the following steps: (1) in a dual-time iteration process, periodically collecting physical quantity data of an unsteady flow field at different time steps, and stacking the data at different time steps according to columns to form a unified space-time hybrid data matrix; (2) executing intrinsic orthogonal decomposition on the space-time mixed data matrix, and constructing a low-dimensional feature space; (3) dynamically updating the low-dimensional feature space; (4) reconstructing flow field data based on the low-dimensional feature space to serve as an initial value of next double-time iteration; and (5) in the internal iteration process, judging the acceleration effect of initial value improvement by calculating the residual error so as to reduce the calculation amount and accelerate the solution. According to the method, on the premise that the calculation precision is guaranteed, the number of iteration steps in calculation and the calculation cost of the unsteady flow field are reduced, the data dimension can be effectively reduced, the use of calculation resources is optimized, and the method has high practical application value.

Description

Technical Field

[0001] The present invention relates to an accelerated calculation method for an unsteady flow field based on time-space mixing matrix POD dimension reduction, which is suitable for fast calculation of an aircraft's unsteady flow field. Background Art

[0002] In the fields of modern engineering and scientific research, the simulation and calculation of unsteady flow fields face severe challenges of high computational complexity and long computational time consumption. Especially when it comes to large-scale flow field evolution and unsteady aerodynamic problems, traditional numerical methods often have extremely high computational costs and cannot meet the real-time and high efficiency requirements of engineering practice. Although the continuous advancement of computational fluid dynamics (CFD) technology has improved the computing speed, in practical applications, problems such as aerodynamic optimization, flow control and environmental simulation still require huge computing resources and long computing cycles. Therefore, how to meet the requirements of efficient flow field simulation in engineering practice and accelerate the calculation of unsteady flow fields has become an important research direction today.

[0003] Against this backdrop, Proper Orthogonal Decomposition (POD), a data-driven order reduction method, offers new possibilities for rapid flow field computation. By extracting the key characteristic patterns in flow data, POD reduces high-dimensional flow data to a lower-dimensional space, effectively reducing the computational effort while preserving key flow characteristics. The POD method significantly improves computational efficiency and has become a popular research topic for rapid flow prediction and optimization.

[0004] However, the existing POD method is still primarily used in steady flow fields, while its application in unsteady computations remains significantly limited. The traditional POD reduction process is typically based on pre-set modes, extracting modes from the statistical average state of the flow field or within a fixed time window. This approach is effective for steady flows, but for unsteady flows, the fixed modes cannot adapt to transient changes due to the time-dependent evolution of the flow field modes, resulting in large reconstruction errors and impacting computational stability and convergence. Furthermore, the computational model after POD reduction is often fixed and cannot adjust modal weights and energy distribution in real time, resulting in poor adaptability to complex unsteady flows. Although some studies have attempted to extend POD to unsteady flows, for example, Audouze et al. established a two-layer POD for unsteady variable parameter problems, constructing spatial and temporal basis functions with special properties so that the reduced-order model satisfies boundary and initial conditions. Cao et al. combined POD with an optimization algorithm, optimizing the POD basis coefficients through gradient optimization to reconstruct a flow field with small residuals, achieving a speedup of 2-3 times. Furthermore, Rapun et al. proposed a local POD method to accelerate the solution of parabolic unsteady problems. This method solves the equations by alternating between numerical schemes and reduced-order models, and determines the reduced-order solution time based on a priori estimates of the approximate error of the reduced-order model, thereby improving the robustness and efficiency of the solution process. In summary, the application of existing POD methods to unsteady flow calculations still has gaps or limitations. Summary of the Invention

[0005] The purpose of the present invention is to provide an accelerated calculation method for unsteady flow fields based on the POD dimension reduction of the spatiotemporal mixing matrix to save the cost of numerical calculation. The present invention constructs a reduced-order model through data-driven historical flow field snapshots and searches for flow fields with smaller residuals in the reduced-order space. Compared with the application of existing POD methods in flow fields, existing PODs are usually reduced based on preset modes and cannot fully adapt to the dynamic changes of the flow field over time, resulting in large reconstruction errors in unsteady flows, affecting the calculation stability and convergence speed. The present invention uses historical flow field snapshots to guide the dimensionality reduction process, breaking through the limitations of traditional POD static modes, dynamically updating and optimizing the reduced-order model so that it can reflect the latest state of the flow field at every moment, and adjusting the energy distribution and modal weights in real time according to historical flow field data; in addition, the present invention maximizes energy dimensionality reduction in the reduced-order space, giving priority to retaining the modes that contribute the most to the total disturbance energy, ensuring that the important dynamic characteristics of the flow field can be effectively retained while the computational efficiency is significantly improved; more importantly, the reconstructed flow field is used as the initial state for calculation at the next moment, so that each iteration can start from a more accurate starting point, thereby effectively accelerating the convergence process of the numerical solution and providing a new accelerated computing method for real-time simulation and prediction of unsteady flows.

[0006] The method and technical solutions proposed by the present invention are as follows:

[0007] A method for accelerating the calculation of unsteady flow fields based on POD dimensionality reduction of a spatiotemporal mixing matrix comprises the following steps:

[0008] Step (1) In the dual-time iteration process, the physical quantity data of the unsteady flow field at different time steps are periodically collected, and the data at different time steps are stacked in columns to form a unified time-space hybrid data matrix;

[0009] Step (2) performing proper orthogonal decomposition (POD) on the data matrix, sorting the eigenvectors according to the eigenvalues, prioritizing the modes that contribute the most to the total disturbance energy, determining the energy retention ratio, dynamically selecting the number of modes, and constructing a low-dimensional feature space;

[0010] Step (3) dynamically updates the low-dimensional feature space, i.e., repeatedly executing steps (1) and (2) within a preset iteration period to adapt to the characteristics of the flow field evolving over time;

[0011] Step (4) reconstructing the flow field data based on the low-dimensional feature space, and using the reconstructed flow field as the initial value of the next dual-time iteration;

[0012] In step (5), during the inner iteration process, the acceleration effect of the initial value improvement is judged by calculating the residual. When the inner iteration residual is lower than the preset threshold, some or all inner iteration steps of the current time step are skipped and the calculation of the next physical time step is directly entered to reduce the amount of calculation and speed up the solution.

[0013] The intrinsic orthogonal decomposition in step (2) uses the energy maximization principle to perform dimensionality reduction processing on the snapshot data, and the process includes:

[0014] a. Assume that the space-time mixed perturbation matrix Its construction method is as follows: in several consecutive time steps, disturbance data including physical quantities such as velocity, pressure, and density are collected, and the spatial field data at each moment are stacked in the form of column vectors to form a column vector sequence arranged in time, thereby forming a unified space-time mixed disturbance matrix, which is subjected to POD decomposition to obtain:

[0015]

[0016] Where x represents the position quantity, and the dimension changes with the dimension of the problem, such as the x coordinate in one dimension, the (x, y) coordinate in two dimensions, and the (x, y, z) coordinate in three dimensions; t represents the time quantity; m is the number of spatial grid points; N is the number of snapshots collected; B = [φ1,φ2,...,φ r ] is the orthogonal space mode set, S=diag(σ1,σ2,...,σ r) is a diagonal matrix of singular values ​​sorted by size, and V is the time coefficient matrix;

[0017] b. Sort the modes according to the contribution of disturbance kinetic energy. The energy contribution is determined by the singular value. Each mode φ i The contribution to the total energy is:

[0018]

[0019] Where r is the rank of the snapshot matrix, σ i For each mode, the corresponding singular value, E i The corresponding energy contribution for each mode.

[0020] c. Select the first N main modes from largest to smallest according to the cumulative energy contribution rate, and satisfy the following adaptive conditions:

[0021]

[0022] Wherein, η is the preset energy retention threshold, which is usually set to 0.95 or 0.98, and is used to ensure that the low-dimensional modal subspace retains the main energy of the disturbance as much as possible;

[0023] d. The selected main mode constructs the low-dimensional feature space B N =[φ1,φ2,...,φ N ], used for subsequent flow field reconstruction and accelerated calculation.

[0024] The low-dimensional feature space in step (3) adopts a dynamic update mechanism, that is, repeatedly performing snapshot acquisition and intrinsic orthogonal decomposition operations within a preset physical time step interval or iteration cycle, including:

[0025] a. During the flow field calculation process, the latest sets of snapshot data are collected again at set intervals (e.g., every several physical time steps) to form an updated perturbation matrix;

[0026] b. Perform singular value decomposition on the disturbance matrix to re-obtain eigenvalues ​​and modes;

[0027] c. Based on the updated energy spectrum, the main mode is filtered according to the preset energy threshold to construct a new low-dimensional feature space;

[0028] d. The updated low-dimensional feature space is used for flow field reconstruction and initial value generation in subsequent physical time steps, thereby adapting to the time-dependent evolution of the main modes in unsteady flows.

[0029] The residual threshold in step (5) is set to a fixed value determined according to the residual convergence trend of the original CFD iteration, preferably 10 -1 to 10 -3 , the specific value is adjusted according to the complexity of the problem.

[0030] The present invention has the following beneficial effects:

[0031] The proposed method applies Proper Orthogonal Decomposition (POD) to a compressible unsteady computational fluid dynamics (CFD) solver. Low-dimensional modeling is performed on flow field snapshot data. The approximate flow field for the next physical time step is predicted by selecting the main POD fundamental modes, which are used as the initial values ​​for the dual-time-step iterative process. Compared with traditional methods that use the solution of the previous time step, this method can effectively shorten the convergence time of the inner iteration, reducing the number of internal iterations by up to approximately 60%. Furthermore, during implementation, this method can dynamically adjust the number of POD modes based on a residual feedback mechanism, further achieving an adaptive balance between error control and computational acceleration. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is a schematic flow chart of the method of the present invention.

[0033] Figure 2 It is to select the model geometry and mesh schematic.

[0034] Figure 3 This is the comparison result between the convergence process of the flow field residual calculated by the subsonic example of the airfoil of the present invention and the original iterative method.

[0035] Figure 4(a)-Figure 4(d) The results of the comparison between the aerodynamic force and surface pressure calculated by the subsonic calculation example of the airfoil of the present invention and the original iterative method are shown;

[0036] Among them, Figure 4(a) shows the change of aerodynamic coefficient Cx with angle of attack; Figure 4(b) shows the change of aerodynamic coefficient Cy with angle of attack; Figure 4(c) shows the distribution of surface pressure in the x direction; Figure 4(d) shows the distribution of surface pressure in the y direction. Figure 5 It is the comparison result between the convergence process of the flow field residual calculated by the hypersonic example of the airfoil of the present invention and the original iterative method.

[0037] Figure 6(a)-Figure 6(d) The results are compared between the aerodynamic forces and surface pressures calculated by the hypersonic calculation example of the airfoil of the present invention and the original iterative method;

[0038] Among them, Figure 6(a) shows the change of aerodynamic coefficient Cx with angle of attack; Figure 6(b) shows the change of aerodynamic coefficient Cy with angle of attack; Figure 6(c) shows the distribution of surface pressure in the x direction; Figure 6(d) shows the distribution of surface pressure in the y direction. Figure 7 This is the comparison result between the flow field residual convergence process calculated by the low pitch speed example of the airfoil of the present invention and the original iterative method.

[0039] Figure 8(a)-Figure 8(d) The results of the comparison between the aerodynamic force and surface pressure calculated by the low pitch speed example of the airfoil of the present invention and the original iterative method are shown;

[0040] Among them, Figure 8(a) shows the change of aerodynamic coefficient Cx with angle of attack; Figure 8(b) shows the change of aerodynamic coefficient Cy with angle of attack; Figure 8(c) shows the distribution of surface pressure in the x direction; Figure 8(d) shows the distribution of surface pressure in the y direction. DETAILED DESCRIPTION

[0041] The specific embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0042] The present invention aims to effectively improve the computational efficiency of unsteady flow fields through the intrinsic orthogonal decomposition technology. By maximizing energy and reducing dimensionality, the original high-dimensional flow field data is compressed into a low-dimensional space, thereby accelerating the numerical simulation and analysis process. This method can improve computational efficiency and reduce computing resource consumption. It is suitable for large-scale flow field simulation and real-time prediction. The method process is as follows: Figure 1 shown.

[0043] First, the principle of the POD method is given:

[0044]

[0045] Where U(x,t) is a physical quantity, t is the time representation, and x is the space representation; is the average value; U'(x,t) is the pulsation value; a i (t) is the time coefficient related to time, is the spatial orthogonal basis function or POD spatial basis mode. The POD method is to take the first N modes with the highest energy and discard some unimportant information.

[0046]

[0047] in is the approximate deviation matrix.

[0048] Based on this principle, the following execution steps are determined in flow field calculation:

[0049] Step 1: Flow field snapshot acquisition

[0050] In the dual time-step iterative solution process, the state data of the unsteady flow field at different time points are periodically collected every several physical time steps. i , extract the disturbing physical quantities (such as velocity, pressure, density, etc.) in the flow field to form a spatial distribution vector U'(x,t i The spatial vectors collected at each time step are stacked in time order as column vectors to construct a unified space-time mixed disturbance matrix. x represents only the position quantity, and the dimension changes with the dimension of the problem, such as the x coordinate in one dimension, the (x, y) coordinate in two dimensions, and the (x, y, z) coordinate in three dimensions.

[0051]

[0052] in

[0053] Step 2: Flow field disturbance extraction and POD decomposition

[0054] After subtracting the mean, the variables are separated, converting the disturbance U'(x,t) into a set of relatively independent modes in space and time. When the number of samples reaches the preset requirement, proper orthogonal decomposition (POD) is performed. By performing eigenvalue decomposition on the covariance matrix of the flow field data, the flow field changes are decomposed into a set of spatial modes (spatial modes) and temporal modes (temporal modes) [B,S,V] = svd(U'), where B represents the spatial modal matrix, which represents the spatial distribution structure of the disturbance; S represents the singular value matrix, which contains all singular values. These values ​​are a measure of the correlation between the corresponding spatial modes and temporal modes, generally representing the contribution of each mode to the total change. The larger the singular value, the more important the mode; and V represents the temporal modal matrix, which reflects the evolution of each spatial mode over time.

[0055] Step 3: Main mode extraction

[0056] In unsteady flows, the main flow structures (such as shear layers, unstable modes, and main vortices) usually carry the vast majority of the energy, while small perturbations and numerical noise have very little energy. By retaining only high-energy modes, the main physical characteristics of the flow field can be preserved to a great extent, while compressing the data size and improving computational efficiency. By sorting the eigenvalues ​​and calculating the cumulative contribution rate, the main modes that meet the preset energy threshold (above 95%) are selected in descending order, the low-energy modes are truncated, and the number of modes is dynamically selected to form a dimensionality reduction representation such as formula (2).

[0057] Step 4: Flow field reconstruction and initial value update

[0058] The disturbance flow field is reconstructed using the selected main mode and time coefficient, and the approximate flow field solution is obtained by adding the time mean. The reconstructed flow field is used as the initial field of the new physical time step in the dual-time iteration process, replacing the traditional initial value of the previous time step, reducing the cumulative error, and thus accelerating the convergence of subsequent inner iterations; if local numerical oscillation problems occur, weighted averaging can be used, that is, the reconstructed field and the physical field of the previous step are weightedly mixed in a certain proportion (such as 80%-20%) to improve the smoothness and stability of the initial field.

[0059] Step 5: Inner iteration residual judgment and skip-step acceleration

[0060] During the inner iteration process, the residual changes are monitored in real time. When the iterative residual drops below a preset threshold (determined based on the actual original CFD pre-calculated residual), the current flow field is considered to be close to convergence, and subsequent redundant inner iteration steps can be skipped, advancing directly to the next physical time step, further improving overall computational efficiency.

[0061] Step 6: Dynamically update the feature space

[0062] To adapt to the time-dependent changes in the unsteady flow field, a new flow snapshot is collected at each step, a new POD decomposition is performed, and the base mode library is dynamically updated. By updating the feature space, the cumulative error problem caused by traditional fixed mode methods is avoided.

[0063] The proposed method achieves an adaptive balance between computational acceleration and error control through periodic sampling and dynamic POD decomposition, combined with residual feedback control. Compared with traditional unsteady flow field calculation methods, it significantly reduces the number of iterations and computation time while maintaining the accuracy of predictions of key physical quantities (such as aerodynamic forces and pressure distribution), making it suitable for large-scale engineering flow problems and real-time simulation scenarios.

[0064] To verify the effectiveness of this method, the NACA0012 airfoil inviscid pitching motion example is used. The variation of the airfoil angle of attack with time is as follows:

[0065] α(t)=α0+α m sin(ωt) (4)

[0066] Among them, the initial angle of attack α0 = 0.016°, the amplitude of the angle of attack α m =2.51°, ω is the pitch angular velocity.

[0067] Model geometry and mesh diagram are shown in Figure 2 , using an unstructured grid with a total grid number of 3,127,431 cells. The grid is denser in the boundary layer area, and the minimum boundary layer spacing is 10 -4 The specific operations of this example correspond to the above-mentioned method, including disturbance extraction, POD decomposition, main mode screening, flow field reconstruction and initial value update, residual monitoring skipping strategy, and dynamic update of the modal library. To verify the effectiveness of this patented method, simulations were performed for different working conditions. The verification conditions are shown in Table 1.

[0068] Table 1 Verification conditions

[0069]

[0070] In different working case examples, the method of the present invention achieves a reduction in the number of internal iteration steps while maintaining the prediction accuracy of key physical quantities such as aerodynamic force and pressure distribution. Specifically, the acceleration effect and accuracy retention results of working case 1 are as follows: Figure 3 and Figure 4(a)-Figure 4(d) As shown in the figure, the lift coefficient reconstructed by the proposed method is in good agreement with the full-order model result. The surface pressure difference at the same time and position does not exceed 10%, and more than 40% of the internal iteration steps are reduced. The acceleration effect and accuracy retention results of working condition 2 are shown in the figure. Figure 5 and Figure 6(a)-Figure 6(d) As shown in the figure, the lift coefficient reconstructed by the proposed method is in good agreement with the full-order model result. The surface pressure difference at the same time and position does not exceed 6%, and about 20% of the inner iteration steps are reduced. The acceleration effect and accuracy retention results of working condition 3 are shown in the figure. Figure 7 and 8(a) As shown in Figure 8(d), the lift coefficient reconstructed by the proposed method agrees well with the full-order model results. The surface pressure difference at the same time and location is approximately 2%, while the maximum reduction in internal iterations is over 60%. This demonstrates the method's ability to significantly accelerate computations and its effectiveness and practicality for widespread application in unsteady flow simulations.

[0071] The embodiments described above may be further combined or replaced, and the embodiments are merely descriptions of preferred embodiments of the present invention and do not limit the concept and scope of the present invention. Various changes and improvements made to the technical solutions of the present invention by persons of ordinary skill in the art without departing from the design concept of the present invention are within the scope of protection of the present invention. The scope of protection of the present invention is defined by the appended claims and any equivalent technical solutions.

Claims

1. A method for accelerating the calculation of unsteady flow fields based on POD dimensionality reduction of space-time mixing matrix, characterized in that: The following steps are involved: (1) In the dual-time iteration process, the physical quantity data of the unsteady flow field at different time steps are periodically collected, and the data at different time steps are stacked in columns to form a unified spatiotemporal mixed data matrix; (2) performing proper orthogonal decomposition (POD) on the spatiotemporal mixed data matrix, sorting the eigenvectors according to the eigenvalues, prioritizing the modes that contribute the most to the total disturbance energy, determining the energy retention ratio, dynamically selecting the number of modes, and constructing a low-dimensional feature space; (3) dynamically updating the low-dimensional feature space, i.e., repeatedly executing steps (1) and (2) within a preset iteration period to adapt to the characteristics of the flow field evolving over time; (4) reconstructing the flow field data based on the low-dimensional feature space, and using the reconstructed flow field as the initial value of the next dual-time iteration; (5) During the inner iteration process, the acceleration effect of the initial value improvement is judged by calculating the residual. When the inner iteration residual is lower than the preset threshold, some or all inner iteration steps of the current time step are skipped and the calculation of the next physical time step is directly entered to reduce the amount of calculation and speed up the solution.

2. The method according to claim 1, characterized in that The intrinsic orthogonal decomposition in step (2) uses the energy maximization principle to perform dimensionality reduction processing on the snapshot data, and the process includes: 2a. Assume that the space-time mixing perturbation matrix The construction method is as follows: in several consecutive time steps, disturbance data including velocity, pressure, density and other physical quantities are collected, and the spatial field data at each moment are stacked in the form of column vectors to form a time-arranged column vector sequence, thereby forming a unified space-time mixed disturbance matrix, and performing POD decomposition to obtain: Where x represents the position quantity, and its dimension varies with the dimension of the problem. It is the x coordinate in one dimension, the (x, y) coordinate in two dimensions, and the (x, y, z) coordinate in three dimensions. t represents the time quantity. m is the number of spatial grid points. N is the number of snapshots collected. B = [φ1,φ2,...,φ r ] is the orthogonal space mode set, S=diag(σ1,σ2,...,σ r ) is a diagonal matrix of singular values ​​sorted by size, and V is the time coefficient matrix; 2b. Sort the modes according to the contribution of disturbance kinetic energy. The energy contribution is determined by the singular value. Each mode φ i The contribution to the total energy is: Where r is the rank of the snapshot matrix, σ i For each mode, the corresponding singular value, E i Corresponding energy contribution for each mode; 2c. Select the first N main modes from largest to smallest according to the cumulative energy contribution rate, and satisfy the following adaptive conditions: Wherein, η is the preset energy retention threshold, which is set to 0.95 or 0.98, and is used to ensure that the low-dimensional modal subspace retains the main energy of the disturbance as much as possible; 2d. The selected main mode constructs the low-dimensional feature space B N =[φ1,φ2,...,φ N ], used for subsequent flow field reconstruction and accelerated calculation.

3. The method according to claim 1, characterized in that The low-dimensional feature space in step (3) adopts a dynamic update mechanism, that is, repeatedly performing snapshot acquisition and intrinsic orthogonal decomposition operations within a preset physical time step interval or iteration cycle, including: 3a. During the flow field calculation process, the latest snapshot data sets are collected again at set intervals (e.g., every several physical time steps) to form an updated perturbation matrix. 3b. Perform singular value decomposition on the disturbance matrix to re-obtain eigenvalues ​​and modes; 3c. Based on the updated energy spectrum, the main mode is filtered according to the preset energy threshold to construct a new low-dimensional feature space; 3d. The updated low-dimensional feature space is used for flow field reconstruction and initial value generation in subsequent physical time steps, thereby adapting to the time-dependent evolution of the main modes in the unsteady flow.

4. The method according to claim 1, wherein The residual threshold in step (5) is set to a fixed value determined according to the residual convergence trend of the original CFD iteration.

5. The method according to claim 4, characterized in that The stated value is 10 -1 to 10 -3 .

6. The method according to claim 4, characterized in that The stated values ​​are adjusted according to the complexity of the problem.

Images