Method for predicting PM2.5 concentration field of road surface area

Abstract

The invention relates to a method for predicting a PM2.5 concentration field of a pavement area. The method comprises the following steps: using simulation software FLUENT to obtain two-dimensional flow field data when a vehicle moves under a plurality of different working conditions, and extracting target physical quantities of a plurality of sampling points on a plane under the different working conditions through numerical simulation post-processing software so as to construct a multi-working-condition vehicle movement flow field snapshot matrix; carrying out order reduction operation on the snapshot matrix by adopting an intrinsic orthogonal decomposition method to obtain main modes of the flow field and mode coefficients of the main modes; constructing a mapping relation between a modal coefficient and a working condition through a least square support vector machine model to quickly obtain a modal coefficient corresponding to a target working condition, and linearly combining the modal coefficient with the main modal to obtain flow field data corresponding to the target working condition; and verifying the accuracy of the reduced-order model by calculating an error between the reduced-order model and the full-order model. The method can be suitable for a dynamic mesh model, and provides an order reduction method for dynamically changing flow field calculation.

Description

Technical Field

[0001] The present invention belongs to the field of numerical simulation and relates to a method for predicting the PM2.5 concentration field in a road surface area. Background Art

[0002] Vehicle movement has a significant impact on the airflow and pollutant diffusion on the road. Studying the pollutant diffusion pattern on the road under vehicle movement is crucial for avoiding inhaling pollutants and protecting the health and safety of the population. Although the traditional dynamic mesh model calculation method can provide an accurate flow field structure and pollutant distribution near the vehicle body during vehicle movement, it often consumes a large amount of computing resources and time costs. This problem is more prominent especially in the case of multiple working conditions or large differences in grid sizes within the computational domain.

[0003] Proper orthogonal decomposition has been widely applied in the field of numerical simulation. It simplifies the snapshot matrix into a linear formula by extracting the main modes of the flow field, achieving a high approximation of the original flow field. However, proper orthogonal decomposition is often only applied to static models. This is because the dynamic mesh model can cause situations such as changes in the positions of grid nodes and different numbers of grids under different working conditions, and the snapshot matrix of the traditional method cannot be constructed. Therefore, it is often difficult to simplify the flow field calculation involving dynamic models through a reduced-order model. How to design a comprehensive and efficient analysis method has become a difficult problem to be solved urgently. Summary of the Invention

[0004] In view of the deficiencies of the prior art, a method for predicting the PM2.5 concentration field in a road surface area provided by this application includes the following steps:

[0005] Step S1: Use the numerical simulation software FLUENT to calculate the PM2.5 concentration field data near the vehicle body on the simulated road surface area during a certain period of time under several different working conditions. In the post-processing software CFD-POST, use multiple groups of sampling lines to extract the PM2.5 concentration field data for each working condition. The sampling lines are parallel to each other and of equal length, the interval between adjacent two sampling lines is less than 0.05 m, several sampling points are provided on each sampling line, the interval between sampling points is less than 0.05 m, and the plane composed of all sampling lines completely covers the road surface area;

[0006] Step S2: For each working condition, export the data sequences on all sampling lines to a CSV format file by the post-processing software CFD-POST, and through data conversion, make the data sequence groups on all sampling lines for each working condition form a data column vector:

[0007] The data conversion includes: marking serial numbers for the data in the data sequence of each sampling line: "o-p", indicating the p-th sampling point in the sampling line with the initial number "o", so that the data sequences of all sampling lines under each working condition are arranged in ascending order of serial numbers from top to bottom in the data column vector;

[0008] Step S3: By repeating Step S2, data column vectors composed of data sequences of all sampling lines under several different working conditions are obtained. All the data column vectors are horizontally combined to obtain a matrix containing PM2.5 concentration field data near the vehicle body in the simulated road surface area under different working conditions, named the snapshot matrix;

[0009] Step S4: Import the snapshot matrix into MATLAB in the form of a numerical matrix, and perform POD reduction calculation on the snapshot matrix, including:

[0010] Construct the snapshot matrix S = [x1, x2,..., x m , where S is an n×m matrix, each column is a column vector, n is the dimension of the column vector, and m is the number of column vectors;

[0011] Use the built-in function SVD in MATLAB to perform singular value decomposition on the snapshot matrix S in Step 51, and obtain three matrices U, Σ, V T ;

[0012] S = UΣV T

[0013] U is an n×n orthogonal matrix, and the column vectors u1, u2,..., u n are the modes;

[0014] Σ is an n×m diagonal matrix, and the diagonal elements σ1, σ2,..., σ r are the singular values;

[0015] V T is an m×m orthogonal matrix, containing the weights of the modes in terms of time or operating conditions;

[0016] Extract the first r modes with larger singular values from the matrix U as the dominant modes of the system:

[0017] U r = [u1, u2,..., u r

[0018] The sum of the energies of the r modes needs to be greater than 90% of the total modal energy, and the square σ d of the singular value σ d 2 represents the energy contribution of the d-th mode;

[0019] ​Calculate the modal coefficient matrix A according to the product of Σ and V, A = Σ·V T ;

[0020] The modal coefficient a of the d-th mode in the modal coefficient matrix A d = σ d ·v d , v d is the corresponding column vector in matrix V T ;

[0021] Step S5: Use the least squares support vector machine model LS-SVM to construct the mapping relationship between the modal coefficients and the operating condition parameters, and then obtain the modal coefficients of the un-simulated operating conditions;

[0022] Step S6: According to the modal coefficients A r = [a1, a2,..., a r of the un-simulated operating conditions obtained in Step S5, substitute the modal coefficients of the un-simulated operating conditions into the simplified linear combination reconstruction formula to obtain the concentration field data matrix of the un-simulated operating conditions:

[0023] c = c0 + u1a1 + u2a2 + … + u r a r ;

[0024] c0 is the average concentration value of all concentration fields in Step S1.

[0025] In the said Step S1, the method for simulating different operating conditions is:

[0026] Step 11. When using the CFD method for numerical simulation to obtain the flow field data, in order to simulate the wind profile at the velocity inlet, the magnitude of the environmental wind is set to:

[0027]

[0028] where, U z represents the velocity magnitude at a height Z from the road surface, and U ref represents the wind speed magnitude at the reference height z ref , and its value depends on the average wind speed magnitude observed by the local meteorological station at a height of z ref = 10m on the local open and flat ground within 10 minutes; α is the ground roughness index, taking 0.22;

[0029] Step 12. Activate the DO radiation model and the solar ray tracing method to simulate the solar radiation at different times of the day. In the solar radiation calculator, set the time zone, longitude and latitude of the experimental location, as well as the sunshine time, and determine the direct and diffuse radiation intensities at the selected time;

[0030] Step 13. Simulate the process of particulate matter diffusion through the discrete phase model in the numerical simulation software FLUENT. At the same time, set the incident parameters of PM2.5 particulate matter in this model, including incident flow rate, injection velocity, and incident duration;

[0031] Among them, the incident flow rate is obtained according to the motor vehicle emission factor EF. The PM2.5 mass flow rate Q of a vehicle that takes x seconds to travel 1 km is expressed as:

[0032]

[0033] The tail gas emission velocity at different driving speeds refers to the experimental results of on-vehicle tail gas collection;

[0034] The temperature of the incident particulate matter is 380K, and the duration needs to be greater than the dynamic mesh movement time;

[0035] Step 14. To simulate the vehicle movement, activate the dynamic mesh (Dynamic meshing) model in the CFD model.

[0036] In step S5, use the least squares support vector machine model (LS-SVM) to construct the mapping relationship between the modal coefficients and the operating conditions parameters, and then obtain the modal coefficients of the un-simulated operating conditions, including:

[0037] Optimize the weight vector of the least squares support vector machine algorithm, and construct a multi-output regression model using the weight vector, Gaussian radial basis kernel function, and error vector;

[0038] Construct a training set with the r modal coefficients obtained in step S4 and input it into the multi-output regression model for processing, and output the predicted values of the r modal coefficients.

[0039] The beneficial effects produced by the present invention are: reducing the calculation cost, accelerating the calculation speed, providing an efficient and accurate analysis method, which can significantly improve the calculation efficiency in the case where the dynamic mesh model will cause the transformation of the grid node positions and the different numbers of grids under different operating conditions, making the POD reduced-order model applicable to the CFD calculation containing the dynamic mesh model, and solving the problem of high requirements for calculation time and calculation resources. Description of the Drawings

[0040] Figure 1 It is the setting interface diagram of the velocity inlet in the CFD model;

[0041] Figure 2 It is the setting interface diagram of the radiation model in the CFD software;

[0042] Figure 3 It is the setting interface diagram of the injection source attributes in the CFD software;

[0043] Figure 4Graph showing the correspondence between the driving speed and the exhaust emission speed of a vehicle;

[0044] Figure 5 Interface diagram for setting the properties of the dynamic mesh region in CFD software;

[0045] Figure 6 Interface diagram for setting the parameters of mesh re - meshing in CFD software;

[0046] Figure 7 Schematic diagram of using a sampling line in CFD - POST software;

[0047] Figure 8 Schematic diagram for exporting the data of 500 sampling lines using POST export work;

[0048] Figure 9 For Figure 8 Schematic diagram of the exported data;

[0049] Figure 10 For Figure 9 Schematic diagram after data conversion in;

[0050] Figure 11 Schematic diagram after the snapshot matrix is imported into MATLAB in the form of a numerical matrix;

[0051] Figure 12 Schematic diagram of the error between the prediction result and the true value;

[0052] Figure 13 Flow schematic diagram of the method for reduced - order calculation of the PM2.5 concentration field near a moving vehicle of the present invention. Detailed implementation manners

[0053] In order to make the objectives, technical solutions, and advantages of the embodiments of the present disclosure clearer, the technical solutions of the embodiments of the present disclosure will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present disclosure. Obviously, the described embodiments are part of the embodiments of the present disclosure, rather than all of the embodiments. Based on the described embodiments of the present disclosure, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the scope of protection of the present disclosure.

[0054] The methods or software used in this embodiment:

[0055] Method 1: Computational Fluid Dynamics (CFD) is a technology that simulates and analyzes physical phenomena such as fluid flow, heat transfer, and chemical reactions through numerical methods and algorithms. Its core is to use a computer to solve the mathematical equations describing fluid motion (such as the Navier-Stokes equations). Through steps such as geometric modeling, mesh generation, physical modeling, numerical solution, and post-processing, detailed information about the flow field is obtained. CFD is widely used in fields such as aerospace, automotive manufacturing, energy, architecture, and environmental engineering, providing an efficient and economical means for optimizing designs and solving problems.

[0056] Method 2: POD (Proper Orthogonal Decomposition) is a method for data dimensionality reduction and pattern extraction, widely used in fields such as fluid mechanics, structural mechanics, and signal processing. By decomposing complex high-dimensional data into a set of orthogonal basis functions and corresponding time coefficients, POD can extract the dominant dynamic characteristics of the system, capture the main energy distribution, and key patterns. Its core is to project the data into a low-dimensional space through eigenvalue decomposition or singular value decomposition (SVD), thereby reducing the computational complexity and improving the interpretability of physical problems. In this embodiment, first, numerical calculations are performed through Fluent to obtain PM2.5 concentration field data near a moving vehicle under different working conditions, thereby forming a snapshot matrix; then, by calculating the covariance matrix and solving for eigenvalues and eigenvectors, a set of POD principal modes are extracted; finally, the original data is projected onto these POD principal modes to obtain modal coefficients.

[0057] Method 3: Support Vector Machine (SVM) is a classical supervised learning method that finds the optimal hyperplane in the feature space for classification or regression, with good non-linear fitting ability and generalization performance. In this embodiment, this method is used to construct the mapping relationship between working condition parameters and modal coefficients.

[0058] Software 1: ANSYS Fluent is a powerful commercial computational fluid dynamics (CFD) software used to simulate complex problems such as fluid flow, heat transfer, turbulence, and chemical reactions. It helps users efficiently solve fluid dynamics problems in industry and academia by supporting multiple physical models (such as turbulence, two-phase flow, and combustion), powerful mesh generation tools (structured and unstructured meshes), and the extended ability of user-defined functions (UDF). In this embodiment, ANSYS Fluent 2021R1 is selected as the CFD numerical calculation solver, and the dynamic mesh model, discrete phase model, and DO radiation model are activated to simulate vehicle motion, particulate matter diffusion, and solar radiation respectively. By solving under the state of activating these models, the particulate matter diffusion situation near the moving vehicle under different working conditions including different wind, heat weather conditions, and vehicle speeds can be obtained, thereby obtaining the original CFD dataset that can be used for POD dimensionality reduction.

[0059] Software 2: CFD-POST is a professional tool for post-processing the results of computational fluid dynamics (CFD) simulations. It can visualize, analyze, and generate reports on simulation data, helping users deeply understand the flow field characteristics and design performance. In this embodiment, CFD-POST 2021R1 is selected for post-processing the calculation results to extract the concentration field data in the street canyon under different working conditions.

[0060] Software 3: MATLAB (Matrix Laboratory) is a data processing and calculation tool centered around matrices. Its programming language and operating environment are optimized for matrix operations, capable of efficiently handling matrix creation, operations, decomposition, solution, and visualization, supporting a wide range of mathematical functions such as linear algebra, Fourier transform, and matrix decomposition. It is an ideal platform for matrix data analysis, simulation modeling, and algorithm development, especially suitable for scientific research and engineering calculation scenarios that require frequent operation of large-scale matrices. In this embodiment, MATLAB 2021b is selected for post-processing the calculation results to establish the mapping relationship between working condition parameters and modal coefficients.

[0061] Embodiment 1

[0062] A method for predicting the PM2.5 concentration field in the road surface area, and the implementation steps are as follows:

[0063] Step 1. First, use the FLUENT calculation software to calculate the PM2.5 concentration data near the vehicle body within a certain period under 75 different working conditions in the road surface area, that is, perform 75 groups of CFD calculations with different initial conditions. The boundary conditions include:

[0064] 1) Ambient gradient wind

[0065] Due to the influence of ground roughness, environmental wind usually occurs in the form of gradient wind. In the velocity inlet setting interface of the CFD model, a Profile file is used to define the incoming flow velocity at the velocity inlet boundary. As Figure 1 shown, in the velocity inlet setting interface, the magnitude of the velocity is defined along the normal direction of the boundary. The velocity is defined relative to a fixed reference point (absolute reference system). The magnitude of the velocity is determined by the following formula (1). The supersonic / initial gauge pressure is set to 0, the turbulence intensity is set to 5%, and the ratio of turbulent viscosity to laminar viscosity is set to 10:

[0066]

[0067] where U z represents the magnitude of the velocity at a height Z from the road surface. U ref represents the magnitude of the velocity at the reference height z ref from the road surface. In wind energy engineering, 10 meters is usually taken as the general reference height z ref for meteorological observations. Different U ref at this reference height are used to characterize different levels of environmental wind speed in the street canyon. U ref = 2.0, 3.3, 5.4, 7.5, 10.8 m / s, respectively representing different wind force levels (level 2 - 6) in the street canyon. α is the ground roughness index. According to the research results of Wang [1] et al., α = 0.22.

[0068] 2) Solar radiation

[0069] The DO radiation model (Discrete Ordinates Model Theory) and the solar ray tracing method (Raytracing approach) are used to simulate the solar radiation at different times of a day. As Figure 2 shown, in the setting interface of the radiation model in the CFD software, the discrete coordinate model (DO) is used to simulate the radiation transfer, and the solar ray tracing method is used to simulate the solar radiation. In the solar radiation calculator, the time zone, longitude and latitude of the experimental site, and the sunshine time are set. The solar direction vector is calculated using the direction of the solar energy calculator. The time step / solar radiation energy update is set to 10. This setting means that the direction and energy intensity of the solar radiation are recalculated and updated every 10 time steps, which can significantly reduce the calculation amount while ensuring the simulation accuracy. According to the research of Mao Qianjun [2] , the direct and diffuse radiation intensities at the selected time are determined.

[0070] The direct radiation intensity represents the radiation intensity directly from the sun.

[0071] The scattered radiation intensity represents the solar radiation intensity after atmospheric scattering.

[0072] Spectral fraction [V / (V + IR)]: Set to 0.5, which may represent the ratio of visible light to infrared light.

[0073] Among them, after the radiation model is set to DO, it is applicable to radiation problems in all optical depth ranges. At the same time, the radiative heat transfer between gas and particles can be considered, and the solar ray tracing method can efficiently add solar radiation as a heat source to the computational domain.

[0074] 3) Setting of PM2.5 incident parameters

[0075] In the interface for setting the injection source properties in the CFD software as shown in Figure 3 , set the exhaust port in the vehicle model as the particulate injection surface and set its injection parameters for simulating the injection of particulate matter (such as PM2.5). The particle type is set to inert, indicating that the particles do not change their chemical composition or participate in chemical reactions during the flow process; the particle material is set to anthracite; the particle diameter distribution is set to uniform, and the discrete phase domain is set to none; the particle diameter is set to 2.5e -6 m, the particle temperature is set to 380 K, and the start time and stop time are set to determine the injection duration of PM2.5 so that the injection duration is sufficient to cover the dynamic mesh movement time. The dynamic mesh movement time refers to the time required for the mesh to be dynamically adjusted to adapt to the flow changes during the simulation; Figure 3 The total flow rate in can be calculated according to the emission factor (EF) of light-duty vehicles (LDV) when driving under normal weather conditions. Existing studies have shown that the average PM2.5 emission factor of light-duty vehicles (LDV) under normal weather conditions is (7.402) mg / (veh·km). Then, when the time required for the vehicle to travel 1 km is x seconds, the PM2.5 mass flow rate Q of the vehicle is expressed as:

[0076]

[0077] Figure 3 The magnitude of the velocity in represents the set velocity of the particulate matter at the initial injection moment, and the real-time velocity during subsequent movement will change dynamically under the influence of external forces such as fluid drag force and gravity. This value is obtained according to the on-vehicle exhaust gas collection experiment [3] as shown in Figure 4As shown, when the vehicle is traveling at speeds of 10, 20, 30, 40, and 50 km / h, the exhaust gas emission speeds are 1.99, 2.29, 2.59, 3.36, and 4.45 m / s respectively.

[0078] 4) Dynamic mesh parameters

[0079] In the interface for setting the properties of the dynamic mesh region in the CFD software as shown Figure 5 below, the vehicle motion is realized using the dynamic mesh model. First, set the rigid body motion properties for the vehicle shell at different positions within the dynamic mesh model, as well as the properties of the relative motion relationship between the rigid body and other regions within the dynamic mesh, and use the user-defined function UDF to implement the process of the vehicle traveling at a constant speed in a certain direction:

[0080] The parameters set by the user-defined function UDF include: vehicle speed, speed direction, and simulation motion time;

[0081] Taking the example of traveling in the y direction at a speed of 5.55 m / s for 12 s, the user-defined function UDF is:

[0082] ((move transient 12 0)

[0083] (time 1 2 3 4 5 6 7 8 9 10 11 12)

[0084] (v_y 5.55 5.55 5.55 5.55 5.55 5.55 5.55 5.55 5.55 5.55 5.55 5.55))

[0085] Subsequently, as shown Figure 6 below, in the interface for setting the mesh rezoning parameters in the CFD software, select the unified mesh rezoning method (Unified Remeshing) to achieve mesh reconstruction. This method is optional in FLUENT 2021R1 and later versions. It should be noted that in order to ensure the mesh quality during the motion as much as possible, the minimum cell size needs to be set to 0.4 times the minimum size in the initial mesh, while the maximum cell skewness needs to be set to 1.4 times the maximum size in the initial mesh.

[0086] In the above step 1, concentration field data under 75 different working conditions composed of the intersection of 3 different radiation intensities, 5 environmental wind speeds, and 5 different vehicle speeds were obtained through numerical calculations.

[0087] Step 2. As shown Figure 7As shown, in the CFD-POST software, sampling lines are used to extract the data of the PM2.5 concentration field on the pavement area simulated in step 1. The sampling lines include 500 parallel and equal-length sampling lines, 1000 sampling points are evenly arranged on each sampling line, the spacing between adjacent sampling lines is equal, and all sampling lines cover the pavement area. That is, a total of 500,000 sampling points are set for the extraction of PM2.5 concentration field data on this plane. Taking the pavement area as an example, in order to capture the distribution characteristics of the concentration field as much as possible and avoid the decrease in calculation efficiency due to too many sampling points, the distance between sampling points should be less than 0.05m;

[0088] Step 3. Use CFD-POST to export PM2.5 concentration data at 1000 sampling points on 500 sampling lines under 75 different working conditions. At this time, the PM2.5 concentration data exists in Excel as a column vector and contains a header (such as "Line 1"). The data on the 500 sampling lines are arranged from top to bottom according to the number of the sampling line header from small to large. Then, use Excel to convert the data in the original column vector, including removing the header number, and deleting the non-PM2.5 concentration data in the original column vector without changing the order of the PM2.5 concentration data. The number of rows in the converted data set is equal to 500,000, and the column vector does not contain the sampling line number that comes with CFD-POST when exported.

[0089] The specific conversion process is as follows: Figure 8 As shown in the figure, use the POST export function ("EXPORT") to export the data of 500 sampling lines. First, select the sampling lines of the data to be exported. In the post-processing software POST, the 500 sampling lines are named Line 1 to Line 500. Next, uncheck "Export geometric information". This step will make the exported data contain only physical quantity information, not geometric information. Finally, select the type of data you want to obtain as "DPM Concentration", that is, PM2.5 concentration data.

[0090] like Figure 9 As shown in the figure, by exporting the data in csv format, a 523120-dimensional data column vector is obtained, which contains the sampling line header number: [Name] and [Data]. In addition, it should be noted that the column vector headers: [Name] and [Data] are deleted at this time, making it a continuous data column vector composed of PM2.5 concentration data.

[0091] The above conversion uses the built-in functions of EXCEL to complete the task, such as Figure 10 As shown, column B is the original data, column F is the converted results, and the functions used in columns A, E, and F are:

[0092] Function I in column A: =IFERROR(IF(LEFT(B3,4)="Line",RIGHT(B3,LEN(B3)-5)&"-0",IF(B6<>"",LEFT(A5,FIND("-",A5)-1)&"-"&RIGHT(A5,LEN(A5)-FIND("-

[0093] ",A5))+1,"")),"") When using Function I to process the original data in column B, when it is detected that the original data in column B contains the string of the sampling line header and the sampling line number, the string containing the sampling line header and the sampling line number is converted into a numerical string connected by the sampling line number and "-0", and is stored in the same row of column A as the new number of the first non-empty data row detected after the string in the original data of column B. The row number of the non-empty data row in the original data of column B is c. Then, the original data of the (c + 1)-th row is detected row by row. If the original data of the (c + 1)-th row is not empty, the number on the right side of "-" in the new number of the c-th row of the original data is incremented by 1 and stored in the (c + 1)-th row of column A as the new number of the original data of the (c + 1)-th row. If the original data of the (c + 1)-th row is empty, an empty string is stored in the (c + 1)-th row of column A as the new number of the original data of the (c + 1)-th row. Function I is repeatedly used to traverse and process the remaining original data in column B to re-number the original data in column B.

[0094] Function II in column E: =(INT((ROW()-1) / 1000)+1)&"-"&ROW()-INT((ROW()-1) / 1000)*1000. Use Function II to create a serial number in the format of "X-Y" row by row in column E, where "X" is an integer starting from 1, and "Y" is a consecutive integer starting from 1. When Y increments from 1 to 1001, X advances one digit until it reaches 500, and the serial numbers are arranged in the order of natural numbers.

[0095] Function III in column F: =VLOOKUP(E1,$A$6:$B$40600,2,FALSE). Use Function III to find the new number in column A that exactly matches the serial number in each row of column E. If found, return the original data in column B that is in the same row as the found new number in column F. If no match is found, or the searched new number is invalid in column A, VLOOKUP will return the error value #N / A.

[0096] The obtained column F is a continuous column vector composed of PM2.5 concentration data.

[0097] Step 4. Arrange the concentration field data under 75 working conditions in order from left to right to construct a concentration field snapshot matrix, where each column is a column vector containing 500,000 data, representing the concentration field information under a specific working condition.

[0098] Step 5. Import the snapshot matrix into MATLAB in the form of a numerical matrix, as Figure 11 shown, and perform a reduction operation on the snapshot matrix through the POD model to obtain the dominant modes and modal coefficients of this series of concentration fields obtained in Step 4. Its essence is to simplify the data in matrix form into a linear combination of POD modes and modal coefficients, that is, multiple vectors are weighted and summed to obtain a new vector, where each vector has a corresponding scalar weight.

[0099] Specifically, Step 5 includes:

[0100] Step 51: Construct the snapshot matrix S = [x1, x2,..., x m , where S is an n×m matrix, each column is a snapshot (column vector), n is the dimension of the system state (column vector), and m is the number of snapshots (column vectors).

[0101] Step 52: Use the built-in function SVD in MATLAB to perform singular value decomposition on the snapshot matrix S in Step 51 [4] to obtain three matrices U, Σ, V T .

[0102] S = UΣV T

[0103] U is an n×n orthogonal matrix, and the column vectors u1, u2,..., u n are the modes.

[0104] Σ is an n×m diagonal matrix, and the diagonal elements σ1, σ2,..., σ r are the singular values.

[0105] V T is an m×m orthogonal matrix, containing the weights of the modes in terms of time or operating conditions.

[0106] The square of the singular value σ d σ d 2 represents the energy contribution of the d-th mode. Based on the magnitudes of the singular values, the first r modes with larger singular values can be selected as the dominant modes of the system. Extract the first r modes from matrix U:

[0107] U r = [u1, u2,..., u r

[0108] ​Generally, the sum of the energies of the retained r modes needs to be greater than 90% of the total modal energy.

[0109] Step 53: Calculate the modal coefficients based on the matrix Σ and the matrix V in Step 52. The modal coefficients describe the weight of each mode over time or operating conditions.

[0110] Calculate the modal coefficient matrix A:

[0111] X = Σ·V T

[0112] Calculate the modal coefficient a of the d-th mode d :

[0113] a d = σ d ·v d

[0114] v d is the corresponding column vector in the matrix V T among them.

[0115] Step 6. Use the modal coefficient data of known operating conditions to construct a model that can predict the modal coefficients of the system under un-simulated operating conditions.

[0116] Through Step 5, different dominant dynamic modes in this series of flow fields, namely POD modes, and the modal coefficients showing the variation laws of these dominant dynamic modes with operating condition parameters are obtained. In order to complete the reconstruction of the un-simulated operating condition flow field with the help of POD modes and the modal coefficients of un-simulated operating conditions, a support vector machine model (SVM) is used to construct the mapping relationship between the modal coefficients and the operating condition parameters.

[0117] In this embodiment, first, 3 groups of operating condition parameters are used as input features, and 30 groups of modal coefficients are used as target variables. Then, the LS-SVM toolbox of MATLAB is used to preprocess the data, including normalization to eliminate the influence of different dimensions on the training process. In the model training stage, the radial basis kernel function (RBF) is selected as the kernel function, and the grid search and cross-validation methods are combined to optimize the kernel parameters and the penalty coefficient to ensure the optimal performance of the model. Through training, the non-linear mapping relationship between the operating condition parameters and the modal coefficients is obtained, and then the modal coefficients of un-simulated operating conditions are obtained, and the corresponding flow field data are obtained through the linear combination described in Step 5;

[0118] Specifically, referring to ZL201910562179.2, the least squares support vector machine model LS-SVM is used to construct the mapping relationship between the modal coefficients and the operating condition parameters, and then the modal coefficients of un-simulated operating conditions are obtained. Step 6 includes:

[0119] Step 61: Construct a multi-output regression model of the least squares support vector machine LS-SVM in the following form:

[0120]

[0121] where, f j (x) is the predicted value of the output modal coefficient at the j-th working condition; γ i (j) is the Lagrange multiplier, i.e., the weight vector, used for the optimization of the j-th output; K(x, x i ) is the kernel function, which calculates the similarity between the input working condition sample x and all training working condition samples x i , and takes this similarity as the weight, implicitly mapping the original data to a high-dimensional space, so as to provide a basis for the non-linear weighted combination of the predicted values f j (x) of each output modal coefficient, enabling the model to adaptively fuse the information of the training samples based on the similarity degree between samples; b (j) is the bias term; l is the number of training working condition samples; k is the output dimension (the number of modal coefficients).

[0122] Step 62: Optimize the weight vector. By solving the linear equations composed of the optimization objective function and the constraint conditions, γ and b can be obtained.

[0123] The optimization objective function and the constraint conditions can be expressed as:

[0124] The optimization objective function of LS-SVM is:

[0125]

[0126] where, γ = [γ1, γ2, …, γ l T is the weight vector; is the error vector, representing the error of the i-th sample in each output dimension; C is the regularization parameter, controlling the trade-off between the error penalty and the model complexity, ∥γ∥ 2 represents the square of the L2 norm of the weight vector. Using the square of the L2 norm (∥γ∥ 2 ) as the regularization term is mainly to simplify the optimization process (such as derivative calculation and closed-form solution derivation) and keep consistent with the structure of the least squares loss function, so as to effectively suppress overfitting while ensuring the solvability of the model. The optimization objective function aims to determine the model parameters γ and b by minimizing the prediction error and constraining the model complexity to prevent overfitting.

[0127] The square of the L2 norm, also known as the square of the Euclidean norm, is the square of the Euclidean length of the vector.

[0128] ​For a vector γ = [γ1, γ2, …, γ l T , its squared L2 norm is defined as:

[0129] ∥γ∥ 2 = γ1 2 + γ2 2 + … + γ l 2 ;

[0130] The constraint condition is:

[0131]

[0132] For the convenience of solving, introduce the Lagrange multiplier λ and construct an unconstrained optimization function:

[0133]

[0134] According to the KKT (Karush-Kuhn-Tucker) conditions, the optimized matrix form is:

[0135]

[0136] where G = (1, 1, …,) T , Ω ∈ R 1×1 is a 1×1 dimensional vector, and Ω is the kernel matrix, defined as:

[0137] Ω ij = K(x i , x j )

[0138] Step 63: Select the Gaussian radial basis kernel function (RBF kernel) to capture the nonlinear relationship between the input (operating condition parameters) and the output (modal coefficients).

[0139]

[0140] where x i , x j represent the input operating condition parameter combinations of the i-th and j-th operating conditions (such as different combinations of environmental wind speed, solar radiation intensity, vehicle speed, etc.) respectively. Its role is to implicitly map the low-dimensional nonlinear relationship to a high-dimensional linearly separable feature space by measuring the geometric distance between the two operating condition parameter combinations, thereby establishing a global nonlinear mapping model between the input operating condition parameters and the modal coefficients. ε is the kernel width parameter, which controls the radial action range of the function.

[0141] ​Based on the grid search method, all parameter combinations are generated within the candidate ranges of the predefined regularization parameter C and the kernel function width ε. For each set of parameters, its average validation error is calculated through ten-fold cross-validation. Finally, the parameter combination with the minimum validation error is selected as the optimal hyperparameter configuration. The parameter combination is the combination of the regularization parameter C and the kernel function width ε.

[0142] Step 7. Since r modes and coefficients that have a greater influence on the flow field are retained according to the extraction of the main modes in Step 5 (r << n), and the modal coefficients A of the un-simulated working conditions are obtained according to Step 6 r =[a1, a2,..., a r , combined with the average concentration value c0 of the original concentration field, the concentration field data matrix of the un-simulated working conditions can be constructed:

[0143] c = c0 + u1a1 + u2a2 + … + u r a r

[0144] Therefore, the dominant modes and the modal coefficients of the un-simulated working conditions are obtained through Step 5 to Step 6. By substituting the modal coefficients of the un-simulated working conditions into the simplified linear combination reconstruction formula, the concentration field data matrix of the un-simulated working conditions can be obtained.

[0145] Furthermore, the concentration field data matrix (in column vector form) of the un-simulated working conditions is restored to the two-dimensional grid data format consistent with the original CFD simulation according to the spatial arrangement rule of the sampling lines described in Step 2 (500 sampling lines, 1000 sampling points per line), and the grid data is visualized in MATLAB through color mapping (such as the colormap function) to intuitively display the distribution characteristics of PM2.5 around the vehicle.

[0146] The accuracy of the reduced-order model is verified by calculating the error with the full-order model. As Figure 12 shown, the error between the prediction result and the true value is less than 15%. Among them, the true value is the concentration distribution result obtained from the CFD numerical simulation of the un-simulated working conditions through Ansys Fluent 2021R1 in Step 1. First, in CFD-Post, the concentration data at the position of a specified line on the flow field is extracted by drawing a line at a specified position and exported as a data file (such as CSV format). Then, the "Regression" function in the "Data Analysis" tool of Excel is used to analyze this set of data, and the regression equation is output. The regression equation is multiplied by 0.85 and 1.15 respectively to obtain the 15% error interval, and the accuracy of the model can be intuitively evaluated through the regression curve graph.

[0147] Through the above process, a multi-output LS-SVM model can be trained to complete the non-linear mapping from the input features to the output modal coefficients.

[0148] Through the CFD-POD coupling method that can handle CFD calculations with moving mesh models, in the case of high fidelity, the calculation time of the PM2.5 concentration field near the moving vehicle model is reduced to 1 / 1000 of the original time.

[0149] Using the method for reduced-order calculation of the PM2.5 concentration field near a moving vehicle of the present invention, the PM2.5 concentration field data in a special area of the multi-functional complex tunnel shown in CN202310796176.1 can be simulated, and the required ventilation volume of the area can be obtained based on the PM2.5 concentration field data in subsequent calculations.

[0150] Although the embodiments of the present invention have been shown and described, for those of ordinary skill in the art, it can be understood that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

[0151] References:

[0152] [1] Wang, H., & Chen, Q. (2012). A new empirical model for predicting single-sided, wind-driven natural ventilation in buildings. Energy and Buildings, 54, 386 - 394.

[0153] [2] Mao, Q., Jin, S., Wu, X., et al. Dynamic variation characteristics of direct solar radiation intensity in Wuhan area under typical weather conditions [J]. Low Temperature Architecture Technology, 2018, 40(11): 108 - 111.

[0154] [3] Li, Y. Simulation study on the emission and diffusion of PM2.5 at the rear of automobiles [D]. Chang'an University, 2015.

[0155] [4] Chen, X., Zhong, W., & Li, T. (2023). Fast prediction of temperature and chemical species distributions in pulverized coal boiler using POD reduced-order modeling for CFD. Energy.

Claims

1. A method for predicting PM2.5 concentration field in a road area, characterized in that: The steps include: Step S1: using numerical simulation software FLUENT to calculate the PM2.5 concentration field data near the vehicle body on the simulated road area under several different working conditions within a certain period of time, and in the post-processing software CFD-POST, using multiple groups of sampling lines to extract the PM2.5 concentration field data of each working condition, the sampling lines are parallel to each other and of equal length, the interval between two adjacent sampling lines is less than 0.05m, each sampling line is provided with a number of sampling points, the spacing between the sampling points is less than 0.05m, and the surface formed by all sampling lines completely covers the road area; Step S2: For each working condition, the data sequences on all sampling lines are exported into CSV format files by the post-processing software CFD-POST, and through data conversion, the data sequences on all sampling lines of each working condition are combined into a data column vector: The data conversion includes: marking the data in the data sequence on each sampling line with a serial number: "op", indicating the "p"th sampling point in the sampling line with the initial number "o", so that the data sequences on all sampling lines of each working condition are arranged in the data column vector from top to bottom in ascending order of serial numbers; Step S3: by repeating step S2, a data column vector consisting of data sequences on all sampling lines under several different working conditions is obtained, and all data column vectors are horizontally combined to obtain a matrix containing PM2.5 concentration field data near the vehicle body on the simulated road area under different working conditions, which is named as a snapshot matrix; Step S4: Import the snapshot matrix into MATLAB in the form of a numerical matrix, and perform POD reduction calculation on the snapshot matrix to obtain the dominant mode, including: Construct snapshot matrix S = [x1, x2, ..., x m ], where S is an n×m matrix, each column is a column vector, n is the dimension of the column vector, and m is the number of column vectors; Use the built-in MATLAB function SVD to perform singular value decomposition on the snapshot matrix S in step 51 to obtain three matrices U, Σ, and V T ; S=UΣV T U is an n×n orthogonal matrix, in which the column vectors u1,u2,...,u n That is the mode; Σ is an n×m diagonal matrix with diagonal elements σ1,σ2,...,σ r are singular values; V T is an m×m orthogonal matrix containing the weights of the modes in terms of time or operating conditions; Extract the first r modes with larger singular values ​​from the matrix U as the dominant modes: U r =[u1,u2,...,u r ] The energy of the r modes needs to be greater than 90% of the total modal energy, and the singular value σ d The square of σ d 2 represents the energy contribution of the dth mode; The modal coefficient matrix A is calculated based on the product of Σ and V, A = Σ·V T ; The modal coefficient a of the dth mode in the modal coefficient matrix A d =σ d ·v d , v d is the matrix V T The corresponding column vector in ; Step S5: Use the least squares support vector machine model LS-SVM to construct a mapping relationship between the modal coefficients and the working condition parameters, and then obtain the modal coefficients of the unsimulated working condition; Step S6: The modal coefficient A of the unsimulated working condition obtained in step S5 r =[a1,a2,...,a r ], the modal coefficients of the unsimulated working condition are brought into the simplified linear combination reconstruction formula to obtain the concentration field data matrix of the unsimulated working condition: c=c0+u1a1+u2a2+…+u r a r ; c0 is the average concentration value of all concentration fields in step S1.

2. The method according to claim 1, characterized in that In step S1, the method for simulating different working conditions is: Step 11. When using CFD method to perform numerical simulation to obtain flow field data, in order to simulate the wind profile at the velocity inlet, the size of the ambient wind is set to: Among them, U z Represents the speed at the Z height from the road surface, U ref Indicates the reference height z ref The wind speed depends on the local meteorological station's wind speed on the open and flat ground. ref = the average wind speed observed at a height of 10m within 10 minutes; α is the ground roughness index, which is 0.22; Step 12. Activate the DO radiation model and the solar ray tracing method to simulate solar radiation at different times of the day. In the solar radiation calculator, set the time zone, longitude and latitude of the experimental site, and sunshine time to determine the direct and diffuse radiation intensities at the selected time. Step 13. The process of particle diffusion is simulated by a discrete term model in the numerical simulation software FLUENT, and the incident parameters of PM2.5 particles are set in the model, including incident flow rate, injection velocity, and incident duration; Among them, the incident flow is obtained according to the motor vehicle emission factor EF, and the PM2.5 mass flow Q of a vehicle that takes x seconds to travel 1 km is expressed as: The exhaust emission rates at different driving speeds refer to the vehicle exhaust collection test results; The temperature of the incident particles is 380K, and the duration must be greater than the moving grid movement time; Step 14. To simulate the vehicle motion, activate the Dynamic meshing model in the CFD model.

3. The method according to claim 1, characterized in that In step S5, the least squares support vector machine model LS-SVM is used to construct a mapping relationship between the modal coefficients and the operating condition parameters, and then the modal coefficients of the unsimulated operating condition are obtained, including: Optimize the weight vector of the least squares support vector machine algorithm and construct a multi-output regression model using the weight vector, Gaussian radial basis kernel function and error vector; The r modal coefficients obtained in step S4 are used to construct a training set and input into a multi-output regression model for processing, and the predicted values ​​of the r modal coefficients are output.

Images