A high-precision fluid-induced vibration simulation method based on FETI

Abstract

The present invention discloses a high-precision flow-induced vibration simulation method based on FETI. First, the corresponding example data is read in, and the NewMark method is used to numerically discretize the flow-induced vibration dynamic process. Secondly, the FETI method is used to perform parallel decomposition on the discretized equations, and the subdomains are glued together by Lagrange multipliers at the partition boundaries. A graph bipartition algorithm for domain boundary balance is proposed to balance the number of elements and the computational amount in each subdomain, ensuring load balance among processes. Finally, the preconditioned conjugate gradient method is used for iterative solution. After iterative solution, according to the obtained corresponding displacements, the NewMark method is used to update the time steps of the flow-induced vibration process. The present invention has completed the solution of grid data on the scale of hundreds of millions, improved the efficiency of solving large-scale flow-induced vibration problems, realized the fast and efficient simulation of flow-induced vibration, and ensured load balance among processes.

Description

Technical Field

[0001] The present invention belongs to the field of numerical simulation implementation of reactor flow-induced vibration, and particularly relates to a high-precision simulation method for large-scale flow-induced vibration based on FETI. Background Art

[0002] The reactor core components are diverse and complex, and the computing resources and storage required for full-core high-fidelity structural mechanics simulation are difficult to achieve in a small cluster computing system. Relying on high-performance computing technology to achieve high-fidelity numerical calculations has become a necessary way to study reactor safety. Relying on supercomputers and using mathematical and physical calculation methods to perform refined numerical simulations on physical processes such as reactor core neutron physics, structural mechanics, and fuel performance has become a forefront research hotspot internationally.

[0003] Fuel rod assemblies operate under high temperature, high pressure, and strong irradiation, and are subject to the erosion of coolant for a long time. The cyclic flow of coolant will cause long-term micro-amplitude vibrations of the fuel rod assemblies, and these vibrations always exist along with the operation of the reactor. These vibrations cause relative displacements at the interfaces where the fuel rods contact the grids and cause cladding wear at the supports (grid-to-rod fretting, GTRF). This progressive wear of the fuel rods is a key factor affecting the structural integrity of the fuel assemblies under normal operation and accident conditions. According to the statistical data of the worldwide fuel rod failure rate, the GTRF wear failure accidents account for up to 78%. To ensure the integrity of the structure during its service life, it is only possible to control the structural vibration response caused by flow-induced vibration within an acceptable range by optimizing the structural design and adjusting the flow rate. Therefore, numerical analysis and calculation of flow-induced vibration of fuel assemblies are an essential and indispensable part of safety analysis.

[0004] High-precision analysis and simulation of the fluid-induced vibration of fuel rods, such as the high-precision numerical simulation of the wire-wrapped fuel rods in the fourth-generation fast reactor, often require generating mesh data with a scale of more than tens of millions of elements. When considering the turbulence model, the required mesh scale reaches hundreds of millions or even billions of mesh nodes, which poses severe challenges to the hardware configuration of the computer, the computational algorithm, and the spatio-temporal efficiency of its preprocessing process. Solving problems using serial algorithms becomes unacceptable in terms of memory space and time. Taking the open-source software Code_Aster as an example, for the fluid-induced vibration simulation of a single wire-wrapped fuel rod with one pitch (about 100,000 mesh elements), the single-core solution takes 2 hours. Estimated from this, for the serial solution of cases with hundreds of millions of elements, it becomes unacceptable in terms of memory space and time. Most current research focuses on small-scale cases and does not analyze cases with a scale of more than tens of millions of elements. Utilizing high-performance computers to develop efficient parallel multi-field numerical simulation solution algorithms and corresponding parallel pre- and post-processing algorithms is the only way to meet this challenge. Conducting research on parallel solution methods for fluid-induced vibration is the key to improving the solution efficiency of fluid-induced vibration.

[0005] The finite element method is an approximate numerical solution method for solving complex engineering problems and has now been widely applied to various disciplines such as mechanics, thermotics, and electromagnetics. The traditional finite element method generally solves problems in a serial manner, known as the serial finite element method. When faced with large and complex engineering problems, the serial finite element method often exhibits phenomena such as excessive solution time and high memory occupancy. With the advent of parallel computers that support parallel computing, three techniques suitable for parallelizing the serial finite element method have emerged, namely parallel solution of equations, the element-by-element (EBE) finite element method, and the finite element domain decomposition method. The finite element domain decomposition method based on the Schur Complement was first proposed by Farhat et al. and is called the Finite Element Tearing and Interconnecting method (FETI), which is commonly used to solve the numerical solution of large linear systems that appear in linear engineering problems. In the FETI method, a main body is decomposed into several non-overlapping subdomains, and the continuity between these subdomains is enforced by Lagrange multipliers. The FETI method uses duality theory, can derive a smaller and relatively well-conditioned dual problem, and effectively solves it through a suitable variant of the conjugate gradient algorithm. Summary of the Invention

[0006] To solve the above problems, the present invention proposes a high-precision fluid-induced vibration simulation method based on FETI. The present invention relates to the simulation of high-precision fluid-induced vibration using the NewMark integration method based on the FETI method and the domain boundary balance applicable to the FETI method. Figure 2The sub-algorithm divides the process of solving the high-precision fluid-induced vibration problem by the Newmark integration method of the FETI method into two stages: in the first stage, the large-scale parallel FETI method is used to divide the solution domain of hundreds of millions of grid cells into multiple non-overlapping sub-domains for parallel solution; in the second stage, the dynamic Newmark method is used to update the time step of the fluid-induced vibration process.

[0007] The specific steps of the method of the present invention are as follows:

[0008] Step 1. Read in the corresponding example data, mainly grid data, as well as material parameters, etc.

[0009] Step 2. Numerically discretize the fluid-induced vibration dynamics process using the Newmark method.

[0010] Step 3. Use the FETI method to perform parallel decomposition on the discretized equations, divide the sub-domains, and the sub-domains are glued together by the Lagrange multiplier λ at the division boundary and satisfy the corresponding boundary conditions. According to the characteristics of the FETI sub-domain division boundary, a sub-algorithm for domain boundary balance is proposed Figure 2 to balance the number of elements and the amount of calculation in each sub-domain, and ensure load balance between processes.

[0011] Step 4. Use the preconditioned conjugate gradient method for iterative solution to obtain the Lagrange multiplier λ through iterative solution. Obtain the corresponding displacement according to the Lagrange multiplier λ and record the displacement value.

[0012] Step 5. According to the displacement value, the Newmark method realizes the time step update of the fluid-induced vibration process.

[0013] The beneficial effects of the present invention: The present invention combines the FETI method with the Newmark integration method to complete the solution of grid data on the order of hundreds of millions, improves the efficiency of solving large-scale fluid-induced vibration problems, and realizes the fast and efficient simulation of fluid-induced vibration. Using the sub-algorithm for domain boundary balance Figure 2 for regional decomposition ensures load balance between processes and further improves the efficiency of solving large-scale fluid-induced vibration problems. Description of the Drawings

[0014] Figure 1 is the specific process of solving the high-precision fluid-induced vibration problem by the Newmark integration method based on the FETI method;

[0015] Figure 2 is the description of halo weights;

[0016] Figure 3 is the time ratio of a single iterative solution;

[0017] Figure 4 is the matrix-vector multiplication for multi-threaded parallel operation to solve;

[0018] Figure 5 is a weak scalability test. Specific implementation manner

[0019] The present invention will be further described below in conjunction with the accompanying drawings. Please refer to Figure 1 ; Figure 1 The specific process for solving the high-precision fluid-induced vibration problem of the present invention is given. It is mainly divided into two stages: the first stage is the FETI solution stage, which can be mainly divided into two steps. The first step is to use the domain boundary balance Figure 2 subdomain algorithm for subdomain division; the second step is to use the preconditioned conjugate gradient method for iterative solution. The second stage is to update the time step using the NewMark method.

[0020] The above implementation steps will be described in detail below:

[0021] Step 1: Take the CEFR single-box fuel rod assembly as an example for numerical simulation. The bottom of the fuel rod bundle is fixed, and the simulation acts on the outer surfaces of the fuel rod and the wire winding in the form of fluid pressure. Read the corresponding grid parameters distributedly. The grid file format is VTK format, mainly tetrahedral grids and hexahedral grids. Set the corresponding material parameters to be read in, mainly material density, elastic modulus, Poisson's ratio, time step, number of subdomains divided, etc. The material density in this simulation is 7890 kg / m 3 , Poisson's ratio is 0.3, elastic modulus is 2.1 * 10 5 MPa, simulate 12 time steps, each time step is 0.01 second, and the number of subdomains divided is 512 subdomains.

[0022] Step 2: The fluid-induced vibration problem is derived to obtain the overall dynamic equation by using the principle of virtual work for the equilibrium equation, mass matrix, damping matrix, geometric equation, physical equation, and displacement boundary conditions. Use a numerical algorithm (NewMark algorithm) to establish the dynamic response equation for the next moment based on the dynamic response values of the system at the current moment and the previous few moments, and solve it. The following overall dynamic equation is obtained for the fuel rod fluid-induced vibration problem modeling:

[0023]

[0024] Among them, M is the mass matrix, C is the damping matrix, K is the stiffness matrix, f(t) is the nodal load, u is the displacement vector of the fuel rod to be solved, the first derivative of u is the velocity, and the second derivative of u is the acceleration. The NewMark method is used for numerical discretization of this dynamic process. The NewMark method makes corrections for the linear acceleration assumption, and two parameters η and γ are introduced into the velocity and displacement expressions at the Δ t +t moment to obtain the discrete format of the NewMark algorithm:

[0025]

[0026]

[0027] Substituting equations (2) and (3) into the original second-order differential equation (1), we get:

[0028]

[0029] Simplify equation (4) into the following linear equation form for solving the unknown u n+1 , where n is time:

[0030] Au n+1 = H n (5)

[0031] where matrix A and matrix H n are respectively

[0032]

[0033]

[0034] Step 3: Use the FETI method to perform parallel decomposition on the discretized equation and divide the subdomains. When using the FETI method to divide the subdomains, use the sub-algorithm based on domain boundary balance proposed by the present invention Figure 2 to divide the solution domain into subdomains. The solution domain is decomposed into I independent non-overlapping subdomains. These subdomains are glued together by the Lagrange multiplier λ at the division boundary and satisfy the boundary conditions. Thus, formula 5 can be transformed into the following form:

[0035]

[0036] where,

[0037]

[0038] represents the displacement vector of each subdomain after decomposition. B I is a signed boolean matrix that defines the connectivity between subdomains and is composed of 1, -1, and 0. If subdomains X and Y are connected, set the coefficient from subdomain X to subdomain Y to 1 and the coefficient from subdomain Y to subdomain X to -1; if subdomains X and Y are not connected, set the coefficient to 0. O is the zero matrix.

[0039] The special boundary treatment method of the FETI method will cause serious unbalanced edge division and edge loss in the traditional point division method. The number of nodes on the subdomain boundary affects the iterative solution efficiency, and the computational cost of the constructed solution matrix increases with the increase of the edge computational cost. Therefore, it is crucial to ensure the balance of the number of nodes on each subdomain boundary. When performing domain decomposition in FETI, the Figure 2 sub-algorithm based on domain boundary balance can not only ensure the balance of internal nodes in each subdomain, but also ensure the balance of boundary nodes of this subdomain, thus improving the corresponding solution efficiency.

[0040] Based on the traditional greedy graph growth algorithm, two weights are added to the vertices: for the vertices belonging to the boundary nodes, the halo weight is set to 1, for the vertices that are not boundary nodes, the halo weight is set to 0, for the vertices of non-boundary nodes inside the subgraph, the halo weight is set to 1, and for the vertices belonging to the boundary nodes, the halo weight is set to 0. Add the vertices with halo weight of 1 to the halo point set V h . Please refer to Figure 2 , the black dots in the left figure are separators, which divide the initial graph into two subgraphs, and the halo weights are both 4. The right figure shows that the two subgraphs are divided into 4 subgraphs, and the halo weights of the subgraphs are 4, 4, 4, 4 respectively.

[0041] The following pseudocode gives the specific process of the Figure 2 sub-algorithm for domain boundary balance:

[0042]

[0043] First, select two seed vertices ω h and ω 1 in the halo point set V 2 , so that ω 1 and ω 2 maintain the distance as much as possible to ensure that ω 1 and ω 2 , the nodes with later growth meet. Divide the initial graph into two subgraphs, and the point sets of the two subgraphs are represented by P 1 and P 2 respectively. P 1 and P 2 are initially empty and grow from the seed vertices ω 1 and ω 2 . Starting from ω 1 and ω 2 , divide the remaining vertices respectively, and adjust according to the halo balance in the two point sets P 1 and P 2 . If the number of vertices with halo weight of 1 in P 1 (or P 2 ) is more than that in P 2(or P 1 ), then select the target vertex with halo weight 1 as much as possible to add to P 1 (or P 2 ). If there are more vertices with halo weight 1, vertices with non-halo weight 1 are preferred to join P 1 (or P 2 ). After considering the halo equilibrium, P 1 Select the distance 1 Recently and ω 2 The farthest vertex (dist(V,ω 1 )-dist(V,ω 2 ) minimum value vertex). This process needs to satisfy the following formula:

[0044]

[0045]

[0046] Where S represents the edge separator, V h Represents the set of halo points in the solution domain. Maintaining the minimum edge separator and halo node balance can ensure node balance at the subdomain boundary and load balance of edges between subdomains. Select different seed vertices for multiple growths and split according to the minimum edge separator and halo node balance.

[0047] As shown in Table 1, Table 2, and Table 3, it can be found that the domain boundary balance based on the present invention Figure 2 Comparison between the sub-domain division algorithm and the open source software metis. Table 1 divides a hexahedral mesh with a total number of units of 29026600. The Avg value in the table is the number of divisions under ideal conditions, and Max and Min represent the maximum number of meshes and the minimum number of meshes in the subdomain, respectively. It can be seen from Table 1 that the maximum and minimum number of subdomains divided by the method proposed in the present invention are closer to the ideal value. Table 2 divides a tetrahedral mesh with a total number of units of 102853760, which also proves that the method proposed in the present invention has a significantly better effect in sub-domain division. The time taken by this method to divide the tetrahedral mesh example into 512, 256, and 128 sub-domains is compared with the time taken by the open source graph division software metis. The time taken by the two is relatively close, and the sub-domain division efficiency is close.

[0048] Table 1 Comparison of domain division in hexahedral mesh example

[0049]

[0050] Table 2 Comparison of domain division in tetrahedral mesh example

[0051]

[0052] Table 3 Comparison of Sub - domain Partitioning Efficiency

[0053]

[0054] Step 4: Both matrix A and matrix B are non - singular matrices, so that the unknown λ has a unique solution. According to the generalized inverse of the matrix, we have:

[0055] u n+1 =A + (H n -B T λ)+Rα s.t.AR=O or R T A=O (12)

[0056] where the matrix R is the basis of the null space of matrix A, and A + is the pseudo - inverse matrix of A. Substituting the boundary condition Bu = 0 into equation (12), we get:

[0057]

[0058] In this way, the original problem is transformed into a numerical solution of the system of equations (13) about the unknowns λ and α. Here, the pre - conditioned conjugate gradient method is used to iteratively solve equation (13). First, a pre - conditioning matrix P is introduced, and both sides of the first equation of equation (13) are left - multiplied by the matrix P to eliminate α. Thus, equation (14) is obtained:

[0059]

[0060] where,

[0061] P=I - G(G T G) -1 G T (15)

[0062] By using iterative solution for equation (14), λ is obtained, and then it is substituted back into equation (11) to solve for α. Substituting λ and α into equation (12), the displacement vector at the current time step of each solution domain is obtained, and the calculation result u n+1 of the current time step is saved.

[0063] Please refer to Figure 3 , Figure 3 which gives the proportion of the single - iteration solution time. The calculation time of the program is mainly determined by the number of iterations and the single - iteration time of the program. During the iterative solution process, the matrix - solving calculation time accounts for 95% of the single - iteration time. Each iteration will perform the following steps: projection, calculation of beta, matrix - vector multiplication, calculation of alpha, update of the result vector, update of the residual vector, and judgment of the iteration number and iteration accuracy.

[0064] The pseudo - code of the algorithm for solving fluid - induced vibration based on the FETI method is as follows:

[0065]

[0066] Lines 8 - 10 in the pseudocode are the projection steps of the iterative solution process of the preconditioned conjugate gradient method, and there are a large number of matrix-vector multiplication operations in the projection. Calculate a in each subdomain k , λ k and r k , and it is necessary to calculate the matrix-vector operation F pk first. The matrix-vector multiplications are independent of each other and can be processed in parallel. Steps 8 - 10 and 17 in the pseudocode can be transplanted to domestic accelerators to improve the solution efficiency. For the Dawning advanced computing platform architecture (the specific experimental environment is shown in Table 4), heterogeneous transplantation is performed on the matrix-vector multiplication part, and the HIP parallel programming model is used to have the domestic accelerator and the CPU jointly complete the calculation of dense matrix-vector multiplication with a high computational density. The OpenMP parallel programming model is used for multi-threaded parallel operations. Figure 4 The process of solving the matrix-vector multiplication by multi-threaded parallel operations is given. One thread is used to control the domestic accelerator and call the rocBLAS linear algebra library for dense matrix calculation. The hipMalloc function is called by the thread to apply for video memory space, and the data stored on the CPU is transferred to the video memory of the domestic accelerator through hipMemcpy. The kernel function is started by calling the hipLaunchKernel function, and then the hipBLAS linear algebra library is called for dense matrix calculation. After the calculation on the domestic accelerator is completed, the result is transferred to the CPU by hipMemcpy, waiting for all devices to complete the calculation, synchronizing the data, and then entering the next iteration. Other threads control the CPU cores and calculate the remaining dense matrix-vector through calling the Intel MKL math library. During the iteration process, in order to avoid the time overhead caused by frequent data interaction between the CPU and the domestic accelerator, all dense matrices are transferred to the domestic accelerator before the solution starts. During the solution process, the CPU and the domestic accelerator execute corresponding tasks according to their respective task arrangements.

[0067] Table 4 Experimental Environment

[0068]

[0069] When performing parallel computing on multiple computing nodes, domestic accelerators need to read data from memory into video memory, and after the computing is completed, the computing results need to be written back to memory. As the problem scale increases, the time overhead of such memory access operations will gradually increase. To make full use of the computing performance of domestic accelerators, the present invention utilizes the double-buffering technology between memory and video memory to reduce the time of memory access. In addition, to make reasonable use of the video memory of domestic accelerators, a data alignment method is adopted to reduce the unnecessary occupation of video memory. Data is created at the alignment boundary to ensure that the cache line size is aligned with the vector data, avoiding additional calculations by the compiler, which may lead to performance loss. And the memory coalescing access technology is utilized to obtain higher throughput to improve the computing performance.

[0070] Step 5: Substitute the displacement value of the next time step obtained by iterative solution into equation (7) to update H n , and thus enter the solution of the next time step. The calculation accuracy of the NewMark method depends on the size of the time step Δt.

[0071] Repeat the above steps 2 to 5 until the set number of time steps is reached (note: the domain boundary balancing Figure 2 sub-algorithm is divided in the first time step and is not divided in subsequent steps). As Figure 5 shown, the present invention is tested on the Dawning advanced computing platform architecture. The computing volume is expanded from 48 million to 760 million, the number of computing nodes is expanded from 128 to 2048, each node processes 64 sub-domains, and the parallel efficiency reaches 54%. It can be clearly found that the present invention has good scalability, and the solution time is all at the minute level. The specific experimental data is shown in Table 5.

[0072] Table 5 Weak scalability experiment

[0073]

Claims

1. A high-precision fluid-induced vibration simulation method based on FETI, characterized in that it includes the following specific steps: Step 1. Read in the corresponding example data, including mesh data and material parameters; Step 2. Numerically discretize the fluid-induced vibration dynamic process using the NewMark method; Step 3. Use the FETI method to perform parallel decomposition on the discretized equations, divide the subdomains, and the subdomains are glued together by the Lagrange multiplier λ at the division boundary and satisfy the boundary conditions; According to the characteristics of the FETI subdomain boundary, use the domain boundary balanced graph bipartition algorithm to decompose the solution domain into I independent non-overlapping subdomains; Step 4. Use the preconditioned conjugate gradient method for iterative solution to obtain the Lagrange multiplier λ through iterative solution; Obtain the displacement according to the Lagrange multiplier λ and record the displacement value; Step 5. According to the displacement value, the NewMark method realizes the time-step update of the fluid-induced vibration process.

2. The high-precision fluid-induced vibration simulation method based on FETI according to claim 1, characterized in that: The mesh described in Step 1 includes tetrahedral meshes and hexahedral meshes; The material parameters include material density, elastic modulus, Poisson's ratio, time step, and the number of divided subdomains.

3. The high-precision fluid-induced vibration simulation method based on FETI according to claim 1, characterized in that: The graph bipartition algorithm described in Step 3 is based on the traditional greedy graph growth algorithm, and halo weights are added to the vertices.

4. The high-precision fluid-induced vibration simulation method based on FETI according to claim 1, characterized in that: For each iteration in Step 4, the following steps are executed: projection, calculation of beta, matrix-vector multiplication, calculation of alpha, update of the result vector, update of the residual vector, and judgment of the iteration number and iteration accuracy.

5. The high-precision fluid-induced vibration simulation method based on FETI according to claim 1, characterized in that: In Step 5, the calculation accuracy of the NewMark method depends on the size of the time step Δt.

Images