An S-type method for shielding surface connection N -MC coupling calculation method

Abstract

An S-type method for shielding surface connection N MC coupling calculation method, in the reactor shielding calculation, the shape of some complex geometric areas is spherical or ellipsoidal, generally S N The computational grid of the method cannot accurately describe this type of coupling surface, and an approximate solution needs to be obtained by using intersecting grid mapping. Taking the sphere as an example, in S N After completing the transport calculation for a slightly larger area than the surface source area, according to S N The z-direction grid is divided into multiple layers of circular rings, and the projection of the rings and S N The intersection of the grid surfaces is determined using the coordinates of the intersection point and the height of the annulus to determine the normal vector at that location on the surface. The angular fluence rate at the center of the SN grid where the intersection point is located is considered the outgoing angular fluence rate for the annulus mapping area. The angular fluence rate can be used to determine the outgoing flux based on particle type, surface source, energy group, surface layer, grid, and direction. This outgoing flux is then converted into a cumulative distribution function, which serves as the source term input to the MC program for continued calculations beyond the surface.

Description

Technical Field

[0001] The present invention relates to the field of nuclear reactor shielding calculation, and is a discrete ordinate-Monte Carlo coupling calculation method for curved surface geometry (including spherical and ellipsoidal surfaces). Background Art

[0002] In order to combine the discrete vertical scale method (hereinafter referred to as S N The advantages of the Monte Carlo method (hereinafter referred to as the MC method) are fast solution speed and accurate geometric modeling. In the 1970s, Emmet et al. of the Oak Ridge National Laboratory in the United States proposed the two-dimensional S N -MC coupling method, based on two-dimensional S N The intermediate interface DOMINO was developed between the differential program DOT and the three-dimensional MC program Morse. The core idea of ​​this method is to use S N The angular fluence rate provided by the calculation is converted into the cumulative distribution function of energy, space and angle as the source term input of the MC calculation. In the subsequent development of this method, three-dimensional S N The coupled calculation of plane and cylindrical surfaces in xyz, rθz coordinate systems by the program and MC program. By 2011, Kulesza of the University of Tennessee in the United States developed the DISCO software package, which can realize the coupled calculation of different geometries of two types of coordinate systems in two-dimensional and three-dimensional, indicating that the method has been developed more maturely in common geometries. Since the geometric structure of the reactor is usually dominated by cylinders, prisms, etc., the plane and cylindrical surface sources of common geometries can meet the coupling calculation requirements in most cases, but for the spherical or ellipsoidal parts of some special components, spherical and ellipsoidal surfaces need to be used as the surface sources for continued calculation. At the same time, taking into account S N The xyz, rθz coordinate system grid in this method cannot accurately describe the geometry of the two surfaces mentioned above. Therefore, it is necessary to develop a method to approximately solve the surface source in the xyz coordinate system. Summary of the Invention

[0003] In order to overcome the problems of the prior art, the present invention aims to provide an S-type shielding system for shielding calculation surface connection. N -MC coupling calculation method, taking the sphere as an example, the present invention processes the coupling calculation problem of the surface in S N After completing the transport calculation for a slightly larger area than the surface source area, according to S N The z-direction grid interval is calculated to divide the spherical surface source into multiple layers of circular rings, and the circular projection and S N Divide the intersection of the xy plane grid and determine the S through which the sphere passes according to the intersection position NThe angular fluence rate of the center points of these grids is approximately regarded as the angular fluence rate on the arc segment passing through the grid. At the same time, the cosine of the angle between the normal vector of the arc segment and the xyz coordinate system is determined according to the intersection point and the height in the z direction, so as to obtain the surface source normal vector emitted at different positions. m,q,g,s The particle outflow on particle type (p), surface source (s), energy group (g), annular layer (k), grid (q), and direction (m) can be obtained, and then the outflow is converted into a cumulative distribution function (CDF) as the source term probability distribution input of the MC program to perform source sampling and continued calculation of areas outside the surface source. Compared with using S N The coupling method directly calculates the entire area, and the coupling method uses MC modeling to more accurately describe the local surface area of ​​the reactor by defining the surface source, and the calculation result is more accurate; compared with the MC method that directly calculates the entire area, the coupling method divides the regular geometric area into S N By reducing the MC calculation area, the calculation time of the MC method is reduced and the calculation efficiency is improved.

[0004] In order to achieve the above objectives, the present invention adopts the following technical solutions to be implemented:

[0005] Step 1: Discrete ordinate method S N The transport calculations yield the grid angle flux rate ψ for the reactor geometry including the curved surface source region. m,q,g,s ;

[0006] Step 2: According to S N The z-direction layering divides the surface into annular zones, and calculates the annular zone projection and S N The intersection of the grid surface, according to the number of intersection points, the ring belt is divided into the same number of ring surface grids, and then the grid on the ring surface and S are obtained. N The mapping relationship of the grid; the specific mapping relationship is the torus grid number and S N Grid numbers correspond to:

[0007] Torus grid numbering format [Number, Layer], S N The grid number format is [i, j, k], where Layer corresponds to k one by one, and the number is obtained by mapping the formula:

[0008] Number=(j-jsuridx1)*(isuridx2-isuridx1+1)+i-isuridx1 (1)

[0009] Where:

[0010] Layer——the layer number corresponding to a certain ring belt;

[0011] Number——the torus grid number on a certain annulus;

[0012] i,j,k——S N The grid numbers for transport calculation division in the x, y, and z directions;

[0013] isuridx1,isuridx2——S N The lower and upper boundary grid numbers for transport calculation in the x-direction;

[0014] jsuridx1——S N The lower boundary grid number for transport calculation in the y direction;

[0015] Step 3: After obtaining the grid mapping relationship, the corresponding S N The angular flux rate at the center of the grid is regarded as the angular flux rate emitted from the annular grid; at the same time, since only the outflow along the positive normal direction of the surface source is considered, the normal vectors at different positions on the surface need to be calculated: according to the annular projection calculated in step 2 and S N The intersection of the mesh surface and the center coordinates of the surface get the normal vector of the intersection:

[0016] Normal vector of a point on the sphere: [2(x-r_x), 2(y-r_y), 2(z-r_z)],

[0017] Normal vector of a point on the ellipsoid:

[0018] in:

[0019] x,y,z——Surface and S N The coordinates of the grid intersections;

[0020] r_x, r_y, r_z——the coordinates of the center of the sphere and ellipsoid;

[0021] a, b, c - the axis lengths of the ellipsoid in three directions;

[0022] The normal vector and S N The calculated discrete direction m is subjected to vector dot product operation, and if the result is positive, it is determined to be forward emission;

[0023] Step 4: According to the angular fluence rate, the outflow on direction m, grid q, layer k, energy group g, surface s, and particle type p is obtained in sequence;

[0024] The process of calculating the jet from the surface is as follows:

[0025] (1) Obtain the particle flow of a certain layer by summing the outgoing flows of the grids on the same layer;

[0026] (2) For curved surfaces, the normal vector changes with the grid position and needs to be calculated based on the particle flight direction Ω mThe dot product with the normal vector at different positions determines whether it is forward emission;

[0027] (3) When calculating the particle flow emitted from a certain grid, the grid area should be the grid area on the curved surface. The derived formula is as follows:

[0028] Spherical mesh:

[0029]

[0030] Ellipsoid mesh:

[0031]

[0032] Where:

[0033] karea(q,k)——the area of ​​a grid on a certain layer of the surface;

[0034] s(k)——surface area of ​​ellipsoidal shell;

[0035] h(k)——the height of a certain layer;

[0036] a, b, c - axis lengths in three directions of the ellipsoid;

[0037] count(k)——the number of surface meshes on a certain layer;

[0038] ——the circumference of the ellipse of the transverse section;

[0039] — the circumference of the longitudinal ellipse perpendicular to the transverse section;

[0040] z(k)——scaling coefficient of the ellipsoid at different heights;

[0041] ——The polar angle of a point in the polar coordinate system;

[0042] Step 5: Based on the particle outflow, calculate the cumulative distribution function (CDF) of direction m, grid q, layer k, energy group g, surface s, and particle type p;

[0043] Step 6: Use the cumulative distribution function CDF to perform source sampling and perform MC calculation;

[0044] To ensure that the geometric positions of the particles obtained by sampling can be evenly distributed on the surface of the sphere or ellipsoid, the height z is first randomly sampled according to the cumulative distribution function of the layer k. The cross-sectional equation is determined based on the height z. For the sphere, its distribution at the boundary of the cross-sectional circular grid is proportional to the central angle θ. For the ellipsoid, in order to simplify the sampling process, it is approximately considered that its distribution at the boundary of the cross-sectional elliptical grid is proportional to the central angle. The sampling formula for the specific geometric position is as follows:

[0045] Spherical:

[0046]

[0047] x0=r z,s cos(θ) (11)

[0048] y0=r z,s sin(θ) (12)

[0049] Ellipsoid:

[0050]

[0051] Where:

[0052] x0, y0, z0——particle position coordinates when sampling geometric positions;

[0053] a, b, c - axis lengths in three directions of the ellipsoid;

[0054] h z，s ——the height of the cross section at a certain layer z of the sphere or ellipsoid;

[0055] ——the lower boundary or upper boundary of a certain height layer of a sphere or ellipsoid;

[0056] a z，s 、b z，s ——the major axis and minor axis of an ellipse of a certain cross section;

[0057] ——Random number generated during sampling, between [0,1];

[0058] θ——center angle, ∈[0,2π];

[0059] ——The corresponding polar angle lower boundary and upper boundary in the polar coordinate system of the surface grid;

[0060] r——radius of the sphere;

[0061] r z,s ——The radius of a certain cross-section circle. Compared with the prior art, the present invention has the following advantages:

[0062] In the method of the present invention, the surface grid and S N The mapping relationship of the calculation grid is used to approximately solve the cumulative distribution function of the surface source. In the shielding calculation of some special geometric structures, the surface coupling calculation can take into account S N High computational efficiency and high precision of MC for complex geometric calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 This is the process of the present invention.

[0064] Figure 2 Schematic diagram of surface mesh decomposition.

[0065] Figure 3 For spherical coupling calculation and S N Compare the results directly. DETAILED DESCRIPTION

[0066] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0067] An S-type method for shielding surface connection N -MC coupling calculation method, in the reactor shielding calculation, the shape of some complex geometric areas is spherical or ellipsoidal, generally S N The computational grid of the method cannot accurately describe this type of coupling surface, and an approximate solution needs to be obtained by using intersecting grid mapping. Taking the sphere as an example, in S N After completing the transport calculation for a slightly larger area than the surface source area, according to S N The z-direction grid is divided into multiple layers of circular rings, and the projection of the rings and S N The intersection of the grid surface, using the intersection coordinates and the height of the ring to obtain the normal vector of a certain position on the surface, and at the same time, the angular flux rate of the SN grid center point where the intersection is located is regarded as the outgoing angular flux rate of the ring mapping area. m,q,g,s The outflow of particle type, surface source, energy group, surface layer, grid, and direction can be obtained, and then the outflow can be converted into a cumulative distribution function as the source input of the MC program to continue the calculation of the area outside the surface.

[0068] like Figure 1 As shown, the method of the present invention comprises the following steps:

[0069] Step 1: Discrete ordinate method S N The transport calculations yield the grid angle flux rate ψ for the reactor geometry including the curved surface source region. m,q,g,s ;

[0070] Step 2: If Figure 2 As shown, Figure 2 The four pictures are sorted by arrows as ①, ②, ③, and ④. N The z-direction layering divides the surface into circular rings; the rings are projected onto S N On the calculated xy plane, as shown in Figure ③, the dotted area is the S where the annular projection passes through. N Grid, while calculating the ring projection and S NThe intersection of the grid surface, according to the number of intersections, the ring belt is divided into the same number of ring grids, according to the intersection points and S N The corresponding relationship between the grid center points is established to obtain the torus grid and S N The mapping relationship of the grid. The specific mapping relationship can be expressed by the ring grid number and S N Grid numbers correspond to:

[0071] Surface mesh numbering format [Number, Layer], S N The grid number format is [i, j, k], where Layer corresponds to k one by one. The number can be obtained by mapping the formula:

[0072] Number=(j-jsuridx1)*(isuridx2-isuridx1+1)+i-isuridx1 (1)

[0073] Where:

[0074] Layer——the layer number corresponding to a certain ring belt;

[0075] Number——surface mesh number on this ring;

[0076] i,j,k——S N The grid numbers used for transport calculations in the x, y, and z directions.

[0077] isuridx1,isuridx2——S N The lower and upper boundary grid numbers for transport calculation in the x-direction;

[0078] jsuridx1——S N The lower boundary grid number for transport calculation in the y direction;

[0079] Step 3: After obtaining the grid mapping relationship, the corresponding S N The angular fluence rate at the center of the grid is regarded as the angular fluence rate emitted from the annular grid, such as Figure 2 As shown in the figure ④; at the same time, since only the outflow along the positive normal direction of the surface source is considered, it is necessary to calculate the normal vectors at different positions on the surface: according to the annular projection calculated in step 2 and S N The intersection of the mesh surface and the center coordinates of the defined surface are used to obtain the normal vector of the intersection point:

[0080] Normal vector of a point on the sphere: [2(x-r_x), 2(y-r_y), 2(z-r_z)],

[0081] Normal vector of a point on the ellipsoid:

[0082] in:

[0083] x,y,z——Surface and S N The coordinates of the grid intersections;

[0084] r_x, r_y, r_z——the coordinates of the center of the sphere and ellipsoid;

[0085] a, b, c - the axis lengths of the ellipsoid in three directions;

[0086] The normal vector and S N The calculated discrete direction m is subjected to vector dot product operation, and if the result is positive, it is determined to be forward emission;

[0087] Step 4: According to the angular fluence rate, the outflow on direction m, grid q, layer k, energy group g, surface s, and particle type p is obtained in sequence;

[0088] The process of calculating the jet from the surface is as follows:

[0089] (1) Obtain the particle flow of a certain layer by summing the outgoing flows of the grids on the same layer;

[0090] (2) For curved surfaces, the normal vector changes with the grid position and needs to be calculated based on the particle flight direction Ω m The dot product with the normal vector at different positions determines whether it is forward emission;

[0091] (3) When calculating the particle flow emitted from a certain grid, the grid area should be the grid area on the curved surface. The derived formula is as follows:

[0092] Spherical mesh:

[0093]

[0094] Ellipsoid mesh:

[0095]

[0096]

[0097] Where:

[0098] karea(q,k)——the area of ​​a grid on a certain layer of the surface;

[0099] s(k)——surface area of ​​ellipsoidal shell;

[0100] h(k)——the height of a certain layer;

[0101] a, b, c - axis lengths in three directions of the ellipsoid;

[0102] count(k)——the number of surface meshes on a certain layer;

[0103] ——the circumference of the ellipse of the transverse section;

[0104] — the circumference of the longitudinal ellipse perpendicular to the transverse section;

[0105] z(k)——scaling coefficient of the ellipsoid at different heights;

[0106] ——The polar angle of a point in the polar coordinate system;

[0107] Step 5: Based on the particle outflow, calculate the cumulative distribution function (CDF) of direction m, grid q, layer k, energy group g, surface s, and particle type p;

[0108] Step 6: Use the cumulative distribution function CDF to perform source sampling and perform MC calculation. To ensure that the geometric positions of the particles obtained by sampling can be evenly distributed on the surface of the sphere or ellipsoid, first randomly sample the height z according to the cumulative distribution function of layer k, and determine the cross-sectional equation based on the height z. For the sphere, its distribution at the boundary of the cross-sectional circular grid is proportional to the central angle θ. For the ellipsoid, in order to simplify the sampling process, it is approximately considered that its distribution at the boundary of the cross-sectional elliptical grid is proportional to the central angle. The sampling formula for the specific geometric position is as follows:

[0109] Spherical:

[0110]

[0111] x=r z,s cos(θ) (11)

[0112] y=r z,s sin(θ) (12)

[0113] Ellipsoid:

[0114]

[0115] Where:

[0116] x0, y0, z0——particle position coordinates when sampling geometric positions;

[0117] a, b, c - axis lengths in three directions of the ellipsoid;

[0118] h z，s ——the height of the cross section at a certain layer z of the sphere or ellipsoid;

[0119] ——the lower boundary or upper boundary of a certain height layer of a sphere or ellipsoid;

[0120] a z，s、b z，s ——the major axis and minor axis of an ellipse of a certain cross section;

[0121] ——Random number generated during sampling, between [0,1];

[0122] θ——center angle, ∈[0,2π];

[0123] ——The corresponding polar angle lower boundary and upper boundary in the polar coordinate system of the surface grid;

[0124] r——radius of the sphere;

[0125] r z,s ——The radius of a certain cross-section circle.

[0126] To verify the effectiveness of the present invention, Figure 3 The shielding problem of a cubic stainless steel point source is given using S N The results of direct calculation by the program and calculation through spherical surface source coupling show that in the area near the surface source, the neutron injection rate error of the two calculation methods is about 1%, which proves the feasibility of the above solution process.

[0127] When using the present invention to solve the shielding problem of some reactors with the above-mentioned special curved surface structure, S N The program calculates shielding areas with simple structures over a large area. Using the MC program to calculate special structural areas through surface sources can reduce the shielding calculation time while ensuring the calculation accuracy of the surface area.

Claims

1. An S method for shielding surface connection N -MC coupling calculation method, characterized by: The steps include: Step 1: Discrete ordinate method S N The transport calculations yield the grid angle flux rate ψ for the reactor geometry including the curved surface source region. m,q,g,s ; Step 2: According to S N The z-direction layering divides the surface into annular zones, and calculates the annular zone projection and S N The intersection of the grid surface, according to the number of intersection points, the ring belt is divided into the same number of ring surface grids, and then the grid on the ring surface and S are obtained. N The mapping relationship of the grid; the specific mapping relationship is the torus grid number and S N Grid numbers correspond to: Torus grid numbering format [Number, Layer], S N The grid number format is [i, j, k], where Layer corresponds to k one by one, and the number is obtained by mapping the formula: Number=(j-jsuridx1)*(isuridx2-isuridx1+1)+i-isuridx1 (1) Where: Layer——the layer number corresponding to a certain ring belt; Number——the torus grid number on a certain annulus; i,j,k——S N The grid numbers for transport calculation division in the x, y, and z directions; isuridx1,isuridx2——S N The lower and upper boundary grid numbers for transport calculation in the x-direction; jsuridx1——S N The lower boundary grid number for transport calculation in the y direction; Step 3: After obtaining the grid mapping relationship, the corresponding S N The angular flux rate at the center of the grid is regarded as the angular flux rate emitted from the annular grid; at the same time, since only the outflow along the positive normal direction of the surface source is considered, the normal vectors at different positions on the surface need to be calculated: according to the annular projection calculated in step 2 and S N The intersection of the mesh surface and the center coordinates of the surface get the normal vector of the intersection: Normal vector of a point on the sphere: [2(x-r_x), 2(y-r_y), 2(z-r_z)], Normal vector of a point on the ellipsoid: in: x,y,z——Surface and S N The coordinates of the grid intersections; r_x, r_y, r_z——the coordinates of the center of the sphere and ellipsoid; a, b, c - the axis lengths of the ellipsoid in three directions; The normal vector and S N The calculated discrete direction m is subjected to vector dot product operation, and if the result is positive, it is determined to be forward emission; Step 4: According to the angular fluence rate, the outflow on direction m, grid q, layer k, energy group g, surface s, and particle type p is obtained in sequence; The process of calculating the jet from the surface is as follows: (1) Obtain the particle flow of a certain layer by summing the outgoing flows of the grids on the same layer; (2) For curved surfaces, the normal vector changes with the grid position and needs to be calculated based on the particle flight direction Ω m The dot product with the normal vector at different positions determines whether it is forward emission; (3) When calculating the particle flow emitted from a certain grid, the grid area should be the grid area on the curved surface. The derived formula is as follows: Spherical mesh: Ellipsoid mesh: Where: karea(q,k)——the area of ​​a grid on a certain layer of the surface; s(k)——surface area of ​​ellipsoidal shell; h(k)——the height of a certain layer; a, b, c - axis lengths in three directions of the ellipsoid; count(k)——the number of surface meshes on a certain layer; ——the circumference of the ellipse of the transverse section; — the circumference of the longitudinal ellipse perpendicular to the transverse section; z(k)——scaling coefficient of the ellipsoid at different heights; ——The polar angle of a point in the polar coordinate system; Step 5: Based on the particle outflow, calculate the cumulative distribution function (CDF) of direction m, grid q, layer k, energy group g, surface s, and particle type p; Step 6: Use the cumulative distribution function CDF to perform source sampling and perform MC calculation. To ensure that the geometric positions of the particles obtained by sampling can be evenly distributed on the surface of the sphere or ellipsoid, first randomly sample the height z according to the cumulative distribution function of layer k, and determine the cross-sectional equation based on the height z. For the sphere, its distribution at the boundary of the cross-sectional circular grid is proportional to the central angle θ. For the ellipsoid, in order to simplify the sampling process, it is approximately considered that its distribution at the boundary of the cross-sectional elliptical grid is proportional to the central angle. The sampling formula for the specific geometric position is as follows: Spherical: x0=r z,s cos(θ) (11) y0=r z,s sin(θ) (12) Ellipsoid: Where: x0, y0, z0——particle position coordinates when sampling geometric positions; a, b, c - axis lengths in three directions of the ellipsoid; h z，s ——the height of the cross section at a certain layer z of the sphere or ellipsoid; ——the lower boundary or upper boundary of a certain height layer of a sphere or ellipsoid; a z，s 、b z，s ——the major axis and minor axis of an ellipse of a certain cross section; ——Random number generated during sampling, between [0,1]; θ——center angle, ∈[0,2π]; ——The corresponding polar angle lower boundary and upper boundary in the polar coordinate system of the surface grid; r——radius of the sphere; r z,s ——The radius of a certain cross-section circle.

Images