A method for design optimization of core melt collector based on least square particle method

Abstract

The application discloses a reactor core meltdown collector design optimization method based on a least square particle method, which is based on a least square moving particle method in a Lagrange method, has higher simulation accuracy and no space consistency limitation, adds a fluid-solid coupling model (PMS model) for calculating the interaction between a coolant and solid melt, and finally calculates the force between solid particles of the melt in the flow field and between the solid particles of the melt and a solid wall through a solid collision model (DEM method). In addition, an optimized particle slip model is used to relieve the particle particle aggregation effect, the calculation density and dynamic viscosity are processed to improve the calculation stability, and the accuracy of particle type determination is improved through a double detection method of algebraic calculation and geometric calculation. The method is a method with the ability of simulating the fluid-solid coupling problem and the solid collision problem in the reactor core meltdown collector and is used for the design optimization of the reactor core meltdown collector.

Description

Technical Field

[0001] This invention belongs to the field of nuclear reactor safety facility design technology, specifically relating to a core melt-collector design optimization method based on the least squares particle method. Background Technology

[0002] In the event of a severe accident in a sodium-cooled fast reactor (NSCRT), core disintegration is of particular concern due to its high fuel enrichment and coolant characteristics. In such cases, retaining molten material within the reactor core is a crucial mitigation measure to ensure pressure vessel integrity. However, the distribution and thickness of the molten core in the lower chamber significantly impact residual heat removal capacity. Therefore, designing a core melt collector with a rational structure and layout is an effective mitigation method. Consequently, the design of the core melt collector is of paramount importance for severe accident management in NSCRTs.

[0003] The interaction between the core melt collector and the molten core involves strong fluid-structure interaction (FSI) and solid-solid collisions. Traditional Eulerian methods, when using meshes to handle FSI problems, face challenges such as phase interface reconstruction and free interfaces. Furthermore, the collector design requires clear observation of the molten core's accumulation shape on the collector, making the Eulerian method less than optimal. In contrast, Lagrangian methods using particle-based computational fluid dynamics (PVMs) demonstrate unique advantages in addressing these issues. The MPS method, in particular, has developed a series of related models, including the PMS model, surface tension model, and viscous force model, achieving considerable success in solving FSI and solid-solid collision problems. However, traditional MPS methods cannot overcome the drawback of excessive errors in discrete operators under asymmetric particle distributions, and these models cannot simultaneously simulate both solid-liquid coupling and solid-solid collisions in the core melt collector effectively. Summary of the Invention

[0004] The problem this invention aims to solve is to address the shortcomings of existing methods by providing a design optimization method for melt-down reactor cores based on the least squares particle method. This method is capable of simulating fluid-structure interaction and solid collision problems in melt-down reactor cores for core design optimization. This method is based on the least squares moving particle method (LSMPS), a Lagrangian method with higher simulation accuracy and no spatial consistency constraints. A fluid-structure interaction model (PMS model) and a solid collision model (DEM method) are added to form the least squares moving particle-discrete element method (LSMPS-DEM), which simulates the distribution and thickness of molten material in the melt-down reactor core of a sodium-cooled fast reactor.

[0005] The present invention adopts the following technical solution:

[0006] A core melting collector design optimization method based on the least squares particle method includes the following steps:

[0007] S1. Based on the simulated operating conditions of the core melt collector, set the initial position, velocity, pressure, phase, particle boundary type, and solid particle number information of different particles, and preprocess the setting information.

[0008] S2. Every ten time steps, call the Link-list algorithm to retrieve and store the particle number of each particle after preprocessing in step S1 and the particles within its cutoff radius.

[0009] S3. A dual detection method combining geometric and algebraic calculations is used to determine the particle surface type;

[0010] S4. Calculate the small correction displacement necessary for the least squares moving particle method based on the particle surface type determined in step S3, ensuring calculation accuracy while avoiding particle displacement distortion.

[0011] S5. Enter the smoothing process to smooth the density and dynamic viscosity coefficient of all particles, reducing the possibility of instability in pressure calculation.

[0012] S6. Explicitly calculate the gravity, viscous force, and surface tension of all particles, including solid particles; ignore the influence of pressure and use gravity, viscous force, and surface tension to solve for the temporary velocity and temporary position of the particles.

[0013] S7. Use the high-precision discretization operator of the least squares moving particle method to discretize the terms in the pressure Poisson equation, and then use the stable double conjugate gradient method to implicitly solve the pressure Poisson equation to obtain the pressure value at the next moment. Update the position and velocity of the particle according to the pressure value.

[0014] S8. Enter the fluid-structure interaction subroutine to correct the solid particles in the molten core, calculate the center of gravity, velocity, and angular velocity of the solid particles in the molten core; and correct the velocity and position of the solid particles contained in the molten core solid particles based on the calculated center of gravity, velocity, and angular velocity of the solid particles in the molten core solid particles.

[0015] S9. Enter the solid collision subroutine loop. Based on the center of gravity, velocity and angular velocity of the molten solid particles calculated in step S8, calculate the interaction force between the molten solid particles within the solid collision time step. Update the position and velocity of the molten solid particles accordingly. The time step of the solid collision calculation program should be much smaller than the least squares discrete time step. After the same loop time, exit the solid collision calculation program and enter the next least squares discrete time step.

[0016] Specifically, in step S2, the Link-list algorithm is used to improve particle retrieval efficiency. The specific operation is as follows:

[0017] Within each time step, the region containing all particles, including those at the wall boundaries, is divided into grids. Then, the number of particles in each grid and their corresponding particle numbers are counted. If the particles are solid, the coordinates of their centroids are counted. This process yields the grid containing the centroid of each particle or solid particle and its coordinates. Since the search range for a particle is a finite-sized circular domain, each particle only needs to search its own grid and the surrounding grids.

[0018] Specifically, in step S3, a dual detection method combining geometric and algebraic calculations is used to determine the particle surface type to improve efficiency and accuracy, as follows:

[0019] The particle surface types are categorized into internal particles, internal particles near the surface, surface particles, and isolated particles. First, through geometric calculations, N0 quadrants are generated with the discriminant particle as the origin, and the projection angle θ of all neighboring particles of that discriminant particle is calculated. j and its quadrant N j The calculation formula is as follows:

[0020]

[0021]

[0022] Where, r β Let r be the projection of the interparticle distance in the y-direction. α Let θ0 be the projection of the particle spacing onto the x-direction, and let θ0 be the size of the angle contained in each quadrant.

[0023] Let Γ(i) be the number of quadrants in which particle i has neighboring particles. Then the proportion of quadrants in which no neighboring particles exist is:

[0024]

[0025] When Ψ(i) is greater than 0.2, particle i is initially identified as a surface particle; otherwise, it is an internal particle.

[0026] After the geometric calculation steps are completed, algebraic calculations are performed; first, the number of neighboring particles N of the particle is calculated. sum and particle number density n′, i.e.

[0027]

[0028] Where w ij For kernel functions; when N sum When n < 6 or n′ < 0.6n0, the particle is determined to be an isolated particle; where n0 represents the initial particle number density. After determining the surface particles and isolated particles, the particles inside the surface are then determined, using the following criteria:

[0029]

[0030] ifsurf(i) is the particle type determination value for particle i. Its values ​​0, 1, 2, and 3 represent that the particle is an internal particle, an internal particle near the surface, a surface particle, and an isolated particle, respectively.

[0031] Specifically, in step S4, an optimized particle slip model is used to calculate the small correction displacement required for the least squares moving particle method:

[0032]

[0033] Where, δr i C represents the corrected displacement of particle i. shift Values ​​less than 0.5, h is the scan radius, C i Let i be the particle concentration at particle i; This represents the unit normal vector of the concentration gradient direction at particle i.

[0034] Specifically, in step S5, the particle density and dynamic viscosity coefficient of particle i after smoothing are obtained by weighted averaging of the density and dynamic viscosity coefficient of the particles in the neighborhood of particle i. The calculation formula is as follows:

[0035]

[0036]

[0037] in and These represent the particle density and dynamic viscosity coefficient of particle i after smoothing, respectively, where j represents the neighboring particle, ρ j Let j be the density of particle j. The kernel function used in the smoothing subroutine has the following computational model:

[0038]

[0039] Among them, among them, Let α be the distance between particles i and j, and let α be the spatial dimension coefficient. d The one-dimensional, two-dimensional, and three-dimensional conditions have different values, namely 2 / r. e ,60 / (7πr e 2 ) and 12 / (πr e 3 Since the calculation cases are all two-dimensional, therefore α d =60 / (7πr) e 2 ), where r eTo smooth the cutoff radius, considering both computational stability and workload, we choose r. e = 4.1Δl, where Δl is the initial distance between adjacent particles.

[0040] Specifically, in step S6, the gravity, viscous force, and surface tension of all particles, including solid particles, are explicitly calculated; the influence of pressure is ignored, and the temporary velocity and temporary position of the particles are solved using gravity, viscous force, and surface tension.

[0041]

[0042]

[0043] The first item within the parentheses on the right is viscous force. This represents the vector sum of surface tension and gravity, where v is the kinematic viscosity. Let be the temporary velocity of particle i. Let be the temporary position of particle i, and Δt be the calculation time step. Let be the velocity of particle i at time n. Let be the position of particle i at time n.

[0044] Specifically, in step S7, the terms in the pressure Poisson equation are discretized using the high-precision discretization operator of the least squares moving particle method, and then the pressure Poisson equation is implicitly solved using the stable double conjugate gradient method to obtain the pressure value at the next moment. The position and velocity of the particle are then updated based on the pressure value.

[0045]

[0046]

[0047]

[0048] in, ρ represents the change in velocity of particle i caused by the pressure gradient, and ρ represents density. Let n and n represent the pressure, velocity, and position of particle i at time n+1, respectively.

[0049] Specifically, in step S8, the fluid-structure interaction subroutine is used to correct the solid particles of the molten core, calculating the center of gravity, velocity, and angular velocity of the solid particles:

[0050]

[0051]

[0052]

[0053] in, Let be the coordinates of the centroid of solid particle ii at time n. Let be the velocity and angular velocity of solid particle ii at time n+1, respectively. Let be the angular velocity of solid particle ii at time n; N is the number of solid particles contained in the solid particle; m i Let i be the mass of particle i; These represent the velocities of solid particle i at time n and time n+1, respectively. Before proceeding to this step, the solid particles in the solid particle are treated as liquid particles for calculation.

[0054] The velocity and position of particle i in the molten solid particles are corrected based on the centroid position, velocity, and angular velocity of the solid particles in the core. The calculation formula is as follows:

[0055]

[0056]

[0057] The coordinate transformation matrix M is in the form of:

[0058]

[0059] in, Let n be the magnitude of the angular velocity of the molten solid particles at time n.

[0060] Specifically, in step S9, based on the velocity of the molten solid particles calculated in step S8, the molten solid particles will move for one solid collision calculation time step to obtain their position at the next calculation moment; based on this position, the solid collision module will calculate the contact situation between each molten solid particle and the molten solid particles and collector wall in the search domain, and calculate the force of the solid particles colliding with each other based on the spring-damping model:

[0061]

[0062]

[0063]

[0064] in, δ represents the normal and tangential contact forces between molten solid particles i and j, respectively; n δ s These are the relative displacements of the contact points between the solid particles of the molten material in the normal and tangential directions, respectively, k. n k n These are the normal and tangential stiffnesses, respectively; c n c sThese are the normal and tangential damping coefficients, respectively; μ is the friction coefficient between solid particles; t represents time; the molten solid particles after overlapping will generate a pushing force, thus pushing each other apart, and the acceleration of this pushing is:

[0065]

[0066]

[0067] in These are the linear and angular accelerations of the molten solid particles ii at time k. and These are the force and torque acting on the center of mass of the molten solid particle ii at time k during the solid collision process, respectively. ii It is the mass of the solid particles in the molten material ii; I ii Let be the moment of inertia of the solid particles ii in the molten material;

[0068] Using the calculated acceleration and angular acceleration, update the position, angle, velocity, and angular velocity at the next solid collision calculation moment:

[0069]

[0070]

[0071]

[0072]

[0073] in, Let Δt represent the velocity, angular velocity, angle, and displacement of solid particle ii at time k+1, respectively. DEM Indicates the time step for solid collision calculation;

[0074] After determining the position and velocity at the next time step, repeat the above operation. Next, the positions of the final molten solid particles in the core melt collector are obtained, where Δt LSMPS This represents the time step in the least squares moving particle method discrete calculation process.

[0075] Compared with the prior art, the present invention has at least the following beneficial effects:

[0076] This invention is a design optimization method for core melt collectors based on the least squares particle method. By employing a high-precision discretization method using the least squares moving particle method, and adding a fluid-structure interaction model and a solid-state collision model, it calculates the interactions between coolant and molten material, molten material with molten material, and molten material with the collector wall in the collector. This enables it to simulate the distribution and accumulation thickness of molten material in core melt collectors with different structures and arrangements in sodium-cooled fast reactors.

[0077] The least squares moving particle method employs an arbitrary high-precision meshless discretization scheme, derived through Taylor expansions of arbitrary order. It does not assume spatial consistency and is superior to the traditional MPS discretization method in both accuracy and limitations.

[0078] The Link-list algorithm is used to search for each particle and all particles within its cutoff radius. Since the search range for a particle is a finite-sized circular domain, each particle only needs to be searched within its own grid and the surrounding grids, which greatly improves particle search efficiency.

[0079] By employing a dual detection method that combines geometric and algebraic calculations, the particle surface type can be determined more accurately.

[0080] An optimized particle slip model based on Fick's law is adopted to slightly modify the particle displacement at each time step, so as to maintain uniform particle spacing and reduce numerical instability caused by particles being too close together.

[0081] Smoothing is applied to the density and dynamic viscosity coefficients. In solid-liquid two-phase calculations, there are often density or viscosity ratios of tens or hundreds of times on both sides of the phase interface. Such large numerical differences can lead to non-convergence in solving the pressure Poisson equation. Smoothing is performed on the particle density and dynamic viscosity coefficient using a kernel function, which can effectively alleviate the computational instability caused by abrupt changes in physical properties.

[0082] In summary, this method utilizes the latest least-squares moving particle method to discretize the gradient and Laplace models, improving computational accuracy. The application of a linked list retrieval method significantly enhances retrieval efficiency. Simultaneously, an optimized system slip model is employed to mitigate particle aggregation effects, resulting in a more uniform particle distribution. Furthermore, solid-structure collision and fluid-structure interaction (FSI) models are added, enabling the method to simulate the FSI and solid-structure collision processes during the migration of molten material from the core of a sodium-cooled fast reactor to the collector. Attached Figure Description

[0083] Figure 1 This is a flowchart of the calculation process.

[0084] Figure 2This is a schematic diagram of the Link-list algorithm for retrieval.

[0085] Figure 3 This is a schematic diagram of a spring-damped model.

[0086] Figure 4 A schematic diagram of a two-dimensional single chimney operating condition.

[0087] Figure 5a The diagram shows the simulation results for a chimney tilt angle of 15°.

[0088] Figure 5b The diagram shows the simulation results for a chimney tilt angle of 30°.

[0089] Figure 5c The diagram shows the simulation results when the chimney tilt angle is 45°.

[0090] Figure 6a The simulation results are shown in the figure where the working medium of the molten material is zirconium oxide particles.

[0091] Figure 6b The simulation results are shown in the figure where the molten working medium is stainless steel particles.

[0092] Figure 6c The simulation results are shown in the figure where the working medium of the molten material is alumina particles.

[0093] Figure 7a The image shows the simulation results for particles with a diameter of 1 mm.

[0094] Figure 7b The image shows the simulation results for particles with a diameter of 2 mm. Detailed Implementation

[0095] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0096] This invention discloses a core melt collector design optimization method based on the least squares particle method. Please refer to [link to relevant documentation]. Figure 1 The calculation flowchart and specific steps are as follows:

[0097] S1. Based on the simulated operating conditions of the core melt collector, set the initial position, velocity, pressure, phase, particle boundary type, and solid particle number information for three types of particles: liquid water, core melt collector, and core melt. Then, preprocess the set information.

[0098] S2. Every ten time steps, call the Link-list algorithm to retrieve and store the particle number of each particle after preprocessing in step S1 and the particles within its cutoff radius.

[0099] Please see Figure 2The Link-list retrieval diagram shown illustrates how, within each time step, the region containing all particles, including those at the wall boundaries, is divided into grids. Then, the number of particles within each grid and their corresponding particle numbers are counted. For solid particles, the centroid coordinates are calculated. This yields the grid containing the centroid of each particle or solid particle, along with its coordinates. Since the retrieval range of a particle is a finite-sized circular domain, each particle only needs to search its own grid and the surrounding grids. In a two-dimensional scenario, nine grids are searched. The grid size scanned in this invention is 6.1Δl, and the retrieval radius of neighboring particles is also 6.1Δl, which is 1.5 times the cutoff radius of the kernel function. Here, Δl represents the initial spacing between adjacent particles.

[0100] S3. Identify and determine the surface type of particles:

[0101] A dual detection method based on geometric and algebraic computation is employed to classify particle surface types into internal particles, near-surface internal particles, surface particles, and isolated particles. First, through geometric computation, N0 quadrants are generated with the discriminant particle as the origin, and the projection angle θ of all neighboring particles of that discriminant particle is calculated. j and its quadrant N j The calculation formula is as follows:

[0102]

[0103]

[0104] Where, r β Let r be the projection of the interparticle distance in the y-direction. α Let θ0 be the projection of the particle spacing onto the x-direction, and let θ0 be the size of the angle contained in each quadrant.

[0105] Let Γ(i) be the number of quadrants in which particle i has neighboring particles. Then the proportion of quadrants in which no neighboring particles exist is:

[0106]

[0107] When Ψ(i) is greater than 0.2, particle i is initially identified as a surface particle; otherwise, it is an internal particle.

[0108] After the geometric calculation steps are completed, algebraic calculations are performed; first, the number of neighboring particles N of the particle is calculated. sum and particle number density n′, i.e.

[0109]

[0110] Where w ij For kernel functions; when N sumWhen n < 6 or n′ < 0.6n0, the particle is determined to be an isolated particle; where n0 represents the initial particle number density. After determining the surface particles and isolated particles, the particles inside the surface are then determined, using the following criteria:

[0111]

[0112] ifsurf(i) is the particle type determination value for particle i. Its values ​​0, 1, 2, and 3 represent that the particle is an internal particle, an internal particle near the surface, a surface particle, and an isolated particle, respectively.

[0113] S4. Enter the particle displacement correction subroutine. Based on the particle surface type determined in step S3 and the optimized particle slip model, calculate the small correction displacement of the least squares moving particle method to eliminate or weaken the velocity and displacement of the outward expansion of surface particles. This ensures high accuracy while avoiding distortion of particle displacement calculation using the least squares moving particle method. The specific calculation model is as follows:

[0114]

[0115] Where, δr i C represents the corrected displacement of particle i. shift Values ​​less than 0.5, h is the scan radius, C i Let i be the particle concentration at particle i; This represents the unit normal vector of the concentration gradient direction at particle i.

[0116] S5. Enter the smoothing process subroutine. Apply a kernel function to smooth the density and dynamic viscosity coefficient of all calculated particles. For calculated particle i, its density and dynamic viscosity coefficient are weighted averaged based on its neighboring particles. The calculation formula is as follows:

[0117]

[0118]

[0119] in and These represent the particle density and dynamic viscosity coefficient of particle i after smoothing, respectively, where j represents the neighboring particle, ρ j Let j be the density of particle j. The kernel function used in the smoothing subroutine has the following computational model:

[0120]

[0121] in, Let α be the distance between particles i and j, and let α be the spatial dimension coefficient. d The one-dimensional, two-dimensional, and three-dimensional conditions have different values, namely 2 / r.e ,60 / (7πr e 2 ) and 12 / (πr e 3 Since the calculation cases are all two-dimensional, therefore α d =60 / (7πr) e 2 ), where r e To smooth the cutoff radius, considering both computational stability and workload, we choose r. e = 4.1Δl, where Δl is the initial distance between adjacent particles;

[0122] Substituting the smoothed values ​​of particle density and dynamic viscosity coefficient into the momentum equation, we obtain the following expression:

[0123]

[0124] in, Let t be velocity, t be time, and p be pressure. Represents gravity. Indicates surface tension;

[0125] S6. Explicitly calculate the gravity, viscous force, and surface tension of all particles, including the solid particles that make up the solid particles of the melt. Then, based on the gravity, viscous force, and surface tension, explicitly solve for the temporary position and velocity of all particles.

[0126] Furthermore, ignoring the influence of pressure, the effects of gravity, surface tension, and viscous force in the momentum equation are explicitly calculated to obtain the temporary velocities and positions of all particles; the calculation equations are as follows:

[0127]

[0128]

[0129] The first item within the parentheses on the right is viscous force. Let v be the vector sum of surface tension and gravity, and v be the kinematic viscosity. It is the temporary velocity of particle i. Δt is the temporary position of particle i, and Δt is the calculation time step. Let be the velocity of particle i at time n. It is the position of particle i at time n.

[0130] In step S7, the high-precision discretization operator of the least squares moving particle method is used to discretize the terms in the pressure Poisson equation, and then the pressure value at the next moment is obtained by implicitly solving the pressure Poisson equation using the stable double conjugate gradient method. The position and velocity of the particle are then updated based on the pressure value.

[0131]

[0132]

[0133]

[0134] in, ρ represents the change in velocity of particle i caused by the pressure gradient, and ρ represents density. Let represent the pressure, velocity, and position of particle i at time n+1, respectively;

[0135] S8. Enter the fluid-structure interaction subroutine to correct the solid particles in the core molten material; calculate the position and velocity of the center of gravity and angular velocity of the solid particles in the core molten material, and then update the position and velocity of the solid particles in the molten material.

[0136] Since step S7 treats the solid particles that make up the molten core as another component of the fluid and calculates the physical quantities of individual particles, this step restores the molten core solid particles according to the fluid-structure interaction model and calculates the center of gravity, velocity, and angular velocity of the molten core solid particles:

[0137]

[0138]

[0139]

[0140] in, Let be the coordinates of the centroid of solid particle ii at time n. Let be the velocity and angular velocity of solid particle ii at time n+1, respectively. Let be the angular velocity of solid particle ii at time n; N is the number of solid particles contained in the solid particle; m i Let i be the mass of particle i; Let represent the velocities of solid particle i at time n and time n+1, respectively.

[0141] Furthermore, the position and velocity of particle i in the molten solid particles are corrected, and the calculation formula is as follows:

[0142]

[0143]

[0144] The coordinate transformation matrix M is in the form of:

[0145]

[0146] in, Let n be the magnitude of the angular velocity of the solid particles in the molten material at time n;

[0147] S9. Based on the velocity of the molten solid particles calculated in step S8, the molten solid particles will move for one solid collision calculation time step to obtain the position at the next calculation moment. Based on this position, the solid collision module will calculate the contact situation between each molten solid particle and the molten solid particles and collector wall in the search domain, and calculate the collision force between the molten solid particles according to the spring-damping model. (See spring-damping model for details.) Figure 3 Calculate the collision forces of solid particles:

[0148]

[0149]

[0150]

[0151] in, δ represents the normal and tangential contact forces between molten solid particles i and j, respectively; n δ s These are the relative displacements of the contact points between the solid particles of the molten material in the normal and tangential directions, respectively; c n c s These are the normal and tangential damping coefficients, respectively; μ is the friction coefficient between solid particles; t represents time; k n k n These are the normal and tangential stiffnesses, respectively, calculated using the following formula:

[0152]

[0153] in, This represents the equivalent mass of solid particles i and j in the molten material; The dimensionless correction coefficient; Δt represents the time step, taken as... Tangential stiffness coefficient k s and viscosity coefficient c s The calculation formula is as follows:

[0154]

[0155] Where v represents Poisson's ratio;

[0156] Furthermore, the overlapping solid particles of the molten core will generate thrust, pushing them apart. The acceleration of this pushing action is:

[0157]

[0158]

[0159] in These are the linear and angular accelerations of the molten solid particles ii at time k. and These are the force and torque acting on the center of mass of the molten solid particle ii at time k during the solid collision process, respectively. ii It is the mass of the solid particles in the molten material ii; I ii Let be the moment of inertia of the solid particles ii in the molten material;

[0160] Furthermore, using the calculated acceleration and angular acceleration, the position, angle, velocity, and angular velocity of the core molten material particles at the next solid collision calculation moment are updated:

[0161]

[0162]

[0163]

[0164]

[0165] in, Let Δt represent the velocity, angular velocity, angle, and displacement of solid particle ii at time k+1, respectively. DEM Indicates the time step for solid collision calculation;

[0166] After determining the position and velocity of the molten solid particles at the next time step, repeat step S8. Next, the positions of the molten solid particles in the core melt collector are obtained after a least-squares discrete calculation step, where Δt LSMPS This represents the time step for discrete calculations using the least squares moving particle method.

[0167] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be further described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0168] This invention presents a design optimization method for a reactor core meltdown collector based on the least squares particle method. To verify the feasibility of this technical solution, this embodiment employs a two-dimensional modeling of the reactor core meltdown collector under single-chimney operation (see [link]). Figure 4 The reactor core melt and collector are solid, while the cooling water is liquid. First, simulation experiment 1 was conducted, simulating chimney inclination angles of 15°, 30°, and 45° with the melt material being zirconium oxide and a particle size of 1 mm. The simulation results are as follows: Figure 5a , Figure 5b , Figure 5c As shown in Figure 2, three materials—zirconia, stainless steel, and alumina—were selected to simulate solid particles from the reactor core melt. The simulation results for these three materials at a chimney inclination angle of 15° and a particle size of 1 mm are shown below. Figure 6a , Figure 6b , Figure 6c As shown; finally, simulations were performed on alumina particles with diameters of 1 mm and 2 mm at a chimney inclination angle of 15°, and the simulation results are shown below. Figure 7a and 7b As shown in the table below, the stiffness, viscosity, and other coefficients used for different particle sizes are as follows:

[0169]

[0170] After obtaining the distribution of the melt in the collector under three different operating conditions, the average height Z of the particles relative to the bottom of the collector can be analyzed. avg And variance Z var The effect of particle distribution is evaluated, and the core melt collector is then improved and optimized.

[0171] The above content is only for illustrating the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. Any modifications made to the technical solution based on the technical concept proposed in this invention shall fall within the scope of protection of the claims of this invention.

Claims

1. A design optimization method for a core melting collector based on the least squares particle method, characterized in that, Includes the following steps: S1. Based on the simulated operating conditions of the core melt collector, set the initial position, velocity, pressure, phase, particle boundary type, and solid particle number information of different particles, and preprocess the setting information. S2. Every ten time steps, call the Link-list algorithm to retrieve and store the particle number of each particle after preprocessing in step S1 and the particles within its cutoff radius. S3. A dual detection method combining geometric and algebraic calculations is used to determine the particle surface type; S4. Calculate the small correction displacement necessary for the least squares moving particle method based on the particle surface type determined in step S3, ensuring calculation accuracy while avoiding particle displacement distortion. S5. Enter the smoothing process to smooth the density and dynamic viscosity coefficient of all particles, reducing the possibility of instability in pressure calculation. S6. Explicitly calculate the gravity, viscous force, and surface tension of all particles, including solid particles; neglecting the effect of pressure, use gravity, viscous force, and surface tension to solve for the temporary velocity and temporary position of the particles: The first item within the parentheses on the right is viscous force. This represents the vector sum of surface tension and gravity. Kinematic viscosity, For particles Temporary speed, For particles temporary location To calculate the time step, For the particle at time n speed, For the particle at time n Location; S7. Discretize the terms in the pressure Poisson equation using a high-precision discretization operator based on the least squares moving particle method, then implicitly solve the pressure Poisson equation using the stable double conjugate gradient method to obtain the pressure value at the next moment, and update the particle's position and velocity based on the pressure value: in, This indicates particles caused by a pressure gradient. The change in velocity, Indicates density, , , They represent the particles at time n+1, respectively. Pressure, speed, and position; S8. Enter the fluid-structure interaction subroutine to correct for the solid particles in the molten core, and calculate the center of gravity, velocity, and angular velocity of the solid particles in the molten core: in, solid particles The centroid coordinates at time n, , Solid particles at time n+1 Velocity, angular velocity; For solid particles at time n angular velocity; This refers to the number of solid particles contained within a solid particle. For particles The quality; , These represent the velocities of solid particle i at time n and time n+1, respectively. Before proceeding to this step, the solid particles in the solid particle are treated as liquid particles for calculation. Based on the center of gravity, velocity, and angular velocity of the solid particles in the molten core, the particles in the molten core are analyzed. The speed and position are corrected using the following formula: Wherein, coordinate transformation matrix The form is: in, Let n be the magnitude of the angular velocity of the solid particles in the molten material at time n; S9. Based on the velocity of the molten solid particles calculated in step S8, the molten solid particles will move for one solid collision calculation time step to obtain the position at the next calculation moment; based on this position, the solid collision module will calculate the contact situation between each molten solid particle and the molten solid particles and collector wall in the search domain, and calculate the force of the collision between solid particles according to the spring-damping model: in, , These are the normal and tangential contact forces between molten solid particles i and j, respectively. , These represent the relative displacements in the normal and tangential directions at the contact points between the solid particles of the molten material. , These are the normal and tangential stiffness, respectively. , These are the normal and tangential damping coefficients, respectively; The coefficient of friction between solid particles; Indicates time; the overlapping molten solid particles will generate a pushing force, thus pushing each other apart. The acceleration of this pushing is: in , These are the linear and angular accelerations of the molten solid particles ii at time k. and These are the force and torque acting on the center of mass of the molten solid particle ii at time k during the solid collision process, respectively. It is the mass of the solid particles in the molten material (ii); Let be the moment of inertia of the solid particles ii in the molten material; Using the calculated acceleration and angular acceleration, update the position, angle, velocity, and angular velocity at the next solid collision calculation moment: in, , , , Representing solid particles Velocity, angular velocity, angle, and displacement at time k+1; Indicates the time step for solid collision calculation; After determining the position and velocity at the next time step, repeat step S9. Next, the positions of the molten solid particles in the core melt collector are obtained after a least-squares discrete calculation step size, where... This represents the time step for discrete calculations using the least squares moving particle method.

2. The core melting collector design optimization method based on the least squares particle method according to claim 1, characterized in that, In step S2, the Link-list algorithm is used to improve particle retrieval efficiency, as detailed below: Within each time step, the region containing all particles, including those at the wall boundaries, is divided into grids. Then, the number of particles in each grid and their corresponding particle numbers are counted. If the particles are solid, the coordinates of their centroids are counted. This process yields the grid containing the centroid of each particle or solid particle and its coordinates. Since the search range for a particle is a finite-sized circular domain, each particle only needs to search its own grid and the surrounding grids.

3. The core melting collector design optimization method based on the least squares particle method according to claim 1, characterized in that, In step S3, a dual detection method combining geometric and algebraic calculations is used to determine the particle surface type to improve efficiency and accuracy, as detailed below: Particle surface types are categorized into internal particles, internal particles near the surface, surface particles, and isolated particles. First, geometric calculations are performed, with the discriminant particle as the origin for generation. In each quadrant, calculate the projection angles of all neighboring particles of the discriminant particle. and its quadrant The calculation formula is as follows: in, For the interparticle distance at Projection in direction, For the interparticle distance at Projection in direction, The size of the angles contained in each quadrant; Let the particle The number of quadrants in which neighboring particles exist is The proportion of quadrants without neighboring particle distributions is when When it is greater than 0.2, the particles They are initially identified as surface particles; otherwise, they are internal particles. After the geometric calculation steps are completed, algebraic calculations are performed; first, the number of neighboring particles of the particle is calculated. and particle number density ,Right now in For kernel function; when or At that time, the particle was determined to be an isolated particle; among them, This represents the initial particle number density; after judging surface particles and isolated particles, the particles inside the surface are judged, and the criteria used are: For particles The particle type determination value, with values ​​of 0, 1, 2, and 3 representing that the particle is an internal particle, an internal particle near the surface, a surface particle, and an isolated particle, respectively.

4. The core melting collector design optimization method based on the least squares particle method according to claim 1, characterized in that, In step S4, an optimized particle slip model is used to calculate the small correction displacement required for the least squares moving particle method: in, Represents particles Corrected displacement, Values ​​less than 0.5 It is the scan radius. For particles Particle concentration at the location; This represents the unit normal vector of the concentration gradient direction at particle i.

5. The core melting collector design optimization method based on the least squares particle method according to claim 1, characterized in that, In step S5, by calculating the particles The particle density and dynamic viscosity coefficient of the neighboring particles are weighted and averaged to obtain the particle... The particle density and dynamic viscosity coefficient after smoothing are calculated using the following formulas: in and Particles Particle density and dynamic viscosity coefficient after smoothing treatment Represents a neighboring particle. For particles density, The kernel function used in the smoothing subroutine has the following computational model: in, For particles i and particles j Particle spacing, spatial dimension coefficient The conditions in one dimension, two dimensions, and three dimensions have different values, respectively. , and Since the calculation conditions are all two-dimensional, therefore ,in, To smooth the cutoff radius, considering both computational stability and workload, we take... ,in This represents the initial spacing between adjacent particles.

Images