Unsteady flow field accelerated calculation method based on spatiotemporal hybrid matrix POD dimension reduction

Abstract

The application discloses a kind of unsteady flow field accelerated calculation methods based on spatiotemporal hybrid matrix POD dimension reduction, comprising the following steps: (1) in the process of double time iteration, periodically collect the physical quantity data of unsteady flow field at different time steps, and the data of different time steps are stacked according to column, and constitute unified spatiotemporal hybrid data matrix;(2) the spatiotemporal hybrid data matrix is carried out eigenvalue orthogonal decomposition, and the low-dimensional feature space is constructed;(3) the low-dimensional feature space is dynamically updated;(4) reconstruct the flow field data based on the low-dimensional feature space, as the initial value of next step double time iteration;(5) in the process of inner iteration, the acceleration effect of initial value improvement is judged by calculating residual, to reduce the amount of calculation and accelerate solution.The application can reduce the number of inner iteration steps and the cost of unsteady flow field calculation under the premise of ensuring calculation accuracy, can effectively reduce the dimension of data, optimize the use of computing resources, has strong practical application value.

Description

Technical Field

[0001] This invention relates to an accelerated calculation method for unsteady flow fields based on spatiotemporal hybrid matrix POD dimensionality reduction, which is applicable to the rapid calculation of unsteady flow fields of aircraft. Background Technology

[0002] In modern engineering and scientific research, the simulation and computation of unsteady flow fields face severe challenges due to high computational complexity and long computational time. Especially when dealing with large-scale flow field evolution and unsteady aerodynamic problems, traditional numerical methods are often extremely costly, making it difficult to meet the real-time and efficiency requirements of engineering practice. Although continuous advancements in computational fluid dynamics (CFD) technology have improved computational speed, practical applications such as aerodynamic optimization, flow control, and environmental simulation still require enormous computational resources and long computation cycles. Therefore, accelerating the computation of unsteady flow fields to meet the requirements of efficient flow field simulation in engineering practice has become an important research direction.

[0003] Against this backdrop, Proper Orthogonal Decomposition (POD), as a data-driven method for reducing the order of flow fields, offers new possibilities for rapid computation of flow fields. By extracting the main characteristic patterns of flow field data, it reduces the dimensionality of high-dimensional flow field data to a lower-dimensional space, thereby effectively reducing computational costs while preserving the main flow characteristics. The POD method can significantly improve computational efficiency, thus becoming a popular research direction for rapid flow field prediction and optimization.

[0004] However, the main applications of existing POD methods are still concentrated in steady flow fields, while their application in unsteady computation remains significantly limited. Traditional POD order reduction processes are typically based on pre-defined modes, extracting modes from the statistically averaged state of the flow field or within a fixed time window. This method is effective for steady flows, but for unsteady flows, since the flow field modes evolve over time, fixed modes cannot adapt to transient changes, leading to large reconstruction errors and affecting computational stability and convergence speed. Furthermore, the computational model after POD order reduction is usually fixed and cannot adjust mode 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 established a two-layer POD for unsteady variable parameter problems to construct spatial and temporal basis functions with special properties, enabling the order-reduced model to satisfy boundary and initial conditions—Cao et al. combined POD with optimization algorithms, optimizing the POD basis coefficients through gradient optimization to reconstruct flow fields with smaller residuals, thus achieving a speedup of 2-3 times. Furthermore, Rapun et al. proposed the local POD method to accelerate the solution of parabolic unsteady problems. This method improves the robustness and efficiency of the solution process by alternately using numerical schemes and reduced-order models to solve the equations, and by determining the reduction solution time based on a prior estimate of the approximation error of the reduced-order model. In summary, the application of existing POD methods in unsteady flow calculations still has gaps or limitations. Summary of the Invention

[0005] The purpose of this invention is to provide an accelerated computational method for unsteady flow fields based on spatiotemporal hybrid matrix (POD) dimensionality reduction, thereby saving numerical computation costs. This invention constructs a reduced-order model using data-driven historical flow field snapshots and searches for a flow field with smaller residuals within the reduced-order space. Compared to existing POD methods, which typically reduce order based on preset modes and cannot adequately adapt to the dynamic changes in the flow field over time, leading to larger reconstruction errors in unsteady flows and affecting computational stability and convergence speed. This invention uses historical flow field snapshots to guide the dimensionality reduction process, breaking through the limitations of traditional POD static modes. It dynamically updates and optimizes the reduced-order model, ensuring that it reflects the latest state of the flow field at every moment, and adjusts the energy distribution and mode weights in real time based on historical flow field data. Furthermore, by maximizing energy dimensionality reduction within the reduced-order space, this invention prioritizes retaining the modes that contribute the most to the total disturbance energy, ensuring that important dynamic characteristics of the flow field are effectively preserved while significantly improving computational efficiency. More importantly, the reconstructed flow field is used as the initial state for the next moment's calculation, allowing each iteration to start from a more precise starting point, thereby effectively accelerating the convergence process of the numerical solution and providing a novel accelerated computational method for real-time simulation and prediction of unsteady flows.

[0006] The technical solution proposed in this invention is as follows:

[0007] An accelerated computation method for unsteady flow fields based on spatiotemporal hybrid matrix POD dimensionality reduction includes the following steps:

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

[0009] Step (2) Perform intrinsic orthogonal decomposition (POD) on the data matrix, sort the eigenvectors according to the eigenvalues, prioritize the retention of the modes that contribute the most to the total perturbation energy, determine the energy retention ratio, dynamically select the number of modes, and construct a low-dimensional feature space;

[0010] Step (3) Dynamically update the low-dimensional feature space, that is, repeat steps (1) and (2) within a preset iteration period to adapt to the characteristics of the flow field evolving over time;

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

[0012] In step (5), during the internal iteration process, the acceleration effect of the initial value improvement is judged by calculating the residual. When the internal iteration residual is lower than the preset threshold, some or all of the internal 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 accelerate the solution.

[0013] The intrinsic orthogonal decomposition in step (2) uses the energy maximization principle to reduce the dimensionality of the snapshot data. The process includes:

[0014] a. Let the spatiotemporal mixing perturbation matrix be... Its construction method is as follows: Perturbation data containing physical quantities such as velocity, pressure, and density are collected over several consecutive time steps. The spatial field data at each moment are stacked as column vectors to form a time-ordered column vector sequence, thus constituting a unified spatiotemporal hybrid perturbation matrix. POD decomposition is then performed on this matrix to obtain:

[0015]

[0016] Where x represents the location quantity, and the dimension varies with the dimension of the problem; for example, it is the x-coordinate in one dimension, (x, y) coordinates in two dimensions, and (x, y, z) coordinates 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 Let S be the orthogonal spatial mode set, where S = diag(σ1, σ2, ..., σ). rV 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 magnitude of their perturbation kinetic energy contribution. The energy contribution is determined by the singular values, and φ is the energy contribution of 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, there is a singular value, E i The energy contribution corresponding to each mode.

[0020] c. Select the top N principal modes in descending order of cumulative energy contribution rate, satisfying the following adaptive condition:

[0021]

[0022] Where η is a preset energy retention threshold, usually set to 0.95 or 0.98, used to ensure that the perturbation principal energy is retained as much as possible in the low-dimensional modal subspace;

[0023] d. The selected principal modes construct a 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, within a preset physical time step interval or iteration period, the snapshot acquisition and intrinsic orthogonal decomposition operations are repeatedly executed, including:

[0025] a. During the flow field calculation, the latest sets of snapshot data are collected again according to the set period (such as every several physical time steps) to form an updated perturbation matrix;

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

[0027] c. Based on the updated energy spectrum, the main modes are selected according to a 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 changes in the main mode of unsteady flow over time.

[0029] The residual threshold in step (5) is set to a fixed value determined based on the residual convergence trend of the original CFD iterations, preferably 10. -1 Up 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 intrinsic orthogonal decomposition (POD) to a compressible unsteady computational fluid dynamics (CFD) solver. It performs low-dimensional modeling on flow field snapshot data, predicts the approximate flow field for the next physical time step by selecting the main POD fundamental modes, and uses these as initial values ​​in the two-time-step iteration process. Compared to traditional methods that use solutions from the previous time step, this method effectively shortens the convergence time of the inner iterations, reducing the number of inner 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. Attached Figure Description

[0032] Figure 1 This is a schematic diagram of the method flow of the present invention.

[0033] Figure 2 This involves selecting the model geometry and mesh diagram.

[0034] Figure 3 This is a comparison of the flow field residual convergence process calculated in the airfoil subsonic example of this invention with the original iterative method.

[0035] Figures 4(a)-4(d) This is a comparison of the aerodynamic forces and surface pressures calculated in the airfoil subsonic example of this invention with the original iterative method;

[0036] Figure 4(a) shows the aerodynamic coefficient Cx as a function of angle of attack; Figure 4(b) shows the aerodynamic coefficient Cy as a function of angle of attack; Figure 4(c) shows the surface pressure distribution in the x-direction; and Figure 4(d) shows the surface pressure distribution in the y-direction. Figure 5 This is a comparison of the flow field residual convergence process calculated by the airfoil hypersonic example of this invention with the original iterative method.

[0037] Figures 6(a)-6(d) This is a comparison of the aerodynamic forces and surface pressures calculated by the airfoil hypersonic calculation example of this invention with the original iterative method;

[0038] Figure 6(a) shows the aerodynamic coefficient Cx as a function of angle of attack; Figure 6(b) shows the aerodynamic coefficient Cy as a function of angle of attack; Figure 6(c) shows the surface pressure distribution in the x-direction; and Figure 6(d) shows the surface pressure distribution in the y-direction. Figure 7 This is a comparison of the flow field residual convergence process calculated in the airfoil low pitch velocity example of this invention with the original iterative method.

[0039] Figures 8(a)-8(d) This is a comparison of the aerodynamic forces and surface pressures calculated in the airfoil low pitching speed example of this invention with the original iterative method;

[0040] Figure 8(a) shows the aerodynamic coefficient Cx as a function of angle of attack; Figure 8(b) shows the aerodynamic coefficient Cy as a function of angle of attack; Figure 8(c) shows the surface pressure distribution in the x-direction; and Figure 8(d) shows the surface pressure distribution in the y-direction. Detailed Implementation

[0041] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0042] This invention aims to effectively improve the computational efficiency of unsteady flow fields through intrinsic orthogonal decomposition technology. By maximizing energy for dimensionality reduction, the originally 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, reduce computational resource consumption, and is suitable for large-scale flow field simulation and real-time prediction. The method flow is as follows: Figure 1 As shown.

[0043] First, let's explain the principle of the POD method:

[0044]

[0045] Where U(x,t) is a physical quantity, t is a time representation, and x is a spatial representation; U is the average value; U'(x,t) is the fluctuation value; a i (t) is a time-dependent coefficient. These are spatially orthogonal basis functions or POD spatial basis modes. The POD method involves selecting the top N highest-energy modes and discarding some unimportant information.

[0046]

[0047] in This is an approximate deviation matrix.

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

[0049] Step 1: Flow field snapshot acquisition

[0050] During the dual-time-step iterative solution process, the state data of the unsteady flow field at different time points are periodically collected at several physical time steps. For each sampling time t... i Extract disturbing physical quantities (such as velocity, pressure, density, etc.) from the flow field to form a spatial distribution vector U'(x,t) i The spatial vectors collected at each time step are stacked in chronological order as column vectors to construct a unified spatiotemporal hybrid perturbation matrix. x represents only the position quantity, and the dimension varies with the dimension of the problem. For example, it is the x coordinate in one dimension, (x,y) coordinate in two dimensions, and (x,y,z) coordinate in three dimensions.

[0051]

[0052] in

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

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

[0055] Step 3: Principal Mode Extraction

[0056] In unsteady flows, the main flow structures (such as shear layers, unstable modes, and main vortices) typically carry the majority of the energy, while small disturbances 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 of energy, low-energy modes are truncated, and the number of modes is dynamically selected to form a dimensionality-reduced representation as shown in formula (2).

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

[0058] The perturbation flow field is reconstructed using the selected dominant mode and time coefficient, and the approximate flow field solution is obtained by adding the time average. The reconstructed flow field is used as the initial field for the new physical time step in the dual-time iteration process, replacing the traditional initial value of the previous time step, reducing accumulated errors and thus accelerating the convergence of subsequent internal iterations. If local numerical oscillations are encountered, a weighted average can be used, that is, the reconstructed field and the previous physical field are weighted and mixed in a certain proportion (such as 80%-20%) to improve the smoothness and stability of the initial field.

[0059] Step 5: Internal Iteration Residual Judgment and Skip-Step Acceleration

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

[0061] Step 6: Dynamically update the feature space

[0062] To adapt to the changes in the unsteady flow field over time, a new flow field 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 method of this invention 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 can significantly reduce the number of iterations and computation time while maintaining the prediction accuracy 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, an inviscid pitch motion example of the NACA0012 airfoil is used. The variation of the airfoil angle of attack with time is as follows:

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

[0066] Wherein, the initial angle of attack α0 = 0.016°, and the magnitude of the angle of attack α m =2.51°, where ω is the pitch angular velocity.

[0067] See the model geometry and mesh diagram. Figure 2 An unstructured mesh is used, with a total mesh size of 3,127,431 elements. The mesh is refined in the boundary layer region, with a minimum boundary layer spacing of 10. -4 The specific operations of this example correspond to the method described above, including: disturbance extraction, POD decomposition, principal mode selection, flow field reconstruction and initial value update, residual monitoring skipping strategy, and dynamic update of the mode library. To verify the effectiveness of this patented method, simulations were performed under different working conditions, as shown in Table 1.

[0068] Table 1 Verification Conditions

[0069]

[0070] In various operational scenarios, the method of this invention reduced the number of internal iteration steps while maintaining the prediction accuracy of key physical quantities such as aerodynamic forces and pressure distribution. Specifically, the acceleration effect and accuracy preservation results for operational scenario 1 are as follows: Figure 3 and Figures 4(a)-4(d) As shown, the reconstructed lift coefficient obtained by the proposed method agrees well with the results of the full-order model, while the surface pressure difference at the same time and location does not exceed 10%, and the internal iteration steps are reduced by more than 40%; the acceleration effect and accuracy retention results for condition 2 are as follows. Figure 5 and Figures 6(a)-6(d) As shown, the reconstructed lift coefficient obtained by the proposed method agrees well with the full-order model results, while the surface pressure difference at the same time and location does not exceed 6%, and an internal iteration step reduction of approximately 20% is achieved; the acceleration effect and accuracy retention results for condition 3 are as follows. Figure 7 and 8(a) As shown in Figure 8(d), the lift coefficient reconstructed by the proposed method agrees well with the results of the full-order model, while the surface pressure difference at the same time and location is at most about 2%, achieving a maximum reduction of over 60% in the number of internal iteration steps. This demonstrates the ability of the proposed method to significantly accelerate computation and its effectiveness and practicality in unsteady flow simulation.

[0071] The embodiments described above can be further combined or replaced, and these embodiments are merely descriptions of preferred embodiments of the present invention, not limitations on the concept and scope of the present invention. Various changes and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the design concept of the present invention are all within the protection scope of the present invention. The protection scope of the present invention is given by the appended claims and any equivalent technical solutions.

Claims

1. A method for accelerating the calculation of unsteady flow fields based on spatiotemporal hybrid matrix POD dimensionality reduction, characterized in that, Includes the following steps: (1) During the dual-time iteration process, 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 hybrid data matrix. (2) Perform intrinsic orthogonal decomposition on the spatiotemporal hybrid data matrix, sort the eigenvectors according to the size of the eigenvalues, prioritize the retention of the modes that contribute the most to the total disturbance energy, determine the energy retention ratio, dynamically select the number of modes, and construct a low-dimensional feature space; (3) Dynamically update the low-dimensional feature space, that is, repeat steps (1) and (2) within a preset iteration period to adapt to the characteristics of the flow field evolving over time. (4) Reconstruct the flow field data based on the low-dimensional feature space, and use the reconstructed flow field as the initial value for the next dual-time iteration; (5) During the internal iteration process, the acceleration effect of the initial value improvement is judged by calculating the residual. When the internal iteration residual is lower than the preset threshold, some or all of the internal 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 accelerate the solution. The intrinsic orthogonal decomposition in step (2) uses the energy maximization principle to reduce the dimensionality of the snapshot data. The process includes: 2a. Let the spatiotemporal mixing perturbation matrix be... The construction method is as follows: Perturbation data containing velocity, pressure, and density physical quantities are collected over several consecutive time steps. The spatial field data at each moment are stacked as column vectors to form a time-ordered column vector sequence, thus constituting a unified spatiotemporal hybrid perturbation matrix. POD decomposition is then performed to obtain: Where x represents the location quantity, and the dimension varies with the dimension of the problem: x-coordinate in one dimension, (x, y) coordinate in two dimensions, and (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 Let S be the orthogonal spatial mode set, where S = diag(σ1, σ2, ..., σ). r ) is a singular value diagonal matrix sorted by size, and V is the time coefficient matrix; 2b. Sort the modes according to the magnitude of their perturbation kinetic energy contribution. The energy contribution is determined by the singular values. φ for each mode... i The contribution to the total energy is: Where r is the rank of the snapshot matrix, σ i For each mode, there is a singular value, E i Energy contribution corresponding to each mode; 2c. Select the top N principal modes in descending order of cumulative energy contribution rate, satisfying the following adaptive condition: Wherein, η is a preset energy retention threshold, which takes the value of 0.95 or 0.98, and is used to ensure that the perturbation principal energy is retained in the low-dimensional modal subspace; 2d. Construct a low-dimensional feature space B using the selected principal modes. N =[φ1,φ2,...,φ N ], used for subsequent flow field reconstruction and accelerated calculation.

2. The method according to claim 1, characterized in that, The low-dimensional feature space in step (3) adopts a dynamic update mechanism, that is, within a preset physical time step interval or iteration period, the snapshot acquisition and intrinsic orthogonal decomposition operations are repeatedly executed, including: 3a. During the flow field calculation, the latest sets of snapshot data are collected again according to the set period to form an updated perturbation matrix; 3b. Perform singular value decomposition on the perturbation matrix to re-obtain eigenvalues ​​and modes; 3c. Based on the updated energy spectrum, the main modes are selected 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 changes in the main mode of unsteady flow over time.

3. The method according to claim 1, characterized in that, The residual threshold in step (5) is set as a fixed value determined based on the residual convergence trend of the original CFD iteration.

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

5. The method according to claim 3, characterized in that, The specified value is adjusted according to the complexity of the problem.

Images