A three-dimensional lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism

Abstract

The application discloses a three-dimensional Lagrange water droplet impact characteristic numerical simulation method based on MPI task parallelism, and comprises the following steps: a point cloud matching module based on KD-tree is built to realize fast pairing of three-dimensional flow field information; an MPI-based initialization computing environment is built; a Runge-Kutta method is used to iteratively solve and calculate water droplet motion trajectories; and a Monte Carlo method is used to calculate water collection coefficients of object surfaces. The application can greatly improve the calculation speed of three-dimensional Lagrange water droplet impact characteristic solving while ensuring the correctness and calculation precision of the calculation result.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of icing numerical simulation, and particularly relates to a three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism. BACKGROUND

[0002] When an aircraft passes through a cloud layer containing supercooled water droplets, the water droplets will impact and freeze on the windward surface of the aircraft, thereby forming an ice layer. The icing on the surface of the wings can significantly reduce the aerodynamic performance of the aircraft, thereby posing a serious threat to flight safety. The freezing, deformation and fragmentation of supercooled large droplets (SLD), as well as the splashing and rebound characteristics after impact, are significantly different from small droplets, and the abnormal icing of ice crystals in the engine inlet, which makes the icing phenomenon more complex and serious. Therefore, it is of great importance to carry out numerical simulation research on aircraft icing to ensure the safety of aircraft flight.

[0003] The numerical simulation of icing on aircraft mainly consists of three steps: flow field calculation, water droplet field calculation and ice shape calculation. The calculation of water droplet impact characteristics is a key factor affecting the results of icing calculation, and Eulerian method or Lagrangian method is usually used for calculation. Eulerian method regards water droplets as a continuous item. It is easier to realize large-scale parallel computing and three-dimensional water droplet impact characteristic calculation because it can share the same grid and calculation framework with flow field calculation, and has been applied more maturely. The main research results include: the calculation of water droplet impact characteristics on three-dimensional complex surfaces such as multi-element airfoils and engine inlets based on Eulerian two-phase flow method; the calculation of three-dimensional NACA0012 airfoil model using Eulerian method by pre-correction and four-step Runge-Kutta method to discretely solve water droplet motion equation, and the comparison with experimental results; the calculation of ice accumulation on UH-1H rotor blades based on Eulerian method of three-dimensional water droplet impact characteristics on rotating components; the calculation of water droplet impact characteristics on rotating components using rotating coordinate system method, and the ice accumulation on rotating components obtained by combining the Messinger icing model considering the influence of water film flow of centrifugal force; the calculation of water collection coefficients on engine nacelle with average volume diameter (MVD) of 80 μm and 250 μm using Eulerian method, and the icing on engine nacelle obtained by Messinger icing model based on this; the study of the influence of jet on water droplet impact on engine inlet components, the calculation of three-dimensional open double-hole cylinder, and the analysis of the influence of double-hole jet on water droplet impact characteristics, etc. However, Eulerian method has defects in dealing with the breaking problem of particle motion and the movement of ice crystals. With the Lagrangian method focusing on the particle itself, it is easy to deal with the problems of large water droplet motion deformation and breaking, as well as ice crystal motion and breaking, which has gradually become a research hotspot. Some researchers use OpenMP parallel combined with neural network acceleration method to calculate the water droplet collection coefficient of GLC305 swept wing and engine nacelle model. Compared with BP neural network, this method reduces the calculation time by 26.8% at most, and the water droplet impact characteristic calculation takes 2546 seconds and 3896 seconds respectively. Some researchers use serial code combined with constrained radial basis function (RBF) interpolation and adaptive grid refinement to control the number of released water droplets to the order of thousands to tens of thousands; water droplet collection coefficient calculation and verification are carried out on sphere, cylinder, three-element airfoil, M6 airfoil and NACA23012 airfoil with ice angle, etc. Although the calculation results are ideal, the calculation results of this method are strongly dependent on interpolation technology. Reducing the number of released water droplets can indeed speed up the overall calculation speed, but it does not conform to the physical reality.Some researchers have proposed a method aiming at inheriting the distributed memory parallel framework used in the flow field calculation, the core idea of which is to roughly calculate the water droplet trajectory first, then to weight the grid cells passed by the water droplets, and then to redivide the grid to ensure that the area frequently crossed by the water droplets can be allocated to as many computing cores as possible. Although this method speeds up the calculation of the Lagrangian method to some extent, it has obvious limitations: even after the grid weighting and redivision, a large number of computing cores are still in idle state. In addition, even if some areas frequently passed by the water droplets are allocated to multiple computing cores, when there is no water droplet passing through the grid area responsible by some core modules, these core modules are still idle and cannot participate in the calculation.

[0004] The above calculation method cannot completely utilize the advantages of parallel computing to make each core participate in the calculation, so the calculation speed is not greatly improved. SUMMARY

[0005] The purpose of the present application is to provide a three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism, which can ensure that the front-end flow field information is not distorted, quickly solve the trajectories of hundreds of thousands of particles, obtain the water collection coefficient of the surface, and close to the physical fact while ensuring the calculation accuracy and improving the calculation efficiency.

[0006] To achieve the above technical purposes, the technical scheme adopted by the present application is:

[0007] A three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism, the method comprising the following steps:

[0008] Step A, initializing the point calculation environment based on MPI, each core reading a set of cgns grid information and a set of front-end flow field.dat files generated by tecplot software;

[0009] Step B, calling the C++ nanoflann library to build a point cloud matching module based on KD-tree, matching the read-in grid coordinate points with the input front-end flow field information one by one, and then assigning the flow field information to the read-in grid nodes, to quickly match the three-dimensional flow field information;

[0010] Step C, initializing the spatial particle information, iteratively solving the water droplet motion trajectory based on the Runge-Kutta method, obtaining the coordinates of each particle at different times, finding the host cell of each particle at each time according to the directional lookup method, judging whether the particle crosses the current grid and whether it hits the wall surface according to the cross discrimination method, and updating the host cell of the particle or solving the position of the water droplet at the next time according to the judgment result until all particles hit the wall surface; using the MPI_Allreduce function to summarize the number of particles collected by all surface cells;

[0011] Step D, calculating the water collection coefficient of the object surface according to the Monte Carlo method.

[0012] In step B, the process of fast matching the three-dimensional flow field information includes the following steps:

[0013] For a two-dimensional data set, a binary tree construction method is used to first divide according to the first dimension feature, and after the division is completed, the nodes with first dimension data greater than the first data of the root node are divided into the right subtree of the root node, and the nodes with first dimension data less than the root node are divided into the left subtree; then recursion is continued, and the second dimension data is divided on the left and right subtrees to construct a KD-Tree.

[0014] The nearest neighbor point is searched through binary search-backtracking; specifically, binary search is used to search from the root node downward until a leaf node is found, during which the visited nodes are added to the stack and the shortest distance is recorded; after the leaf node is found, backtracking is started, and the previously visited nodes are popped out from the stack one by one, and it is judged whether the hypersphere with the query point as the center and the current shortest distance as the radius intersects with the partition surface, if there is an intersection, the other branch of the node is entered, and binary search is continued until the leaf node is searched; the cycle is repeated until there is no intersection between the hypersphere and the partition surface.

[0015] In step C, the process of initializing the spatial particle information includes the following steps:

[0016] The total number of particles to be released and the total number of iterations are set, the number of initialized spatial particle coordinates is the total number of particles to be released divided by the total number of iterations and then divided by the number of calculation cores; each core passes through the C++ interface to obtain the coordinates of the particles to be released in the current iteration step. <random>The std::uniform_real_distribution class in the library generates the spatial coordinates of the corresponding number of particles in a plane in front of the model; initializes the host cell where the particles are located, and initializes the particle velocity to the velocity of the incoming flow from a distance.

[0017] Step C, the process of iteratively calculating the trajectory of the water droplet based on the Runge-Kutta method, includes the following steps:

[0018] The Lagrange equations for the motion of the water droplet are as follows:

[0019] ;

[0020] in , These are the velocity vectors of the air and the water droplets, respectively. and The densities of air and water droplets are respectively. The drag coefficient, The Reynolds number is... The dynamic viscosity of air. The diameter of the water droplet is the average volume; g represents the acceleration due to gravity.

[0021] The Lagrange equations of motion for a water droplet can be expressed as a first-order ordinary differential equation for the droplet's velocity, and their expression is:

[0022] ;

[0023] In the formula, t represents time and u represents the velocity vector;

[0024] Solving the first-order ordinary differential equation using the fourth-order Runge-Kutta method yields:

[0025] ;

[0026] In the formula, For time step;

[0027] The coefficients are expressed as follows:

[0028] ;

[0029] ;

[0030] ;

[0031] ;

[0032] In obtaining the speed of water droplet movement After that, I received The position of the water droplet at that moment is:

[0033] .

[0034] Step C further comprises:

[0035] Step C1, obtaining the particle next time position coordinate;

[0036] Step C2, connecting the particle n time particle host unit center and n+1 time particle position;

[0037] Step C3, judging whether the line segment intersects with the edge or face of the n time particle host unit by the intersection judgment method, if there is no intersection, the particle is still in the current host unit, the host unit of the particle is updated, and the step C1 is returned; otherwise, if there is intersection, it is transferred to step C4;

[0038] Step C4, judging whether the crossed edge or face is a wall surface, if it is a wall surface, it is determined that the particle hits the wall surface and exits the loop; if the crossed edge or face is not a wall surface, the host unit is updated to the unit number which is based on the crossed edge or face and has the same edge or face as the n time host unit, the n+1 time particle position is connected from the center of the new host unit, and the step C1 is returned.

[0039] Step D further comprises:

[0040] Suppose that the particle concentration in the initial position of the remote water droplet release is , the particle concentration of the wall surface impact area is When the number of particles released at the initial position of the water droplet release is enough, the local water droplet collection rate of the wall surface impact area is calculated according to the following formula:

[0041] ;

[0042] Wherein is the water collection coefficient of the wall surface, is the number of particles impacting on the wall surface, is the number of particles released remotely; is the area of the wall surface unit, is the area of the remote release particle plane, is the equivalent mass of the water droplet particle.

[0043] Compared with the prior art, the beneficial effects of the present application are as follows:

[0044] The purpose of the present application is to provide a three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism, which can ensure that the front flow field information is not distorted, quickly solve the trajectories of hundreds of thousands of particles, obtain the water collection coefficient of the surface, thereby greatly improving the calculation speed of the three-dimensional Lagrangian water droplet impact characteristic solution, and also ensuring the correctness and calculation accuracy of the calculation result. BRIEF DESCRIPTION OF DRAWINGS

[0045] Figure 1 The three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism provided by the embodiment of the application;

[0046] Figure 2 The non-structure grid host unit search process schematic diagram provided by the embodiment of the application;

[0047] Figure 3 The original two-dimensional data set distribution schematic diagram;

[0048] Figure 4 The data distribution schematic diagram after the whole point set is divided into left and right two parts according to the first dimension feature (x-axis feature) and (5, 2) is selected as the root node;

[0049] Figure 5 The data distribution schematic diagram after the left and right two planes are divided into upper and lower two planes according to the second dimension feature (y-axis feature) and (4, 5) and (7, 6) are respectively selected as the root nodes;

[0050] Figure 6 The final KD-Tree structure schematic diagram;

[0051] Figure 7 The original data distribution diagram when the point to be queried is (6, 4.5);

[0052] Figure 8 The data distribution diagram after the first dimension feature of the root node is compared;

[0053] Figure 9 The data distribution diagram after the second dimension feature of the node on the right half plane split line is compared;

[0054] Figure 10 The binary search result schematic diagram;

[0055] Figure 11 The particle trajectory partition schematic diagram provided by the embodiment of the application;

[0056] Figure 12 The calculation flowchart of solving the particle trajectory by using the Lagrangian method provided by the embodiment of the application;

[0057] Figure 13 The solving schematic diagram of the water collection coefficient based on the Monte Carlo method provided by the embodiment of the application;

[0058] Figure 14 The calculation time and acceleration ratio change curve diagram along with the increase of the number of calculation cores provided by the embodiment of the application;

[0059] Figure 15 A schematic diagram of water droplet motion trajectory solved by the NACA0012 wing segment provided by the embodiment of the present application is shown in the figure;

[0060] Figure 16 A water collection coefficient cloud chart of the NACA0012 wing segment surface provided by the embodiment of the present application is shown in the figure;

[0061] Figure 17 A comparison chart of water collection coefficient curve and test value of the NACA0012 wing segment cross section provided by the embodiment of the present application is shown in the figure;

[0062] Figure 18 A schematic diagram of water droplet motion trajectory of the cylinder flow provided by the embodiment of the present application is shown in the figure;

[0063] Figure 19 A water collection coefficient cloud chart of the cylinder flow provided by the embodiment of the present application is shown in the figure;

[0064] Figure 20 A comparison chart of water collection coefficient curve and test value of the cylinder flow cross section provided by the embodiment of the present application is shown in the figure;

[0065] Figure 21 A schematic diagram of water droplet motion trajectory of the round ball flow provided by the embodiment of the present application is shown in the figure;

[0066] Figure 22 A water collection coefficient chart of the round ball surface provided by the embodiment of the present application is shown in the figure;

[0067] Figure 23 A comparison chart of water collection coefficient curve and test value of the round ball cross section provided by the embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0068] The embodiments of the present application are further described in detail below with reference to the accompanying drawings.

[0069] Referring to Figure 1 The present application discloses a three-dimensional Lagrangian water droplet impact characteristic numerical simulation method based on MPI task parallelism, which comprises the following steps:

[0070] Step A: initializing a point calculation environment based on MPI, each core reads a set of cgns grid information and a set of.dat files of front flow field generated by tecplot software;

[0071] Step B: calling C++nanoflann library to build a point cloud matching module based on KD-tree, matching the read-in grid coordinate points with the input front flow field information one by one, and then assigning the flow field information to the read-in grid nodes, so as to quickly match the three-dimensional flow field information;

[0072] Step C, initializing space particle information, based on Runge-Kutta method iterative solution to calculate water droplet trajectory, the coordinate of each particle at different time is obtained, the host unit of each particle at each time is found according to the directional search method, whether the particle crosses the current grid and whether it hits the wall is judged according to the intersection discrimination method, and the host unit of the particle is updated or the position of the water droplet at the next time is solved according to the judgment result, until all particles hit the wall; the particle number collected by all surface elements is summarized by using the MPI_Allreduce function;

[0073] Step D, calculating the water collection coefficient of the object surface according to the Monte Carlo method.

[0074] Step A specifically comprises: opening the front flow field by using tecplot software, and then outputting the.dat file of the flow field by using the tecplot software; initializing the MPI environment by using the MPI_Init(), MPI_Comm_Rank() and MPI_Comm_size() functions; and reading the CGNS grid by each core through the CGNS library function.

[0075] After initializing the MPI environment, each core reads a set of CGNS grid information and a set of.dat flow field information generated by the tecplot software, calls the point cloud matching module, quickly matches the read-in grid coordinate points with the input front flow field information, and then assigns the flow field information to the read-in grid nodes, so as to realize the quick pairing of the three-dimensional flow field information. The following will illustrate the specific pairing principle of the point cloud matching module based on KD-tree of the application by an example.

[0076] Figures 3 to 10 The point cloud matching method based on KD-tree provided by the embodiment of the application is shown in the schematic diagram. Figures 3 to 6 The division process of the original two-dimensional scatter plot according to the binary search method is gradually shown, Figures 7 to 10 The search flowchart of quickly finding the nearest neighbor point of the target point according to the nearest neighbor search method is shown.

[0077] Suppose that there is a two-dimensional data set composed of the following seven two-dimensional data:

[0078] ;

[0079] Firstly, the first dimension feature is divided, the process is consistent with the binary tree construction method, after the division is completed, the nodes with the first dimension data greater than the first bit data of the root node are divided into the right subtree of the root node, and the nodes with the first dimension data less than the root node are divided into the left subtree; the recursion is continued, and the second dimension data is divided on the left and right subtrees. As shown in Figures 3-6 . Figure 3 is the original data distribution, Figure 4 is the first time to divide according to the first dimension feature (x-axis feature), and the whole point set is divided into left and right two parts by selecting (5, 2) as the root node. Figure 5 is the second time to divide according to the second dimension feature (y-axis feature), and the left and right two planes are divided into upper and lower two planes by selecting (4, 5) and (7, 6) as the root nodes. Figure 6 is the final KD-Tree structure, and (1, 3), (2, 8), (8, 7) and (9, 4) are leaf nodes.

[0080] After the KD-Tree is constructed, the nearest neighbor point is searched by a binary search-backtracking method. Firstly, a binary search method is used to search from the root node downwards until a leaf node is found, during which the visited nodes are added to a stack, and the shortest distance is recorded. After the leaf node is found, backtracking is started, and the previously visited nodes are popped out from the stack one by one. It is judged whether the hyper-sphere with the query point as the center and the current shortest distance as the radius intersects with the split plane. If there is an intersection, the other branch of the node is entered, and the binary search is continued until the leaf node is searched. The cycle is repeated until there is no intersection between the hyper-sphere and the split plane.

[0081] For example, it is assumed that the query point is (6, 4.5), and the original data distribution is shown in Figure 7 The binary search is performed, and the first dimension feature of the root node is compared. Since 6>5, the right child is entered, as shown in Figure 8 The distance between the query point and the root node is calculated, and dist≈2.69 is obtained as the current shortest distance. The binary search is performed on the second dimension feature of the nodes on the split line of the right half plane, and 4.5<6, so the left child of the node is entered, as shown in Figure 9 The distance between the query point and the node (7, 6) is calculated, and dist≈1.80 is obtained, which is smaller than the current shortest distance, so the current shortest distance is updated. The circle with the query point as the center and the current shortest distance as the radius has no intersection with other split lines, as shown in Figure 10 Therefore, the node (7, 6) is considered to be the nearest neighbor of the query point, and the shortest distance is 1.80. The search is completed. In the actual grid and flow field information matching process, the searched point is the point in the point set X, not the point outside the point set. Therefore, the nearest neighbor of the searched node is the node itself. The flow field information of the node is assigned to the node information of the grid.

[0082] In step C, the process of initializing the spatial particle information includes: initializing the spatial particle coordinates of each core at a far distance from the incoming flow; the number of initialized spatial particle coordinates is the total number of planned released particles divided by the total number of iteration steps and then divided by the number of calculation cores; initializing the host unit of the particle, and initializing the particle velocity as the far distance incoming flow velocity.

[0083] Specifically, first, the total number of particles to be released and the total number of steps to be iterated are set, and the number of initialized spatial particle coordinates is the total number of particles to be released divided by the total number of iteration steps and then divided by the number of computing cores. Each core passes through C++ <random>The std::uniform_real_distribution class of the library generates the spatial coordinates of the corresponding number of particles at a certain plane in front of the model (for example, at 4-8 times the characteristic scale in front of the model). As shown in FIG. 1, a certain amount of particles are randomly released at a certain plane in front of the model. The host cell in which the particles are located is initialized, and the particle velocity is initialized as the far-field incoming flow velocity. Figure 11

[0084] Secondly, the particle motion equation is solved by using the four-step Runge-Kutta method to obtain the coordinates of each particle at different times.

[0085] Specifically, the equation of particle motion is:

[0086] ;

[0087] wherein , are the velocity vectors of air and water droplets respectively, are the densities of air and water droplets respectively, Cd is the drag coefficient, Re is the Reynolds number, is the dynamic viscosity of air, and D is the average volume diameter of water droplets.

[0088] The Lagrangian water droplet motion equation is expressed as a first-order ordinary differential equation of the water droplet velocity, and the expression is:

[0089] ;

[0090] The four-order Runge-Kutta method is used for solving:

[0091] ;

[0092] wherein each coefficient can be expressed as:

[0093] ;

[0094] ;

[0095] ;

[0096] ;

[0097] After obtaining the water droplet motion velocity , the water droplet position at time t can be expressed as:

[0098] .

[0099] ​​Thirdly, the host unit of each particle at each time is found according to the directional search method, whether the particle crosses the current grid and whether the particle hits the wall are judged according to the intersection discrimination method, and the host unit of the particle is updated or the position of the water droplet at the next time is solved according to the judgment result, until all particles hit the wall.

[0100] The calculation process of solving the water droplet trajectory by the Lagrangian method is shown in Figure 12 . After obtaining the position coordinates of the particle at the next time, the directional search method is used to find the host unit of the particle at the next time. As shown in Figure 2 , c1, c2, c3, c4, c5, and c6 are the centers of the grid units, Pn and Pn+1 are the position coordinates of the water droplet at the n time and the n+1 time respectively, and the dashed lines are the connecting lines from Pn, c1, c2, c3, c4, and c5 to the point Pn+1 respectively; the intersection points of the dashed lines and the grid boundaries are in the box. The process of directional search includes the following steps: connecting the particle host unit at the n time and the particle position at the n+1 time, judging whether the line segment intersects with the edge (face) of the particle host unit at the n time by the intersection discrimination method, if there is no intersection, the particle is still in the current host unit, the host unit of the particle is updated, if there is an intersection, judging whether the crossed edge (face) is a wall, if it is a wall, it is determined that the particle hits the wall and exits the loop, if the crossed edge (face) is not a wall, updating the host unit to the unit number which has the same edge (face) as the n time host unit based on the crossed edge (face), connecting the particle position at the n+1 time from the center of the new host unit, and judging the intersection of the line segment c1-Pn+1 and the edge (face) of the new host unit by the intersection discrimination method again, and so on until the host unit of the water droplet at the n+1 time is found.

[0101] Finally, the MPI_Allreduce function is used to summarize the number of particles collected by all surface units. Since the number of particles calculated by each core is the result after the calculation task is allocated, the particles of different cores may hit the same grid, and the MPI_Allreduce reduction function is needed to accumulate the number of particles collected in each grid, so as to facilitate the subsequent solution of the wall water collection coefficient.

[0102] In step D, the wall water collection coefficient is calculated by the Monte Carlo method. The specific method of solving the water collection coefficient based on the Monte Carlo method is shown in Figure 13 . It is assumed that the initial position of the water droplet released far away is the square area on the left, the particle concentration in the square area is n0 / A0, and the curved area on the right is the wall hitting area, the particle concentration is n wall / A wall . When the number of particles released at the initial position of the water droplet is sufficient, the local water droplet collection rate of the wall hitting area can be calculated as follows:

[0103] .

[0104] wherein is the water collection coefficient of the wall surface, are the number of particles impinging on the wall surface, respectively, and the number of particles released at a distance; are the area of the wall surface element, respectively, and the area of the plane of the particles released at a distance, is the equivalent mass of the water droplet particles.

[0105] Several examples are given below as specific measures of the method disclosed in the present application.

[0106] Example 1, in order to test the advantage of the method disclosed in the present application in improving the calculation efficiency, the water collection coefficient distribution on the NACA0012 wing section is calculated according to the calculation case in the literature with different core numbers, the computer configuration is tested and the calculation results are compared with the Lewice icing calculation software. The chord length of the wing section is 0.5334, the span length is 1 / 10 of the chord length. The attack angle is 4°, the Mach number of the incoming flow is 0.32, the air density is 1.097 kg / m 3 , the static pressure is 89867 Pa, the particle size (MVD) is 40 μm, and the LWC is 1.0 g / m 3 . Figure 14 The calculation time and speed-up ratio curves under the calculation core number of 1, 4, 8, 12, 24, 48, and 60 are given. With the increase of the calculation core number, the calculation time required decreases obviously, and when the core number is 48, the calculation time required is only 132.04 seconds. With the increase of the calculation core number, the speed-up ratio shows a linear growth trend. When the core number is 60, the speed-up ratio reaches 29.78, which highlights the outstanding advantage of the method disclosed in the present application in improving the calculation efficiency. Figure 15 is the water droplet trajectory curve obtained by the calculation of this example. Figure 16 is the wing section surface water collection coefficient cloud chart obtained by the calculation of this example. Figure 17 is the comparison result of the water collection coefficient distribution curve (corresponding to calculate) obtained by the calculation of the wing section spanwise cross section of this example and the calculation result (corresponding to Lewice) of the Lewice calculation software, wherein the vertical axis Beta represents the water collection coefficient, which proves the accuracy of the water collection coefficient calculation obtained by the calculation of this example.

[0107] Example 2, this example is a three-dimensional cylindrical disturbance calculation example, the diameter of the cylinder is 101.6 mm, the incoming flow velocity is 80 m / s, the air density is 1.097 kg / m 3 , the static pressure is 89867 Pa, the attack angle is 0°, the particle size (MVD) is 16 μm, Figure 18 The distribution of part of the trajectories of this example is given, Figure 19 The surface water collection coefficient cloud chart of the embodiment is given, due to the existence of computer errors, in order to ensure that the two end particles can be smoothly found into the host unit, when setting the particle release area, the release limit in the Z direction is set to 0.01 times and 0.99 times of the total length in the Z direction respectively, so there is a certain error on both sides, but the whole 3D cylindrical disturbance is actually a result of stretching in the Z direction, so Figure 20 The water collection coefficient distribution curve (corresponding to calculate) at the section of Z=0.5 times the total length in the Z direction is compared with the test value (corresponding to EXP), due to the fact that the icing wind tunnel test itself has many parameters, and the repeated test results also have a certain deviation, wherein the black solid line part is the collection rate range of different test batches under this working condition.

[0108] In embodiment three, the embodiment is to further verify the calculation accuracy of the water droplet impact characteristics of the disclosed method in pure three dimensions. The calculation working condition of the circular ball disturbance is that the diameter of the circular ball is 15.04 cm, the air density is 1.097 kg / m 3 , the liquid water content (LWC) is 1 g / m 3 , the incoming flow velocity is 75 m / s, the static pressure is 95840 Pa, and the water droplet diameter is 11.5 μm. Figure 21 is the three-dimensional circular ball disturbance part water droplet trajectory, Figure 22 is the three-dimensional circular ball disturbance surface water collection coefficient distribution cloud chart, Figure 23 is the comparison of the water collection coefficient distribution graph (corresponding to calculate) at the section of Z=0 of the three-dimensional cylindrical disturbance and the test value (corresponding to EXP). The water droplet impact limit and the maximum water collection coefficient Beta are in good agreement with the test results.

[0109] In summary, the disclosed method can greatly improve the calculation speed of the three-dimensional Lagrangian water droplet impact characteristics, while ensuring the correctness and calculation accuracy of the calculation results.

[0110] Although the preferred embodiments of the present application have been described, those skilled in the art can make further changes and modifications to the embodiments once they know the basic inventive concept. Therefore, the appended claims are intended to be interpreted as including all changes and modifications falling within the scope of the present application.

[0111] Obviously, those skilled in the art can make various modifications and variations to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalents, the present application also intends to include these modifications and variations.< / random> < / random>

Claims

1. A method for numerical simulation of the characteristics of three-dimensional Lagrangian water droplet impingement based on MPI task parallelism, characterized in that, The method comprises the following steps: Step A: calculating the environment based on MPI initialization points, each core reading a set of cgns grid information and a set of.dat files of the front flow field generated by tecplot software; Step B: calling the C++ nanoflann library to build a point cloud matching module based on KD-tree, matching the read-in grid coordinate points with the input front flow field information one by one, and then assigning the flow field information to the read-in grid nodes to quickly match the three-dimensional flow field information; Step C: initializing the spatial particle information, iteratively solving and calculating the water droplet motion trajectory based on the Runge-Kutta method, obtaining the coordinates of each particle at different times, finding the host cell of each particle at each time according to the directional lookup method, judging whether the particle crosses the current grid and whether it collides with the wall surface according to the cross discrimination method, and updating the host cell of the particle or solving the position of the water droplet at the next time according to the judgment result until all particles collide with the wall surface; using the MPI_Allreduce function to summarize the number of particles collected by all surface cells; Step D: calculating the water collection coefficient of the object surface according to the Monte Carlo method; In step B, the process of quickly matching the three-dimensional flow field information comprises the following steps: For a two-dimensional data set, a binary tree construction method is adopted, first, the data is divided according to the first dimension feature, after the division is completed, the nodes with data greater than the first data of the root node are divided into the right subtree of the root node, and the nodes with data smaller than the root node are divided into the left subtree; then, the recursion is continued, and the data is divided according to the second dimension data on the left and right subtrees to construct a KD-Tree; The nearest neighbor is searched through binary search-backtracking; specifically, binary search is used to search from the root node downward until a leaf node is found, during which the visited nodes are added to the stack, and the shortest distance is recorded; after the leaf node is found, backtracking is started, and the previously visited nodes are popped out from the stack one by one, it is judged whether the hypersphere with the query point as the center and the current shortest distance as the radius intersects with the partition surface, if there is intersection, the other branch of the node is entered, and the binary search is continued, until the leaf node is searched; the cycle is repeated until there is no intersection between the hypersphere and the partition surface. In step C, the process of initializing the spatial particle information comprises the following steps: Set the total number of particles to be released and the total number of steps to be iterated, the number of initialized spatial particle coordinates is the total number of particles to be released divided by the total number of iteration steps and then divided by the number of computing cores; each core passes the C++ <random>The std::uniform_real_distribution class of the library generates the spatial coordinates of the corresponding number of particles at a plane in front of the model; the host cell of the particle is initialized, and the velocity of the particle is initialized as the far-field flow velocity.< / random> 2. The MPI task parallel based three-dimensional Lagrangian water droplet impact characteristic numerical simulation method according to claim 1, characterized in that, In step C, the process of iteratively solving and calculating the water droplet motion trajectory based on the Runge-Kutta method comprises the following steps: The Lagrangian water droplet motion equation is constructed as: where u a d are the velocity vectors of air and water droplets, respectively, p a d are the densities of air and water droplets, respectively, C d is the drag coefficient, Re is the Reynolds number, m a is the dynamic viscosity of air, and D is the average volume diameter of water droplets; g represents the gravitational acceleration;​​ The Lagrangian water droplet motion equation is expressed as a first-order ordinary differential equation of the water droplet velocity, and the expression is: In the formula, t represents time, and u represents velocity vector; The first-order ordinary differential equation is solved by using the fourth-order Runge-Kutta method, and the following formula is obtained: In the formula, Δt is the time step; Wherein, each coefficient is expressed as: After obtaining the water droplet motion speed the water droplet position at time t n+1 is obtained as 3. The MPI task parallel based three-dimensional Lagrangian water droplet impact characteristic numerical simulation method according to claim 1, characterized in that, Step C further comprises: Step C1: obtaining the position coordinates of the particle at the next time; Step C2, connecting the particle host cell center at time n and the particle position at time n+1; Step C3, judging whether the line segment intersects with the edge or face of the particle host cell at time n by cross-discrimination method, if there is no intersection, the particle is still in the current host cell, updating the host cell of the particle, returning to step C1; otherwise, if there is intersection, turning to step C4; Step C4, judging whether the intersected edge or face is a wall, if it is a wall, determining that the particle hits the wall and exiting the loop; if the intersected edge or face is not a wall, updating the host cell to the cell number which shares the edge or face with the host cell at time n based on the intersected edge or face, connecting the particle position at time n+1 from the center of the new host cell, and returning to step C1.

4. The MPI task parallel based three-dimensional Lagrangian water droplet impact characteristic numerical simulation method according to claim 1, characterized in that, Step D further comprises: Assume that the particle concentration in the initial position of the water droplet release is n0 / A0, and the particle concentration in the wall impact area is n wall / A wall When the number of particles released at the initial position of the water droplet release is sufficient, the local water droplet collection rate of the wall impact area is calculated according to the following formula: where β is the water collection coefficient of the wall, n wall is the number of particles impinging on the wall, n0is the number of particles released at a distance; A wall is the area of the wall element, A0is the area of the plane where the particles are released at a distance, m d is the equivalent mass of the water droplet.

Images