Method for analyzing severe accident in nuclear reactor based on advanced particle method

Abstract

A method for analyzing a sever accident in a nuclear reactor based on an advanced particle method includes steps of: 1) performing geometric modeling, setting initial conditions and boundary conditions; 2) updating material physical properties and key parameters; 3) performing mechanical structure module calculation, updating solid particle stress, strain, internal energy, displacement and velocity; 4) performing thermal hydraulic module calculation, updating fluid particle internal energy, position and velocity; 5) performing chemical reaction module calculation, updating particle matter composition and internal energy; 6) performing neutron physics module calculation, updating particle neutron flux density; and 7) outputting data. The method of the present invention is based on the discrete form of the advanced particle method, which is capable of accurately capturing cross-sectional changes, matter changes, and phase changes. Compared with grid method, the present invention can effectively avoid a mesh distortion problem existing in a large deformation.

Claims

1. A method for analyzing a sever accident in a nuclear reactor based on an advanced particle method, comprising steps of: step 1: according to a nuclear reactor severe accident analysis calculation object, constructing a particle geometry model, wherein an initial state of each particle is determined by an actual calculation object; setting particle property parameters, enthalpies, velocities, positions and stress-strain states, and setting boundary conditions according to the actual calculation object, comprising pressure boundaries, temperature boundaries, heat flow boundaries, pressure gradient boundaries, velocity boundaries, load boundaries, symmetry boundaries and period boundaries; wherein the particle geometry model, key parameters, the initial state and the boundary conditions of the nuclear reactor severe accident analysis calculation object are established in the step 1, thereby restoring a real state of an arbitrarily complex nuclear reactor severe accident analysis calculation object; step 2: calculating material property changes for each particle with a nuclear reactor material property library, comprising densities, specific heat capacities, thermal conductivities, melting points, boiling points, latent heat of melting, latent heat of vaporization, thermal conductivities, Young's modulus, Poisson's ratios, thermal expansion coefficients and thermal creep coefficients; and updating the key parameters, comprising particle number densities, particle neighborhood sets and time steps required by the advanced particle method; wherein the material properties and the key parameters of each particle is updated in the step 2; material property update is used to obtain changes in reactor type material properties during the severe accident in real time, and key parameter update ensures necessary conditions required for accuracy of the advanced particle method, which provide conditions and supports for subsequent nuclear reactor severe accident analysis calculations; step 3: considering possible changes in mechanical structures during the severe accident in the nuclear reactor and performing mechanical structure calculations, comprising thermal expansion, elastic deformation, plastic deformation, creep and fracture calculations; calculating a thermal expansion strain according to an equation (1): $\begin{array}{cc}{\left[d{s}^T\right]}_i={\kappa }_i\left(T_i-T_i^{\mathrm{ref}}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)& \mathrm{equation}\left(1\right)\end{array}$ wherein [ds.sup.T].sub.i is a thermal expansion stress increment tensor of a particle i, N/m.sup.2; κ.sub.i is a thermal expansion coefficient of the particle i; T.sub.i is a temperature of the particle i, K; T.sub.i.sup.ref is a reference temperature of the particle i, K; calculating an elastic stress according to an equation (2):
σ.sub.αβ.sup.E=λε.sub.yy.sup.Eδ.sub.αβ+2με.sub.αβ.sup.E   equation (2) wherein σ.sub.αβ.sup.E a component of an elastic stress tensor, N/m.sup.2; λ is a first parameter of a Lamé constant, $\lambda =\frac{E\upsilon }{\left(1+\upsilon \right)\left(1-2\upsilon \right)};$ μ is a second parameter of the Lamé constant, $\mu =\frac{E}{2\left(1+\upsilon \right)};$ wherein E is the Young's modulus and ν is the Poisson's ratio; ε.sub.yy.sup.E is a sum of diagonal terms of the elastic strain tensor, ε.sub.yy.sup.E=ε.sub.11.sup.E+ε.sub.22.sup.E+ε.sub.33.sup.E; wherein ε.sub.11.sup.E, ε.sub.22.sup.E, ε.sub.33.sup.E are diagonal terms of first, second and third rows of the elastic strain tensor; δ.sub.αβ is a Kronecker function, ${\delta }_{\alpha \beta }=\{\begin{array}{cc}1& \alpha =\beta \\ 0& \alpha \ne \beta \end{array};$ ε.sub.αβ.sup.E is a component of the elastic strain tensor; α is α a direction, being any value in x, y, z; β is a β direction, being any value in x, y, z; if an equivalent force of the particle is greater than a yield limit, considering that the plastic deformation occurs, and calculating plastic stress-strain according to an equation (3) and an equation (4). $\begin{array}{cc}{\left[d{\epsilon }^p\right]}_n={\left[d\epsilon \right]}_n-{\left[d{\epsilon }^T\right]}_n-\frac{{\left[d\epsilon \right]}_n-{\left[d{\epsilon }^T\right]}_n-{\left[s\right]}_{n-1}d{\lambda }_n}{1+2\mu d{\lambda }_n}& \mathrm{equation}\left(3\right)\end{array}$ $\begin{array}{cc}{\left[d\epsilon \right]}_n={\left[d{\epsilon }^p\right]}_n+{\left[d{\epsilon }^e\right]}_n+{\left[d{\epsilon }^T\right]}_n& \mathrm{equation}\left(4\right)\end{array}$ wherein [dε.sup.P].sup.n is a plastic strain increment tensor at a time step n, N/m.sup.2; [dε].sub.n is a strain increment tensor at the time step n, N/m.sup.2; [dε.sup.T].sub.n is a thermal expansion strain increment tensor at the time step n, N/m.sup.2; [dε.sup.e].sub.n is an elastic strain increment tensor at the time step n, N/m.sup.2; [s].sub.n−1 is a strain tensor at a time step n−1, N/m.sup.2; dλ.sub.n is an incremental parameter, determined by material mechanical properties; determining a creep calculation by a creep rate, and calculating according to an equation (5)
=ωσϕ  equation (5) wherein is the creep rate; ω is a thermal creep coefficient; σ is a stress tensor, N/m.sup.2; ϕ is a fast neutron flux density, neutron/m.sup.2/s; wherein fracture is determined according to an inter-particle stress size and a particle spacing; when the inter-particle stress is greater than a fracture stress threshold or when the particle spacing is greater than a tensile fracture threshold or when the particle spacing is less than a compression fracture threshold, the fracture is considered to occur; there is no solid stress between fractured particles, and an interaction therebetween is transformed into a collision interaction; wherein stress or strain of each particle under thermal expansion, elasticity, plasticity, creep and fracture is calculated in the step 3, and then the velocities, the positions and energy of the particle are calculated by mass conservation, energy conservation and momentum conservation equations; the three conservation equations are in a same form; a scattering term of the stress is calculated by the advanced particle method in a discrete form, and the advanced particle method in the discrete form is shown in equations (6) to (12): $\begin{array}{cc}D{\varphi }_i=H_{s}M{b}_i& \mathrm{equation}\left(6\right)\end{array}$ $\begin{array}{cc}D={\left[\frac{\partial }{\partial x}\frac{\partial }{\partial y}\frac{\partial }{\partial z}\frac{{\partial }^2}{\partial x^2}\frac{{\partial }^2}{\partial y^2}\frac{{\partial }^2}{\partial z^2}\frac{{\partial }^2}{\partial x\partial y}\frac{{\partial }^2}{\partial x\partial z}\frac{{\partial }^2}{\partial y\partial z}\right]}^T& \mathrm{equation}\left(7\right)\end{array}$ $\begin{array}{cc}{H}_s=\mathrm{diag}\left(\frac{1}{n^0}\frac{1}{n^0}\frac{1}{n^0}\frac{2}{n^0{l}_0}\frac{2}{n^0{l}_0}\frac{2}{n^0{l}_0}\frac{1}{n^0{l}_0}\frac{1}{n^0{l}_0}\frac{1}{n^0{l}_0}\right)& \mathrm{equation}\left(8\right)\end{array}$ $\begin{array}{cc}{M}^{-1}=C=\underset{j\ne i}{.Math.}\frac{w_{\mathrm{ij}}{p}_{\mathrm{ij}}.Math.p_{\mathrm{ij}}}{n^0}& \mathrm{equation}\left(9\right)\end{array}$ $\begin{array}{cc}{b}_i=\underset{j\ne i}{.Math.}{w}_{\mathrm{ij}}{\varphi }_{\mathrm{ij}}{p}_{\mathrm{ij}}& \mathrm{equation}\left(10\right)\end{array}$ $\begin{array}{cc}{p}_{\mathrm{ij}}=\{\begin{array}{cc}{\left[\frac{x_{\mathrm{ij}}}{r_{\mathrm{ij}}}\frac{y_{\mathrm{ij}}}{r_{\mathrm{ij}}}\frac{z_{\mathrm{ij}}}{r_{\mathrm{ij}}}\frac{x_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\frac{y_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\frac{z_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\frac{x_{\mathrm{ij}}{y}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}\frac{x_{\mathrm{ij}}{z}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}\frac{y_{\mathrm{ij}}{z}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}\right]}^T& j\in \mathrm{Internal}\mathrm{or}\mathrm{Dirichlet}\\ {\left[n_x{n}_y{n}_z\frac{2n_x{x}_{\mathrm{ij}}}{l_0}\frac{2n_y{y}_{\mathrm{ij}}}{l_0}\frac{2n_z{z}_{\mathrm{ij}}}{l_0}\frac{\left(n_y{x}_{\mathrm{ij}}+n_x{y}_{\mathrm{ij}}\right)}{\left(l_0\right)}\frac{\left(n_z{x}_{\mathrm{ij}}+n_x{z}_{\mathrm{ij}}\right)}{\left(l_0\right)}\frac{\left(n_z{y}_{\mathrm{ij}}+n_y{z}_{\mathrm{ij}}\right)}{\left(l_0\right)}\right]}^T& j\in \mathrm{Neumann}\end{array}& \mathrm{equation}\left(11\right)\end{array}$ $\begin{array}{cc}{\varphi }_{\mathrm{ij}}=\{\begin{array}{cc}\frac{{\varphi }_j-{\varphi }_i}{r_{\mathrm{ij}}}& j\in \mathrm{Internal}\mathrm{or}\mathrm{Dirichlet}\\ \frac{\partial {\varphi }_j}{\partial n}& j\in \mathrm{Neumann}\end{array}& \mathrm{equation}\left(12\right)\end{array}$ wherein D is a second order differential operator as shown in the equation (7); H.sub.s is a coefficient diagonal matrix as shown in the equation (8); M is a correction matrix as shown in the equation (9); b.sub.i is a correction parameter vector as shown in the equation (10); M.sup.−1 is an inverse matrix of the correction matrix M; C is a coefficient matrix, which is same as the inverse matrix of the correction matrix M; x is an x direction; y is a y direction; z is a z direction; n.sup.0 is an initial particle number density; l.sub.0 is an initial particle spacing, m; diag is diagonal matrix compliance; w.sub.ij is a kernel function value between the particle i and a particle j; x.sub.ij is a distance between the particle i and the particle j in the x direction, m; y.sub.ij is a distance between the particle i and the particle j in the y direction, m; z.sub.ij is a distance between the particle i and the particle j in the z direction, m; r.sub.ij a distance between the particle i and the particle j, m; n.sub.x is a component of a normal vector of the particle j in the x direction; n.sub.y is a component of the normal vector of the particle j in the y direction; n.sub.z is a component of the normal vector of the particle j in the z direction; ϕ.sub.j is an arbitrary parameter value of the particle j, which is a vector or a scalar; ϕ.sub.i is an arbitrary parameter value of the particle i, which is a vector or a scalar; $\frac{\partial {\varphi }_j}{\partial n}$ is a partial derivative of an arbitrary parameter of the particle j in a direction of the normal vector of the particle j; Internal is an internal particle; Dirichlet is is a fixed-value boundary condition, comprising the pressure boundaries, the temperature boundaries and the velocity boundaries; Neumann is a gradient boundary condition, comprising the heat flow boundaries, the pressure gradient boundaries, and the load boundaries; wherein macroscopic state changes in velocity, position and energy of the actual calculation object under the thermal expansion, the elasticity, the plasticity, the creep and the fracture are obtained; step 4: performing a thermal hydraulic calculation, comprising fluid motions under gravity, viscous forces, pressure and surface tensions, and fluid energy changes during heat transfer; calculating the fluid motions according to an equation (13), so as to obtain velocity and position changes of fluid particles; $\begin{array}{cc}\rho \frac{du}{dt}=-\nabla P+{\mu }_f{\nabla }^{2}u+\rho f+\rho g& \mathrm{equation}\left(13\right)\end{array}$ wherein t is time, s; ρ is a density, kg/m.sup.3; u is a velocity vector, m/s; P is a pressure, Pa; μ.sub.f is a dynamic viscosity, Pa.Math.s; f is a surface tension vector, N/kg; g is a gravitational acceleration vector, m/s.sup.2; wherein a gradient term of the pressure and a velocity Laplace term of a viscosity term uses the advanced particle method in the discrete form as shown in the equations (6) to (12), and the pressure is calculated by implicit iterative solution of a pressure Poisson equation, as shown in an equation (14): $\begin{array}{cc}\frac{1}{\rho }{\left({\nabla }^2{P}^{k+1}\right)}_i=\frac{1-\xi }{\Delta t}\nabla .Math.u^{*}-\frac{\xi }{\Delta t^2}\left(\frac{n^{*}-n^0}{n^0}\right)& \mathrm{equation}\left(14\right)\end{array}$ wherein n* is a temporary particle number density, which is calculated from particle position obtained by calculating a gravity term, the viscosity term and a surface tension term for the particle; ξ is a pressure Poisson's equation weighting coefficient, ranging from 0 to 1; Δt is a time step, s; u* is a temporary velocity vector, m/s; P.sup.k+1 is a pressure value at a time step k+1, Pa; ∇is a Hamiltonian arithmetic; wherein during a discretization process of the advanced particle method, a dynamic viscosity in the viscosity term adopts a harmonic average of the dynamic viscosity, as shown in an equation (15) $\begin{array}{cc}{\mu }_{\mathrm{ij}}=\frac{2{\mu }_i{\mu }_j}{{\mu }_i+{\mu }_j}& \mathrm{Equation}\left(15\right)\end{array}$ wherein μ.sub.ij is a dynamic viscosity between the particle i and the particle j, Pa.Math.s; μ.sub.i is a dynamic viscosity of the particle i, Pa.Math.s; μ.sub.j is a dynamic viscosity of the particle j, Pa.Math.s; calculating the surface tension with a free energy model based on a surface tension model, as shown in an equation (16)
f=F(r.sub.ij−r.sub.min)(r.sub.ij−r.sub.e)/m   equation (16) wherein f is a surface tension vector, N/kg; m is a mass, kg; F is a free energy coefficient; r.sub.min is a minimum distance of the particle i from surrounding particles, adopting 1.5 l.sub.0; r.sub.e is a particle radius of action; wherein a fluid heat transfer process calculation comprises radiation, heat conduction, convection, heat flow boundaries and chemical heat as shown in an equation (17), from which energy changes of the fluid particles are calculated; $\begin{array}{cc}\rho \frac{\partial H}{\partial t}={\rho \left(\frac{\partial H}{\partial t}\right)}_{\mathrm{Trans}}+{\rho \left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}+Q_{\mathrm{heatflow}}+Q_{\mathrm{chem}}& \mathrm{equation}\left(17\right)\end{array}$ wherein H is an enthalpy, J/kg; ρ is a density, kg/m.sup.3; Q.sub.heatflow is a volume heat flow boundary, W/m.sup.3; Q.sub.chem is chemical heat, W/m.sup.3; ${\left(\frac{\partial H}{\partial t}\right)}_{Trans}$ is a partial derivative of the enthalpy against time due to the heat conduction and convective heat exchange, W/kg; ${\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}$ is a partial derivative of the enthalpy against time due to radiative heat exchange, W/kg; wherein a partial derivative calculation of the enthalpy against time due to the heat conduction and the convective heat transfer uses a thermal conductivity differential equation, as shown in an equation (18); in the advanced particle method, the heat conduction and the convection use a same set of models; for the heat conduction, the particles do not move; and for the convection, the particles move; a temperature Laplace term in the equation (18) uses the advanced particle method in the discrete form as shown in the equations (6) to (12), the thermal conductivity during a discrete process adopts a harmonic average; $\begin{array}{cc}{\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{Trans}}=\frac{k}{\rho }{\nabla }^{2}T& \mathrm{equation}\left(18\right)\end{array}$ wherein k is the thermal conductivity, W/(m.Math.K); T is the temperature, K; calculating the radiative heat exchange according to an equation (19): $\begin{array}{cc}{\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}=\frac{{\mathrm{\epsilon \sigma }}_s\left(T_i^4-T_j^4\right)}{\rho l_0}& \mathrm{equation}\left(19\right)\end{array}$ wherein ε is an emissivity; σ.sub.s is a Stefan-Boltzmann constant; T.sub.i is the temperature of the particle i, K; T.sub.j is a temperature of the particle j, K; determining the heat flow boundaries according to an actual situation, comprising heat source and the cooling boundaries; determining the chemical heat by chemical reactions; wherein changes of fluid particle velocities, positions and energy are calculated in the step 4, which truly reflect fluid motion and temperature field changes of the calculation object during the serious accident in an actual nuclear reactor; step 5: performing a chemical reaction calculation, comprising redox reactions, eutectic reactions and corrosion phenomena, wherein the redox reactions determine a material change rate based on a reaction rate, the eutectic reactions are determined based on a diffusion rate, and the corrosion phenomena are determined by a corrosion rate; the chemical reaction calculation relays on a nuclear reactor material database; wherein the step 5 calculates the material component changes of a core material under the chemical reactions, and truly reflects chemical reaction effects on the severe accident in the nuclear reactor through the material property changes; step 6: performing a neutron physics calculation using a multigroup approximation S.sub.N difference Boltzmann transport equation as shown in an equation (20) $\begin{array}{cc}\Omega .Math.\nabla \varphi \left(r,\Omega ,E_n\right)+{.Math.}_t\varphi \left(r,\Omega ,E_n\right)=\int \int {.Math.}_s\left(r,{\Omega }^{\prime },E_n^{\prime }.fwdarw.\Omega ,E_n\right)\varphi \left(r,{\Omega }^{\prime },E_n^{\prime }\right)d{\Omega }^{\prime }{\mathrm{dE}}_n^{\prime }+\frac{\chi \left(r,E_n\right)}{4\pi }\int \int v{.Math.}_f\left(r,{\Omega }^{\prime },E_n^{\prime }\right)\varphi \left(r,{\Omega }^{\prime },E_n^{\prime }\right)d{\Omega }^{\prime }d{E}_n^{\prime }+\frac{Q\left(r,E_n\right)}{4\pi }& \mathrm{equation}\left(20\right)\end{array}$ wherein Ω is a direction vector; Ω′ is another direction vector, but different from Ω; ϕ(r,Ω,E.sub.n) is a neutron angular flux density when an input is r,Ω,E.sub.n; ϕ(r,Ω′,E′.sub.n) is a neutron angular flux density when the input is r,Ω′,E′.sub.n; Σ.sub.t is a total neutron cross section; Q(r,E.sub.n) is neutron source intensity; E.sub.n is neutron energy; E′.sub.n is another neutron energy, but different from E.sub.n; Σ.sub.s(r,Ω′,E′.sub.n.fwdarw.Ω,E.sub.n) is a scattering cross section; χ(r,E.sub.n) is a fission spectrum; ν is a quantity of neutrons released per fission; Σ.sub.f(r,Ω′,E′.sub.n) is a neutron fission cross section; selecting a multi-group core cross-section database according to an actual calculation reactor type; arranging a structured grid on a core particle geometry arrangement, and adopting a coarse mesh finite difference method to accelerate solution; using a particle mesh mapping technology to realize grid-particle information interaction; wherein the step 6 calculates the neutron angular flux density in the core, and changes a heat source distribution in the core through the neutron angular flux density to change the material physical properties and material stress strain, thus realizing coupling analysis of neutron physics and thermal engineering hydraulics; and step 7: outputting required data; determining whether an end-of-calculation condition is met, if not, advancing the time step and returning to the step 2, if yes, ending calculation; wherein the steps 1-7 are used to analyze the severe accident in the nuclear reactor, which integrate possible key factors during the severe accident in the nuclear reactor, comprising the mechanical structure changes, the fluid motions, heat transfer phase changes, the chemical reactions and the neutron physics, so as to analyze key phenomena of the severe accident in the nuclear reactor, mainly comprising: core heating transients, core melting, melt migration, debris bed behavior, core melt retention behavior, melt-concrete interaction, and melt-coolant interaction.

2. The method, as recited in claim 1, which is based on the discrete form of the advanced particle method for accurately capturing cross-sectional changes, matter changes, and phase changes.

3. The method, as recited in claim 1, which is based on the discrete form of the advanced particle method to effectively avoid a mesh distortion problem existing in a large deformation.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0146] FIGURE is a flow chart of a method for analyzing a sever accident in a nuclear reactor based on an advanced particle method according to the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

[0147] Referring to the accompanying drawing and embodiment, the present invention will be further illustrated.

[0148] The present invention provides a method for analyzing a sever accident in a nuclear reactor based on an advanced particle method, as shown in FIGURE. High-temperature melting process analysis of a single fuel rod in a typical pressurized water reactor under simplified conditions is taken as an example, comprising steps of: [0149] step 1: according to a single fuel rod in a typical pressurized water reactor, constructing a particle geometry model, wherein a simplified fuel rod considers only a UO.sub.2 core block and a Zr-4 alloy cladding; the fuel rod in a bare state with no additional stress strain; the core block is set as an internal heat source, a size of which is calculated based on a decay power; [0150] step 2: calculating material property changes for the UO.sub.2 core block and the Zr-4 alloy cladding with a nuclear reactor material property library, comprising densities, specific heat capacities, thermal conductivities, melting points, boiling points, latent heat of melting, latent heat of vaporization, thermal conductivities, Young's modulus, Poisson's ratios, thermal expansion coefficients and thermal creep coefficients; and updating the key parameters, comprising particle number densities, particle neighborhood sets and time steps required by the advanced particle method; [0151] step 3: considering possible changes in mechanical structures during the severe accident in the nuclear reactor and performing mechanical structure calculations, comprising thermal expansion, elastic deformation, plastic deformation, creep and fracture calculations; [0152] calculating a thermal expansion strain according to an equation (1):

[00017]$\begin{array}{cc}{\left[d{s}^T\right]}_i={\kappa }_i\left(T_i-T_i^{ref}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)& \mathrm{equation}\left(1\right)\end{array}$ [0153] wherein [0154] [ds.sup.T].sub.i is a thermal expansion stress increment tensor of a particle i, N/m.sup.2; [0155] κ.sub.i is a thermal expansion coefficient of the particle i; [0156] T.sub.i is a temperature of the particle i, K; [0157] T.sub.i.sup.ref is a reference temperature of the particle i, K; [0158] calculating an elastic stress according to an equation (2):

σ.sub.αβ.sup.E=λε.sub.yy.sup.Eδ.sub.αβ+2με.sub.αβ.sup.E   equation (2) [0159] wherein [0160] σ.sub.αβ.sup.E a component of an elastic stress tensor, N/m.sup.2; [0161] λ is a first parameter of a Lamé constant,

[00018]$\lambda =\frac{E\upsilon }{\left(1+\upsilon \right)\left(1-2\upsilon \right)};$ [0162] μ is a second parameter of the Lamé constant,

[00019]$\mu =\frac{E}{2\left(1+\upsilon \right)};$

wherein E is the Young's modulus and ν is the Poisson's ratio; [0163] ε.sub.yy.sup.E is a sum of diagonal terms of the elastic strain tensor, ε.sub.yy.sup.E=ε.sub.11.sup.E+ε.sub.22.sup.E+ε.sub.33.sup.E; wherein ε.sub.11.sup.E, ε.sub.22.sup.E, ε.sub.33.sup.E are diagonal terms of first, second and third rows of the elastic strain tensor; [0164] δ.sub.αβ is a Kronecker function,

[00020]${\delta }_{\alpha \beta }=\{\begin{array}{cc}1& \alpha =\beta \\ 0& \alpha \ne \beta \end{array};$ [0165] ε.sub.αβ.sup.E is a component of the elastic strain tensor; [0166] α is α a direction, being any value in x, y, z; [0167] β is a β direction, being any value in x, y, z; [0168] if an equivalent force of the particle is greater than a yield limit, considering that the plastic deformation occurs, and calculating plastic stress-strain according to an equation (3) and an equation (4).

[00021]$\begin{array}{cc}{\left[d{\epsilon }^p\right]}_n={\left[d\epsilon \right]}_n-{\left[d{\epsilon }^T\right]}_n-\frac{{\left[d\epsilon \right]}_n-{\left[d{\epsilon }^T\right]}_n-{\left[s\right]}_{n-1}d{\lambda }_n}{1+2\mu d{\lambda }_n}& \mathrm{equation}\left(3\right)\end{array}$$\begin{array}{cc}\left[d\epsilon \right]={\left[d{\epsilon }^p\right]}_n+{\left[d{\epsilon }^e\right]}_n+{\left[d{\epsilon }^T\right]}_n& \mathrm{equation}\left(4\right)\end{array}$ [0169] wherein [0170] [dε.sup.P].sup.n is a plastic strain increment tensor at a time step n, N/m.sup.2; [0171] [dε].sub.n is a strain increment tensor at the time step n, N/m.sup.2; [0172] [dε.sup.T].sub.n is a thermal expansion strain increment tensor at the time step n, N/m.sup.2; [0173] [dε.sup.e].sub.n is an elastic strain increment tensor at the time step n, N/m.sup.2; [0174] [s].sub.n−1 is a strain tensor at a time step n−1, N/m.sup.2; [0175] μ is a second parameter of the Lamé constant,

[00022]$\mu =\frac{E}{2\left(1+\upsilon \right)};$ [0176] dλ.sub.n is an incremental parameter, determined by material mechanical properties; [0177] determining a creep calculation by a creep rate, and calculating according to an equation (5)

=ωσϕ  equation (5) [0178] wherein [0179] is the creep rate; [0180] ω is a thermal creep coefficient; [0181] σ is a stress tensor, N/m.sup.2; [0182] ϕ is a fast neutron flux density, neutron/m.sup.2/s; [0183] wherein fracture is determined according to an inter-particle stress size and a particle spacing; when the inter-particle stress is greater than a fracture stress threshold or when the particle spacing is greater than a tensile fracture threshold or when the particle spacing is less than a compression fracture threshold, the fracture is considered to occur; there is no solid stress between fractured particles, and an interaction therebetween is transformed into a collision interaction; [0184] wherein stress or strain of each particle under thermal expansion, elasticity, plasticity, creep and fracture is calculated in the step 3, and then the velocities, the positions and energy of the particle are calculated by mass conservation, energy conservation and momentum conservation equations; the three conservation equations are common and in the same form, which will not be listed here; a scattering term of the stress is calculated by the advanced particle method in a discrete form, and the advanced particle method in the discrete form is shown in equations (6) to (12):

[00023]$\begin{array}{cc}D{\varphi }_i=H_s{\mathrm{Mb}}_i& \mathrm{equation}\left(6\right)\end{array}$$\begin{array}{cc}D=[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\frac{\partial }{\partial x}& \frac{\partial }{\partial y}\end{array}& \frac{\partial }{\partial z}\end{array}& \frac{{\partial }^2}{\partial x^2}\end{array}& \frac{{\partial }^2}{\partial y^2}\end{array}& \frac{{\partial }^2}{\partial z^2}\end{array}& \frac{{\partial }^2}{\partial x\partial y}\end{array}& \frac{{\partial }^2}{\partial x\partial z}\end{array}& {\frac{{\partial }^2}{\partial y\partial z}]}^T\end{array}& \mathrm{equation}\left(7\right)\end{array}$$\begin{array}{cc}{H}_s=\mathrm{diag}\left(\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\frac{1}{n^0}& \frac{1}{n^0}\end{array}& \frac{1}{n^0}\end{array}& \frac{2}{n^0{l}_0}\end{array}& \frac{2}{n^0{l}_0}\end{array}& \frac{2}{n^0{l}_0}\end{array}& \frac{1}{n^0{l}_0}\end{array}& \frac{1}{n^0{l}_0}\end{array}& \frac{1}{n^0{l}_0}\end{array}\right)& \mathrm{equation}\left(8\right)\end{array}$$\begin{array}{cc}{M}^{-1}=C=\underset{j\ne i}{.Math.}\frac{w_{\mathrm{ij}}{p}_{\mathrm{ij}}.Math.p_{\mathrm{ij}}}{n^0}& \mathrm{equation}\left(9\right)\end{array}$$\begin{array}{cc}{b}_i=\underset{j\ne i}{.Math.}{w}_{\mathrm{ij}}{\varphi }_{\mathrm{ij}}{p}_{\mathrm{ij}}& \mathrm{equation}\left(10\right)\end{array}$$\begin{array}{cc}{p}_{\mathrm{ij}}=\{\begin{array}{cc}[\begin{array}{ccc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\frac{x_{\mathrm{ij}}}{r_{\mathrm{ij}}}& \frac{y_{\mathrm{ij}}}{r_{\mathrm{ij}}}\end{array}& \frac{z_{\mathrm{ij}}}{r_{\mathrm{ij}}}\end{array}& \frac{x_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\end{array}& \frac{y_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\end{array}& \frac{z_{\mathrm{ij}}^2}{l_0{r}_{\mathrm{ij}}}\end{array}& \frac{x_{\mathrm{ij}}{y}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}\end{array}& \frac{x_{\mathrm{ij}}{z}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}& {\frac{y_{\mathrm{ij}}{z}_{\mathrm{ij}}}{l_0{r}_{\mathrm{ij}}}]}^T\end{array}& j\in \mathrm{Internal}\mathrm{or}\mathrm{Dirichlet}\\ \left[\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}{n}_x& n_y\end{array}& n_z\end{array}& \frac{2n_x{x}_{\mathrm{ij}}}{l_0}\end{array}& \frac{2n_y{y}_{\mathrm{ij}}}{l_0}\end{array}& \frac{2n_z{z}_{\mathrm{ij}}}{l_0}\end{array}& \frac{\begin{array}{c}{n}_y{x}_{\mathrm{ij}}+\\ n_x{y}_{\mathrm{ij}}\end{array}}{l_0}\end{array}& \frac{\begin{array}{c}{n}_z{x}_{\mathrm{ij}}+\\ n_x{z}_{\mathrm{ij}}\end{array}}{l_0}\end{array}& \frac{\begin{array}{c}{n}_z{x}_{\mathrm{ij}}+\\ n_x{z}_{\mathrm{ij}}\end{array}}{l_0}\end{array}\right]& j\in \mathrm{Neumann}\end{array}& \mathrm{equation}\left(11\right)\end{array}$$\begin{array}{cc}{\varphi }_{\mathrm{ij}}=\{\begin{array}{cc}\frac{{\varphi }_j-{\varphi }_i}{r_{\mathrm{ij}}}& j\in \mathrm{Internal}\mathrm{or}\mathrm{dirichlet}\\ \frac{\partial {\varphi }_j}{\partial n}& j\in \mathrm{Neumann}\end{array}& \mathrm{equation}\left(12\right)\end{array}$ [0185] wherein [0186] D is a second order differential operator as shown in the equation (7); [0187] H.sub.s is a coefficient diagonal matrix as shown in the equation (8); [0188] M is a correction matrix as shown in the equation (9); [0189] b.sub.i is a correction parameter vector as shown in the equation (10); [0190] M.sup.−1 is an inverse matrix of the correction matrix M; [0191] C is a coefficient matrix, which is same as the inverse matrix of the correction matrix M; [0192] x is an x direction; [0193] y is a y direction; [0194] z is a z direction; [0195] n.sup.0 is an initial particle number density; [0196] l.sub.0 is an initial particle spacing, m; [0197] diag is diagonal matrix compliance; [0198] w.sub.ij is a kernel function value between the particle i and a particle j; [0199] x.sub.ij is a distance between the particle i and the particle j in the x direction, m; [0200] y.sub.ij is a distance between the particle i and the particle j in the y direction, m; [0201] z.sub.ij is a distance between the particle i and the particle j in the z direction, m; [0202] r.sub.ij a distance between the particle i and the particle j, m; [0203] n.sub.x is a component of a normal vector of the particle j in the x direction; [0204] n.sub.y is a component of the normal vector of the particle j in the y direction; [0205] n.sub.z is a component of the normal vector of the particle j in the z direction; [0206] ϕ.sub.j is an arbitrary parameter value of the particle j, which is a vector or a scalar; [0207] ϕ.sub.i is an arbitrary parameter value of the particle i, which is a vector or a scalar;

[00024]$\frac{\partial {\varphi }_j}{\partial n}$

is a partial derivative of an arbitrary parameter of the particle j in a direction of the normal vector of the particle j; [0208] Internal is an internal particle; [0209] Dirichlet is is a fixed-value boundary condition, comprising the pressure boundaries, the temperature boundaries and the velocity boundaries; [0210] Neumann is a gradient boundary condition, comprising the heat flow boundaries, the pressure gradient boundaries, and the load boundaries; [0211] wherein macroscopic state changes in velocity, position and energy of the actual calculation object under the thermal expansion, the elasticity, the plasticity, the creep and the fracture are obtained; [0212] step 4: performing a thermal hydraulic calculation, comprising fluid motions under gravity, viscous forces, pressure and surface tensions, and fluid energy changes during heat transfer; [0213] calculating the fluid motions according to an equation (13), so as to obtain velocity and position changes of fluid particles;

[00025]$\begin{array}{cc}\rho \frac{du}{dt}=-\nabla P+{\mu }_f{\nabla }^{2}u+\rho f+\rho g& \mathrm{equation}\left(13\right)\end{array}$ [0214] wherein [0215] t is time, s; [0216] ρ is a density, kg/m.sup.3; [0217] u is a velocity vector, m/s; [0218] P is a pressure, Pa; [0219] μ.sub.f is a dynamic viscosity, Pa.Math.s; [0220] f is a surface tension vector, N/kg; [0221] g is a gravitational acceleration vector, m/s.sup.2; [0222] wherein a gradient term of the pressure and a velocity Laplace term of a viscosity term uses the advanced particle method in the discrete form as shown in the equations (6) to (12), and the pressure is calculated by implicit iterative solution of a pressure Poisson equation, as shown in an equation (14):

[00026]$\begin{array}{cc}\frac{1}{\rho }{\left({\nabla }^2{P}^{k+1}\right)}_i=\frac{1-\xi }{\Delta t}\nabla .Math.u^{*}-\frac{\xi }{\Delta t^2}\left(\frac{n^{*}-n^0}{n^0}\right)& \mathrm{equation}\left(14\right)\end{array}$ [0223] wherein [0224] n is a temporary particle number density, which is calculated from particle position obtained by calculating a gravity term, the viscosity term and a surface tension term for the particle; [0225] ξ is a pressure Poisson's equation weighting coefficient, ranging from 0 to 1; [0226] Δt is a time step, s; [0227] u* is a temporary velocity vector, m/s; [0228] P.sup.k+1 is a pressure value at a time step k+1, Pa; [0229] ∇is a Hamiltonian arithmetic; [0230] wherein during a discretization process of the advanced particle method, a dynamic viscosity in the viscosity term adopts a harmonic average of the dynamic viscosity, as shown in an equation (15)

[00027]$\begin{array}{cc}{\mu }_{\mathrm{ij}}=\frac{2{\mu }_i{\mu }_j}{{\mu }_i+{\mu }_j}& \mathrm{Equation}\left(15\right)\end{array}$ [0231] wherein [0232] μ.sub.ij is a dynamic viscosity between the particle i and the particle j, Pa.Math.s; [0233] μ.sub.i is a dynamic viscosity of the particle i, Pa.Math.s; [0234] μ.sub.j is a dynamic viscosity of the particle j, Pa.Math.s; [0235] calculating the surface tension with a free energy model based on a surface tension model, as shown in an equation (16)

f=F(r.sub.ij−r.sub.min)(r.sub.ij−r.sub.e)/m   equation (16) [0236] wherein [0237] f is a surface tension vector, N/kg; [0238] m is amass, kg; [0239] F is a free energy coefficient; [0240] r.sub.min is a minimum distance of the particle i from surrounding particles, adopting 1.5 l.sub.0; [0241] r.sub.e is a particle radius of action; [0242] wherein a fluid heat transfer process calculation comprises radiation, heat conduction, convection, heat flow boundaries and chemical heat as shown in an equation (17), from which energy changes of the fluid particles are calculated;

[00028]$\begin{array}{cc}\rho \frac{\partial H}{\partial t}={\rho \left(\frac{\partial H}{\partial t}\right)}_{\mathrm{Trans}}+{\rho \left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}+Q_{\mathrm{heatflow}}+Q_{\mathrm{chem}}& \mathrm{equation}\left(17\right)\end{array}$ [0243] wherein [0244] H is an enthalpy, J/kg; [0245] t is time, s; [0246] ρ is a density, kg/m.sup.3; [0247] Q.sub.heatflow is a volume heat flow boundary, W/m.sup.3; [0248] Q.sub.chem is chemical heat, W/m.sup.3;

[00029]${\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{Trans}}$

is a partial derivative of the enthalpy against time due to the heat conduction and convective heat exchange, W/kg;

[00030]${\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}$

is a partial derivative of the enthalpy against time due to radiative heat exchange, W/kg; [0249] wherein a partial derivative calculation of the enthalpy against time due to the heat conduction and the convective heat transfer uses a thermal conductivity differential equation, as shown in an equation (18); in the advanced particle method, the heat conduction and the convection use a same set of models; for the heat conduction, the particles do not move; and for the convection, the particles move; a temperature Laplace term in the equation (18) uses the advanced particle method in the discrete form as shown in the equations (6) to (12), the thermal conductivity during a discrete process adopts a harmonic average;

[00031]$\begin{array}{cc}{\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{Trans}}=\frac{k}{\rho }{\nabla }^{2}T& \mathrm{equation}\left(18\right)\end{array}$ [0250] wherein [0251] k is the thermal conductivity, W/(m.Math.K); [0252] T is the temperature, K; [0253] calculating the radiative heat exchange according to an equation (19):

[00032]$\begin{array}{cc}{\left(\frac{\partial H}{\partial t}\right)}_{\mathrm{radiation}}=\frac{{\mathrm{\epsilon \sigma }}_s\left(T_i^4-T_j^4\right)}{\rho l_0}& \mathrm{equation}\left(19\right)\end{array}$ [0254] wherein [0255] ε is an emissivity; [0256] σ.sub.s is a Stefan-Boltzmann constant; [0257] T.sub.i is the temperature of the particle i, K; [0258] T.sub.j is a temperature of the particle j, K; [0259] determining the heat flow boundaries according to an actual situation, comprising heat source and the cooling boundaries; [0260] determining the chemical heat by chemical reactions; [0261] wherein changes of fluid particle velocities, positions and energy are calculated in the step 4, which truly reflect fluid motion and temperature field changes of the calculation object during the serious accident in an actual nuclear reactor; [0262] step 5: performing a chemical reaction calculation, wherein the embodiment only takes eutectic reactions between the UO.sub.2 core block and the Zr-4 alloy cladding into account; [0263] wherein the step 5 calculates the material component changes of a core material under the eutectic reactions which may lead to premature melting of the core block; [0264] step 6: performing a neutron physics calculation using a multigroup approximation S.sub.N difference Boltzmann transport equation as shown in an equation (20)

[00033]$\begin{array}{cc}\Omega .Math.\nabla \varphi \left(r,\Omega ,E_n\right)+{\Sigma }_t\varphi \left(r,\Omega ,E_n\right)=\int \int {\Sigma }_s\left(r,{\Omega }^{\prime },E_n^{\prime }.fwdarw.\Omega ,E_n\right)\varphi \left(r,{\Omega }^{\prime },E_n^{\prime }\right)d{\Omega }^{\prime }{\mathrm{dE}}_n^{\prime }+\frac{\chi \left(r,E_n\right)}{4\pi }\int \int v{\Sigma }_f\left(r,{\Omega }^{\prime },E_n^{\prime }\right)\varphi \left(r,{\Omega }^{\prime },E_n^{\prime }\right)d{\Omega }^{\prime }{\mathrm{dE}}_n^{\prime }+\frac{Q\left(r,E_n\right)}{4\pi }& \mathrm{equation}\left(20\right)\end{array}$ [0265] wherein [0266] Ω is a direction vector; [0267] Ω′ is another direction vector, but different from Ω; [0268] ϕ(r,Ω,E.sub.n) is a neutron angular flux density when an input is r,Ω,E.sub.n; [0269] ϕ(r,Ω′,E′.sub.n) is a neutron angular flux density when the input is r,Ω′,E′.sub.n; [0270] Σ.sub.t is a total neutron cross section; [0271] Q(r,E.sub.n) is neutron source intensity; [0272] E.sub.n is neutron energy; [0273] E′.sub.n is another neutron energy, but different from E.sub.n; [0274] Σ.sub.s(r,Ω′,E′.sub.n.fwdarw.Ω,E.sub.n) is a scattering cross section; [0275] χ(r,E.sub.n) is a fission spectrum; [0276] ν is a quantity of neutrons released per fission; [0277] Σ.sub.f(r,Ω′,E′.sub.n) is a neutron fission cross section; [0278] selecting a pressurized water reactor multi-group core cross-section database to calculate the core cross section; arranging a structured grid on a core particle geometry arrangement, and adopting a coarse mesh finite difference method to accelerate solution; using a particle mesh mapping technology to realize grid-particle information interaction; [0279] wherein the step 6 calculates the neutron angular flux density in the core, and changes an internal heat source power of the core block through the neutron angular flux density; and [0280] step 7: outputting required data; determining whether an end-of-calculation condition is met, if not, advancing the time step and returning to the step 2, if yes, ending calculation.