Continuum Damage Dynamic Model Combined with Transient Elastic Equation and Heat Conduction Equation to Solve RPV Stress

: The development of the world cannot be separated from energy: the energy crisis has become a major challenge in this era, and nuclear are visualized well. In addition, we have given error estimation for h -type and p -type adaptive meshes. So, our research can provide mechanical theoretical support for nuclear energy safety applications and RPV design.


Research Motivation and Significance
Nuclear energy plays an important role in today's energy system, especially nuclear power generation. Nuclear energy is a safe, clean, and economical energy source [1]. Meanwhile, it has many advantages, such as the small size of reaction equipment, slowing

Related Work
In addition, pressure vessel is the core component of the design and operation of next-generation reactors. Fatigue damage analysis, crack propagation simulation, and pipe opening stress calculation are usually required for RPV. However, many models belong to static mechanical analysis, and the Young's modulus and Poisson's ratio of the model are calculated according to fixed values. A damage model was established to describe the dynamic changes of Young's modulus and Poisson's ratio. The reliability of RPV also includes some uncertain factors, including the existence of the coupling of internal pressure and inertial force, combined with probabilistic fracture mechanics, estimation of stress intensity factors, and in turn, these works can help to analyze the pressure vessel's fracture and reliability analysis [6,7]. The mechanical properties and electromagnetic properties are the external performance under irradiation [8]. These advantages are beneficial to establish the nondestructive evaluation technology of embrittlement. Numerical calculation combined with the local nonlinear dynamics method and this criterion based on the critical splitting stress have greatly improved the global static method to describe the crack propagation.
The structural damage dynamic model is usually used in combination with the fatigue analysis model. In continuum mechanics, the damage is calculated as a post-processing of elastic or elastoplastic macroscopic analysis. However, this important work has not been cited in the mechanical research of RPV, which also reflects the uniqueness of our work. The damage is considered to be isotropic, and there is a micro-defect closure effect on both macro and micro scales [9]. Secondly, the damage evolution equation can explain different damage mechanisms when forging alloy materials. The numerical results show that the damage evolution equation can reflect the anisotropic accelerated creep and creep fracture time under different stress levels and loading directions [10]. Through the numerical simulation of the representative volume element (RVE) of quasi-brittle materials, an anisotropic damage model with the least internal variables can be constructed [11]. The orientation distribution function of the two elastic modulus is numerically determined, and the influence of the nucleation and propagation of microcracks is considered by the phase field method. In reference [12], an energy-based damage model is proposed to simulate the crack propagation of very low cycle fatigue (VLCF). This model can be used to predict the failure period, and the comparison of fracture surfaces also shows good consistency. The above model is only applicable to the damage model of a mechanical single field independent of temperature change, which also reflects that the research work in this paper is different from the current model. We consider the influence of physical information such as temperature and shear strain on the material. The fatigue damage analysis of RPV requires a variety of theories to solve, including prediction of the crack growth of steel RPV based on the maximum main stress propagation standard and combination with the probability direction standard to predict the unstable crack path [13][14][15].
The cross-section FEM has the characteristics of fast calculation, stratification, and local mechanics, and it has been widely used to solve various engineering problems [16][17][18]. However, they are all numerical theories of the static cross-section FEM, and there are few transient cross-section FEM methods. Transient energy studies the mechanical changes in each time period, which is wider than the static practical range. This is also the reason why this paper uses the transient cross-section FEM to study the nuclear pressure vessel. Changsik provides a simple method to estimate the cross-section stress distribution of the nozzle designed according to Section 3 of the ASME code. This method requires the geometric information of pressure vessels and nozzles. The limitation of this method is that the stress distribution in the cross-section needs to use accurate stress concentration factors, and the method discussed is only effective under internal pressure [19]. The error of the RPV stress theoretical analysis method for 2D cross-section FEM analysis is relatively large at the edge. Reference [20] derived the accurate theoretical formulas of radial and axial displacement of cylindrical vessels and pipelines under thermal stress through the fourth-order differential equation. The edge effect has an important influence on the geometric deformation of pressure vessels and pipelines under thermal mechanical load. The maximum relative error of radial displacement at the edge reaches 42.2%, and the maximum relative error of axial displacement reaches 28.5%. Sectional FEM is also used to study the influence of buckling and post-buckling behavior of composite laminates [21]. The results are compared with those of two finite element models. Residual stresses have a significant influence on the buckling and post-buckling behavior of closed-section thin-walled laminated structures.
Recently, there have been some simulation models of RPV. The study of stress intensity of the RPV pipe mouth is generally controlled by parameters such as size, shape, inner radius and thickness of the nozzle, etc. It is concluded that the optimal design of the nozzle can minimize the stress intensity (Tresca yield criterion) and conflict between the quality of RPV [22]. However, the working environment of RPV belongs to multi-physics and needs to consider the interaction of temperature and stress, which is also the difference of the nozzle model established in our paper. In addition, we also compare the axial stress and hoop stress. What is more, through simulation, three-dimensional thermal hydraulic parameter distributions can be obtained; with the increase of the injection rate, the disturbance of the temperature field and the velocity field becomes more intense, and it is more likely to cause thermal fatigue [23,24].
Pressurized thermal shock (PTS) also affects the structural integrity of the pressurized water RPV. The literature [25,26] studies the pressurization-thermal shock phenomenon in pressure vessels (RPV). The results show that the assessment of crack initiation, stopping, and tearing instability in thermal shock (PTS) events (of RPV) is studied. According to the new results, the tearing process of RPV is still stable even for large initial cracks larger than the maximum assumed crack size in the code [27,28]. The most important performance factor is mainly the application of fluid and probabilistic fracture mechanics to comprehensively evaluate the structural integrity of RPV under hypothetical PTS accidents. When the emergency nuclear cooling (ECC) water is injected, a large temperature gradient will be generated, which will lead to a large thermal stress in the RPV wall. Predicting the thermo-mechanical behavior of the pressure vessel can also improve the safe coefficient, optimizing the ECCS performance [29,30]. Other literature work thermal shock force is generally in the form of a graph. In this paper, the fitting formula of thermal shock force and temperature is given, which is convenient for outputting the function value corresponding to any moment during the simulation process.
In addition, regarding the issue of stress prediction, synchronization accelerator X-ray diffraction measures the stress generated by the clad pressure vessel steel during thermal shock. Experimental measurements show that the peak stress intensity factor occurs during thermal shock, rather than in a steady state [31]. The internal structure of the pressure vessel is more complicated and accompanied by robust radiation, so the experimental measurement is relatively difficult. At present, ultrasonic technology can be used to measure the stress of the pressure vessel, and ultrasonic transducers with different frequency ranges are used to evaluate the hoop and axial residual stress [32]. The experiment indicates that it is limited, the calculated stress and ultrasonic measurement results have a high degree of consistency. In addition, the elements added in RPV steel will affect the toughness and crack resistance of the material. Figure 1 below shows the evolution of the main chemical elements of RPV in China in the last 40 years. The key properties of nuclear RPV are high strength, good toughness, corrosion resistance, good compatibility with coolant, stable microstructure, good welding, hot and cold processing performance, and developing to an ultra-high strength direction.  Furthermore, there are also many numerical methods for pressure vessel stress evaluation, such as the extended 3D finite element method (XFEM) to calculate the stress field of the reactor pressure vessel (RPV). The sub-model contains three types of cracks: axial, circumferential, and inclined directions [33]. However, there are a few coupling field simulation models with time term for RPV. In this paper, the transient thermal-mechanical coupling model is studied, which is different from other research work. We use FDM-FEM to discrete the 3D pipeline port area. In addition, there are also multi-scale coupled numerical methods for force analysis of pressure vessels [34], and simulation of thermo-hydraulic phenomena (such as heat, mass, and dissolution transmission in nuclear pressure vessels (RPV)). Four subsystems have been solved; the parameter correlation of RPV can more realistically react the heat transfer simulation of the pressure vessel. Reference [35] FVM is reliable for solving the neutron diffusion equation, and it can obtain an accurate threedimensional distribution of neutron flux and power of the core. For a metal pressure vessel, generally, corrosion-resistant materials are used to prevent the material from becoming fragile due to rust and chemical attack, and the tensile strength will be reduced. Eventually, it will cause bursting under high internal pressure. Actual nuclear pressure vessels are composite materials. There are many finite element methods for pressure vessel structures/components and pipelines [36]. They even include linear and nonlinear, static and dynamic, stress and deflection analysis, thermal problems, fracture mechanics problems, and solid coupling [37,38]. COMSOL, ANSYS, ABAQUES, etc. are commonly used in the finite element analysis software of pressure vessels and pipelines.

Contributions
The main contributions of this article are four points. Firstly, this paper proposed a simplified one-dimensional RPV estimation formulas for axial and hoop stress and introduces the working principle of RPV. The second contribution is that we proposed the continuous damage dynamics model combined with the transient cross-section FEM Method (CDDM-Fractal Fract. 2022, 6, 215 5 of 37 TCFEM), which can adapt to the variable parameter mechanical calculation model in the high-temperature environment. The change trend of Young's modulus and Poisson ratio is calculated and visualized. The third contribution is that we use axisymmetric FEM to analyze the nuclear pressure vessel, which not only improves the calculation speed but also obtains the overall mechanical change cloud map. A variable parameter model is more accurate than traditional fixed parameter calculation. In addition, we found that the stress at both ends of the RPV was significantly greater than that in the middle. The fourth contribution is that we build the physical model of the mechanical and thermal coupling, analyze the mechanical change of the RPV pipe mouth, and calculate the difference between the pipe mouth axial stress and loop stress, and finally the specific temperature and mechanical change cloud diagrams are given. In a word, our work is beneficial to the structural design and the RPV's security assessment.

Structure and Framework of This Paper
The structural arrangement and design of this paper consists of five sections. Section 2 mainly consists of two parts. Part one is a one-dimensional simplified mechanical equilibrium problem, which describes the internal pressure and axial stress, and the equilibrium problem of hoop stress. Part two mainly introduces the continuous damage dynamic model and the numerical theory of cross-section finite element (DDM-TSFEM). We convert it into a two-dimensional section to solve. Section 3 mainly introduces the axisymmetric finite element method. We obtained the three-dimensional stress-strain cloud map of the pressure vessel through thermal shock force, which is more intuitive than the section finite element method. Another feature of this model is the addition of deformation. Section 4 is mainly about the thermal field-force field coupling of the RPV pipe mouth physical model. By establishing three-dimensional transient solid heat transfer and elastic mechanics equations, the axial and radial stress variation trends are finally obtained at different times. The biggest feature of this example is that the stress change at the RPV nozzle is considered. Section 5 is mainly a summary and outlook and provides the relevant research conclusions of this paper and the problems that need to be studied subsequently.

RPV Working Principle of Nuclear Power
Nuclear power plants can convert nuclear energy into electrical energy for life and industrial use. The core component of nuclear power plants is nuclear pressure vessels. Common nuclear power plants can be divided into pressurized water reactor nuclear power plants, heavy water reactor nuclear power plants, boiling water reactor nuclear power plants, and fast reactor nuclear power plants according to different reactor principles. At present, China's main nuclear power plants are composed of pressurized water reactor nuclear power plants and heavy water reactor nuclear power plants. More than 60% of the world 's nuclear power plants are PWR nuclear power plants, which are mainly composed of reactors, steam generators, steam turbines, generators, and related system equipment.
At present, in nuclear power plants, the role of reactors is to conduct nuclear fission and convert nuclear energy into heat energy from water. Water as a coolant absorbs the heat generated by nuclear fission in the reactor, and water at high temperature and high pressure becomes saturated steam. The steam pressure promotes the rotation of the steam turbine, and the heat energy is converted into mechanical energy. Then, the steam turbine drives the generator to rotate and converts mechanical energy into electrical energy. The cooled water is pumped back to the reactor by the main pump and heated again. Thus, the cycle is repeated to form a closed cycle of heat absorption and heat release. The pressure of the loop is controlled by the regulator. Usually, the primary circuit and its auxiliary systems and plants are collectively referred to as nuclear islands (NIs). In summary, the PWR nuclear power plant converts nuclear energy into electrical energy in four steps, which are implemented by four main devices: (a) Reactor-converting nuclear energy into water heat. (b) Steam generator-transferring the heat from the high-temperature and high-pressure water in the first loop to the water in the second loop, so that it becomes saturated steam. (c) Steam turbine-converting the heat energy of saturated steam into the mechanical energy of high-speed rotation of a steam turbine rotor. (d) Generator-converting mechanical energy from the steam turbine into electrical energy. To use power generation, they need to go through multiple complex processes. The working principle diagram of the pressurized water reactor nuclear power plant is shown in Figure 2.

RPV Classification and Internal Structure
The nuclear pressure vessel is an important device of nuclear power plants. Highstrength alloy steel (Fe, Mn, C, Zn, and other elements) is generally used in the vessel. The internal material of the nuclear pressure vessel will encounter thermal impact force, high temperature, strong radiation, crack propagation, chemical corrosion, and other factors when working. It is difficult to measure the internal force by a direct experiment method. Therefore, the finite element method can be used to solve the numerical solution according to the elastic equation and boundary information. Example 1 mainly introduces the geometric two-dimensional transient elastic equation of the continuum damage dynamic model to solve the internal force of the nuclear pressure vessel. Common pressure vessels can be divided into four grades according to temperature and internal pressure. The classification results and scope standards are shown in Table 1. Table 1. Classification of nuclear pressure vessels by temperature and pressure.

Pressure Classification
Pressure Range Temperature Classification Temperature Range The geometric dimensions of common nuclear reactor pressure vessels are generally ellipsoidal spherical vessels. The inner shell of the pressure vessel is made of harder steel materials such as austenitic stainless steel. Pressure vessels are commonly used key equipment in the nuclear power plant, petrochemical, metallurgical, power generation and aerospace sectors. They generally work in high temperature, high pressure, corrosion and radiation environments. Especially nuclear pressure vessels, strong neutron radiation will also cause continuous damage to the material and cause brittle fracture. In severe cases, there is a risk of explosion. Therefore, pressure vessel design and internal load control must be strictly implemented in accordance with the regulations. Currently, nuclear power plants mainly use nuclear fission to release energy. The internal structure of the nuclear pressure vessel is shown in Figure 3.

Four Model Assumptions of RPV
Before establishing the RPV model in this article, we need to give the assumptions of each model, which will help the model to describe the scope of use more accurately.
Model 1: This model assumes that the material is isotropic, the internal pressure is uniform, and the thickness of the container wall is greater than d m ≥ 0.05 m.
Model 2: The CDDM-TCFEM method assumes that the obtained cross-sections are all continuous and uniform, isotropic materials. It satisfies the six assumptions of linear elasticity theory, including continuity, complete elasticity, uniformity, isotropy, slight deformation, and no initial stress. It is assumed that the continuous damage is a small defect, no obvious crack is formed, and the temperature will not cause the creep of the RPV vessel wall.
Model 3: The axisymmetric model assumes that the RPV shell is a symmetrical geometry, and the interior is subjected to a uniform outward pressure P.
Model 4: The thermal-mechanical coupling model of the pipe mouth assumes that the RPV material is isotropic; the process we study only emits hot steam and does not release the cooling liquid, because the release of the cooling liquid requires the addition of the hydrodynamic Navier-Stokes equation. We assume that at each moment, the temperature and pressure values remain relatively stable, and there is no sudden increase or decrease.

RPV Stress by the Simple Mechnical Balance
Model 1: Considering a simplified pressure vessel force analysis model, this vessel with radius r and wall thickness d is subjected to an internal gage pressure or thermal shock p along the longitudinal direction and hoop direction of the vessel to analyze the longitudinal stress σ l and hoop stress σ θ . This model has axisymmetric coordinates; there is no shear stress. When the working condition of the nuclear pressure vessel is stable, we make a cut across the section of the RPV to analyze the longitudinal stress σ l of the spherical pressure vessel. Although the derivation process of the problem borrows the area to represent the axial stress σ l and the surface pressure F generated by the internal pressure P, the final formula shows that the radius R and the thickness d m of the vessel are constant, and the final change of the axial stress depends only on the internal pressure of the RPV. Therefore, the axial stress obtained by the model belongs to a simplified one-dimensional approximation.
Similarly, the circumferential stress is also symmetrical. We cut the vessel along any axis. The tangential direction of the cylinder is the circumferential stress σ θ , and the equilibrium equation is established along the z direction. D is the diameter of the RPV, the pressure P acts on the projection of the half section, and S = D 2 l sin α is balanced with the circumferential stress σ θ acting on the two sections. We can finally get: The above model is only a one-dimensional static analysis on the cross section, and the calculated stress results are rough estimates. In fact, the three-dimensional force analysis of the nuclear RPV cannot be obtained by this method. We also want to get the local stress defects of the pressure vessel and the overall stress changes. When the nuclear vessel reacts, the inside is affected by thermal shock. The following sections will introduce the other two nuclear pressure vessel force analysis methods in detail. They are the cross-section method and the axisymmetric method. The hoop stress and axial stress in static equilibrium are illustrated in Figure 4. To overcome the shortcoming in deterministic approach, probabilistic approach can be adopted where the failure probability of RPV is evaluated using Monte-Carlo simulation (Sun et al., 2018). Evaluation of fracture parameter of an RPV during the PTS event using two different fracture mechanics approaches is discussed in detail in the following sections.

A simplified physical model of pressure vessel:
Consider a simplified pressure vessel force analysis model, this vessel with radius r and wall thickness d subjected to an internal gage pressure or thermal shock p. Along the longitudinal direction and hoop direction of the vessel to analyze the longitudinal stress l  and hoop stress h  .
This model is a axisymmetric coordinates, there is no shear stress.

Continuum Damage Dynamics Model
This section will introduce our proposed method continuum damage dynamics model with transient cross-section FEM in detail, which is referred to as the CDDM-TCFEM method. This model mainly assumes that the nuclear pressure vessel is an isotropic material. A large amount of heat will be released instantaneously during the nuclear fission reaction. The surrounding gas will form thermal shock, long-term erosion, hightemperature effects, and microscopic cracks formed on the inner surface of the pressure vessel. Residual stress cannot be ignored and will have a certain impact on the container itself. In order to accurately describe the magnitude of the impact, we established a continuum dynamic damage model [39,40]. First, the degree of damage of the material needs to be defined, that is, the volume fraction of the part containing microscopic defects on the surface of the material, which can be marked as w(t); then, the effective area of the . Through the stress definition, it is not difficult to obtain the effective stress expression:σ . (3) This model introduces a symmetric fourth-order damage effect tensor M, that can connect Cauchy stress σ with real stressσ. Their relationship is shown in Equation (4).
Considering the symmetry of the stress tensor, the Cauchy stress vector and real stress vector are expressed as: In the matrix form of the damage effect tensor M, the variablesû andd represent the two damage parameters. The variableû represents the damage effect of the Poissondependent transverse shear deformation. The variabled is a loss parameter related to the internal temperature of the RPV [41][42][43]. Of course, we can refer to some of the work of Lemaitre and Chaboche for related damage models. The damage effect tensor M is shown in Equation (7).
In Equation (8),û andd are two damage parameters;û represents the Poisson ratio and transverse shear strain-related damage effect in the initial state. E 0 and v 0 are the undamaged elastic modulus and Poisson's ratio, and the corresponding initial values are E 0 = 206 Gpa and v 0 = 0.3. After the material is damaged, the real elastic modulus and the Poisson ratio are nonlinear functions as follows: The internal temperature T of the nuclear pressure vessel is generally controlled at 20-300 • C. This numerical experiment simulates the temperature change of [20,600]. The change of temperature will affect some physical parameters in the material, including the Young's modulus, yield strength, thermal conductivity, and thermal expansion coefficient. In this model, the influence factor d of the damage dynamics model is modified to a function that is positively correlated with temperature. This model only studies the part of the internal temperature of the pressure vessel that is linearly increased, and it is in a periodic high-temperature state for a long time in the later period. It belongs to a nonlinear change, and the deformation will creep. This change can be described by the nonlinear relationship between strain and time:û = α ln(γT(t) + 1).
The damage caused by temperature to the container material is a nonlinear change process. The functional relationship between the nonlinear part damage factord and temperature is:d In Equation (8), where ε t is a weak white noise, then ε t obeys the standard normal distribution, which can be denoted as ε t ∼ N ξ, σ 2 , and it satisfies the relationship E(ε t ) = ξ, E ε 2 t = σ 2 ; the probability density function is shown in Equation (11): In addition, the random disturbance sequence is added, which is equivalent to a correction of Poisson' s v value, making it closer to the actual real value. The model assumes that the range of random disturbance is ∆ε t = ε t − ε t−1 = 10 −4 . The sources of uncertainty include the increase of martensite integral as well as the influence of uncertain factors such as thermal shock force and crack propagation on the material. The value of this model is to dynamically characterize the changes of Young's modulus E and Poisson's ratio v with temperature T and shear strain γ. The nuclear pressure vessel works in a high-temperature and high-pressure environment for a long time, and it is easy for the material inside the vessel to encounter thermal shock and chemical corrosion as well as the radiation of nuclear fuel and many other effects. Hence, it is essential to establish a dynamic damage model to describe this physical damage.
The traditional method considers that the elastic modulus E and Poisson's ratio v change very little or as a constant value to calculate the stress of RPV. However, in practice, these parameters change with temperature. Based on the continuous damage model, some more appropriate parameter values can be obtained from our proposed model. The influence of temperature on the material structure parameters is primarily considered. The parameters û,d, ε t represent a temperature-dependent variable. The temperature and pressure inside the nuclear pressure vessel indicate a dynamic nonlinear change trend. In the initial stage, the pressure vessel will instantly release a large amount of heat, but with the addition of the coolant system, the temperature will gradually decrease.

Transient Cross-Section FEM Method
In this section, we will introduce a two-dimensional transient cross-section FEM method. In other words, we use the CTFEM method to solve a two-dimensional linear elastic equation with finite difference approximation for the time term and finite element approximation for the space term. The pressure vessel can be divided into different sections according to the radial and axial direction. So, this problem has been simplified to many thin rings and rectangular slices: that is, turning a three-dimensional problem into a twodimensional problem [44,45]. The advantage of this method is that local details can be observed and the solution time is relatively fast. If FEM combined with ARIMA method can also be applied to the variable force prediction of RPV [46]. The time term T > 0 of linear elasticity equation is discretized by the finite difference method, and the spatial term Ω ⊂ R 2 is discretized by the finite element method. The boundary area is denoted as ∂Ω. Partial parameters of the elasticity equation need to be combined with the model of continuum damage dynamics. For the displacement of two-dimensional transient elastic mechanics, there are two degrees of freedom on each mesh node, and the displacement components along the x and y directions can be written in the form of vectors as follows Equation (12).
Here, we use displacement u to represent the two-dimensional strain tensor matrix ε(u).
In the same way, we can use variable displacement u to represent the stress tensor matrix σ(u).
The λ and µ are the Lame coefficients, and the elastic modulusÊ and the possion ratê v represent the parameters solved in the structural damage model.
The two-dimensional transient elastic equation can be simplified as: Then, the solution space of the transient elasticity equation exists in u ∈ [0, T] × Ω; differential Equation (16) can also be written in the form of a component equation.
Boundary conditions include displacement boundary and force boundary conditions; they are shown in Equations (18) and (19). The effect of the two boundary conditions can indicate the initial state of the pressure vessel, and it clearly describes the boundary force position and constraint conditions: The initial conditions corresponding to the displacement and velocity are as follows: In Equation (17),ρ(t) > 0 is the density of the pressure vessel, and the density decreases slightly with the increase of temperature [47,48]. The right end function can be denoted as f = f x , f y T , Ω × [0, T] → R 2 , and the displacement function in the initial boundary condition can be expressed as Under the initial condition, when t = 0, the corresponding displacement term function is g = (g 1 , g 2 ) T . The corresponding boundary initial velocity Then, when combined with the variational principle, the transient elastic equation is discretized. For any u = (u 1 , u 2 ), v = (v 1 , v 2 ) , and it satisfies the spatial relationship a continuous functional, and the following relationship satisfies Equation (23).
Among them, B(u, v) and a(u, v) are bilinear functions, and the specific expressions are as follows in Equation (25).
Of course, it can also be calculated directly, and the results obtained in the two forms of Equations (25) and (26) are equivalent.
The backward Euler scheme is used to discrete the time term. As for the time step ∆τ = T N , N ∈ N + , the time used in step n is t n = n∆t, and the corresponding function The initial conditions of the transient linear elastic equation are divided into displacement and velocity.
The time term is discreted by the central difference scheme, and the terms u n+1 ∈ V, 1 ≤ n ≤ N − 1 are to be solved such that we can get a semi-discrete variational equation as follows in Equation (31).
The basis function of elastic plate displacement constitutes the finite element solution . NF is expressed as the number of displacement components; then, the finite element solution of the displacement component can be written as shown in Equation (33).
We choose the text function By moving terms and sorting equations, ∀v 1h ∈ V h , i = 1, 2, we can obtain the form of the discrete function as follows in Equation (34).
Then, we bring Equation (33) into Equation (34), and we can get a displacement component Equation (35).
Bring Equation (33) into Equation (36), and we can get another displacement component, as shown in Equation (37).
Then, we can integrate Equations (35) and (36), so it is easily to obtain the vector iteration formulas. (38) Using the same method, we can get the sparse matrices A, B, and C. Finally, we transform the elastic differential equation into an algebraic iterative equation.
Further sorting out Equation (40), we can get vector X n−1 .
The initial iteration value can be obtained according to the boundary conditions, such as X 0 = The middle part and both ends of the vessel are the key positions of mechanical analysis. The nuclear pressure vessel cuts n equal parts along the longitudinal direction, and each section is actually a rectangular slice. The inner side is subjected to thermal shock, and the outer side is a free end. The upper and lower sides of the rectangle are fixed. Then, the transient linear elastic equation is discretized according to the FEM-FDM theory [49,50]. According to the material parameters, size, and boundary information of the nuclear pressure vessel, the structural mechanics problem is solved according to the principle of minimum potential energy or the variational method. In this example, the ring area and the axial direction are considered. The rectangular regions are all discretized by triangular elements of the upgraded spectrum. The two-dimensional transient elastic equations can be written as shown in Equation (42).
The displacement boundary conditions and the force boundary conditions corresponding to the numerical examples are as follows: (1) The displacement boundary conditions.
v| x=x 0 ,y=y 0 ,t = u| x=x 0 ,y=y 0 ,t = 0, (x, y, t) ∈ Γ outer . (45) Among them, Γ up and Γ down represent the upper and lower boundaries of the rectangle. d is the thickness, and H is the height of the longitudinal section of the pressure vessel.
(2) The force boundary conditions. The thermal shock force (the exterior force per unit volume) on the inside of the rectangular section is f = f x , f y T ; since the problem is a transient equation, f is a function of t, and the corresponding force boundary conditions are:

Radial Section Solved by CDDM-TCFEM Method
The radial section of the nuclear pressure vessel is a circle, the outer boundary belongs to the free end, and the inner side is subjected to the thermal shock force f , which also satisfies the two-dimensional transient elastic equation.
(1) The displacement boundary conditions. The outer side of the circular section is a fixed end, the two components of the displacement are zero, and the corresponding displacement boundary conditions are: u(x, y, t)| x=x 0 ,y=y 0 ,t =ū = 0 (x, y, t) ∈ Γ 1 .
We denoted R as the outer radius and r as the inner radius; then, the thickness of the pressure vessel can be written as d = R − r. The thermal shock force on the inside of the pressure vessel is uniformly variable force of outward extrusion. For the thermal shock force f = f x , f y T acting on the inner side of the radial cross-section of pressure vessel, the force boundary condition is as follows: p x = f x t j cos θ i ,p y = f y t j sin θ i .
For the two-dimensional transient elastic mechanics problem, the numerical solution can be solved according to the FDM-FEM theory; the time term is approximated by the second-order central difference, and the space term is discretized by a finite element. Items are discretized by triangular elements. For rectangular areas, they are finally divided into 380 domain elements, 58 edge elements, and 220 mesh vertices. As for the circular sections, they are divided into 344 elements and 256 mesh vertices. For rectangular areas and circles, all ring sections use LST elements, and the displacement field function can be expressed in the following Equation (60).
The six-node isoparametric LST element is composed of three vertices of a triangle and the midpoints of three sides, and its shape function is as follows: Among them, L 1 , L 2 are area coordinates, which can be solved according to the relationship between area coordinates and rectangular coordinates.
a i , b i , c i are replaced by the vertex coordinates of the triangular element.
Then, according to Section 4.2, the transient linear elasticity theory of this paper can be solved. This numerical experiment is mainly a dimensionality reduction processing method, which simplifies the three-dimensional problem into a two-dimensional plane problem, which makes the original problem easier to solve and improves the calculation speed. At the same time, the rectangular section and circular section use the LST quadratic element ratio. The accuracy of the approximation of the CST linear element is higher, and some preprocessing of the parameters is required before the numerical solution. The outer radius R = 3 of this numerical experiment, the inner radius r = 2.75 m, and the thickness of the pressure vessel is d = R − r = 0.25 m. The material selected for the pressure vessel is A508 metal, and its density will decrease with the increase of temperature. We fit the densityρ(t) of the material with respect to the temperature T and find that it satisfies the following functional relationship shown in Equation (66).
The fitting coefficient is α 0 = 7865, α 1 = −0.522, α 2 = 4.69 × 10 −4 , α 3 = −2.78 × 10 −7 . The Root Mean Square Error (RMSE) is E rmse = 13.64, and the goodness-of-fit R 2 = 0.991. The fitting error is very small, which also shows that the variable density function relationship established by us is reliable. In fact, the densityρ(t) of our model is a function of change with time t, the varying density values can be obtained from Table 2. Equation (67) obtains the functional relationship expression of temperature.
According to the test data of the change of temperature inside the pressure vessel with time, we found the fitting function of temperature with time t, which is in the form of a piecewise function. t+1.185 The variation trend of the internal parameters of the pressure vessel under the structural damage model is shown in Figure 5 below. Figure 5a shows the relationship between the internal temperature and thermal shock force of the pressure vessel with time. The function relationship of temperature change with time can be used. Equation (67) describes that the internal pressure of the pressure vessel is mainly generated by thermal shock, and we use a nonlinear function to approximate the change of the thermal shock f (t) inside the pressure vessel within 7250-7350s. Figure 5b contains the graph of the Young's modulus of the cladding and the parent material as a function of temperature obtained through experimental tests. According to the data trend, the Young's modulus E will decrease with the increase of temperature. Figure 5c shows (under the continuous damage model) the three-dimensional variation of the elastic modulus with the influencing variables u and d; while u and d are both temperature-related functions, the change of u has a significant effect. The numerical results show that the elastic modulus of the model with damage will decrease with the increase of u. For Figure 5d, it is the change trend of Poisson rate under the model with damage.
In traditional linear elasticity theory, Young's modulus and the Poisson rate are constant. However, in the high-temperature and high-pressure environment, these physical parameters usually change. This paper not only considers the influence of temperature but also combines the damage factor. Therefore, the improved continuous structural damage model satisfies this change rule. Figure 5d can observe that with the increase of u, the Poisson rate increases, and u has a proportional relationship with the temperature function. In other words, when the temperature range is 20-600, the Poisson rate will increase with the increase of temperature. It exists in the range of v ∈ [0.305, 0.357]; from the definition of Poisson rate, it can also show that the change value of transverse strain with the increase of temperature is greater than that of longitudinal strain.

Numerical Example 1 Result Display
When the temperature is from 20 to 600 • C, the specific heat coefficient, heat transfer coefficient, thermal expansion coefficient, and density inside the nuclear pressure vessel will change dynamically [49][50][51]. The experimental data are shown in Table 2. From the table, we can see that the coefficient of thermal expansion α increases with the increase of temperature, and the density, thermal conductivity, and specific heat capacity all increase with temperature. The experimental data are shown in Table 2. For the two-dimensional transient linear elastic equation, through the continuous damage variable parameter model established above combined with the FDM-FEM numerical theory, we can obtain the stress analysis of the two-dimensional section of the pressure vessel along the radial and axial direction [52,53]. The second-order central difference is used in the time term, and the finite element is used in the space term. The mesh generation and the application of boundary load can be referred to Figure 6a,b. The force F 1 (t) and  Figure 6c shows that when t 3 = 7254 s, the rectangular section is subjected to the horizontal right transient thermal shock force F 1 (t 3 ) = 3 × 10 6 N/m 2 , and the stress variation diagram is generated under this force. Figure 6d shows that when t 36 = 7320 s, the rectangular section is subjected to the horizontal right transient thermal shock force F 1 (t 36 ) = 1.77 × 10 7 N/m 2 , and the stress diagram is generated by the right boundary of the rectangular section. In addition, the loading strain of the ring section is different from that of the rectangular section, which is subjected to uniform radiation transient force F 2 . In the actual solution process, it needs to be decomposed F 2 into two horizontal and vertical components: that is, f 2x , f 2y = (F 2 (t) cos θ 2 , F 2 (t) sin θ 2 ). Figure 6e is the strain of the circular section under the boundary force of F 2 (t 3 ) = 3 × 10 6 N/m 2 , respectively. Similarly, when t 36 = 7320 s, the variable force F 2 (t 36 ) = 1.77 × 10 7 N/m 2 is the boundary force loaded on the inner side of the circular section and Figure 6f shows the strain ε xx of the circumferential section of RPV.  In this example, we mainly use the CMMD-TCFEM method to solve the axial and radial stress of RPV. The solution idea belongs to the dimensionality reduction method, and the three-dimensional RPV is divided into a section in the radial and axial directions for mechanical modeling. Then, through the continuous structural dynamic model with damage, we can obtain a more accurate Young modulus and Poisson rate. This numerical solution conclusion is that the elastic modulus E will decrease with the increase of temperature T, while the Poisson ratio will increase with the increase of temperature. At the same time, the fitting function of density with temperature is given in this numerical experiment. The real density corresponding to each time step can be accurately solved. Finally, through numerical comparison, it is found that the stress and strain of the pressure vessel wall material will increase with the increase of the internal thermal shock force, and the position of the boundary fixed connection has the phenomenon of stress concentration. The advantage of this numerical experiment is that it can quickly analyze how the stress of the pressure vessel changes under the transient force. This example transforms the three-dimensional problem into a two-dimensional problem. The calculation speed is improved. The defect is that only the change of local force can obtain the overall stress. In order to make up for the analysis defect, the axisymmetric finite element method is considered in Example 2, and the detailed theoretical and simulation results are in Section 5.

The Theories Axisymmetric FEM Method of RPV
Pressure vessels are similar to capsule vessels, which have the characteristics of geometric symmetry. Therefore, in addition to the section method mentioned above, the axisymmetric finite element method can also be used to solve this problem [54,55]. The stress of three-dimensional pressure vessels can be quickly obtained, since the calculation amount of the axisymmetric method is relatively small, which is more intuitive than the section method to reflect the change of internal mechanical properties of pressure vessels. The displacement function of the axisymmetric problem can be expressed as: Similar to the plane problem, for axisymmetric problems, we take one arbitrary element, and the numbers of three nodes are i, j, m, and the coordinates of nodes are respectively (r i , z i ), r j , z j , (r m , z m ). Let the corresponding node displacements be (u i , w i ), u j , w j , (u m , w m ). Bring the node coordinates and displacements into Equation (69), respectively.
The shape function matrix N and the node displacement vector {q} e can be written as Equations (71) and (72). (74) f l (r, z) = a l + b l r + c l z r (l = i, j, m).
We need to emphasize that the constitutive equation is written as σ = D(ε − ε 0 ) + σ 0 ; This formula takes into account the initial stress σ 0 and initial strain ε 0 , but our model considers the values of initial stress and initial strain is zero. At the same time, it is assumed that the RPV material in the axisymmetric model is isotropic. According to the relationship between stresses and strains, we can bring Equation (73) into {σ} = Dε; then, the element stress matrix can be written as: In Equation (76), S is a stress matrix, andÊ andû are the Young's modulus and Poisson's ratio, which can be calculated from the damage model, respectively. (77) The axisymmetric single element stiffness matrix can be obtained by the principle of virtual work.
The virtual work expression for 3D elasticity includes the volume integrals terms, which can be written as dV = rdθ(drdz) = rdθdA.
Then, the virtual strain of the element becomes Equation (80).
The mechanical equilibrium equation is established based on virtual work principle, and then, we can obtain Equation (81).
The sum of the above three vectors represents the external forces {F}, which consists of body forces, surface tractions, and point loads, respectively [56,57]. The equivalent nodal force of the triangular element is denoted as {F} and the virtual displacement is {δq e }; then, the virtual strain of the element can be written as {δε e } = B{δq} e . At the same time, the virtual displacement {δu e } = N{δq e } is eliminated on both sides, and we can obtain Equation (83).
Among them, the stiffness matrix of the triangular element is: The elements stiffness matrix satisfies the relationship shown in Equation (85).
Then, the shape function N includes the external force terms of the element. After finishing, we can get Equation (86).
Finally, each element matrix is assembled into a total stiffness matrix. Thus, we can get the final algebraic equation Kq = F.

Axisymmetric Numerical Simulation Example
This section will give examples of three-dimensional pressure vessels. The advantages of the axisymmetric method can make up for the defect of the poor overall evaluation effect of the section method, and it is more intuitive to show the stress state of three-dimensional vessels. The advantages of the axisymmetric method can make up for the defects of the poor overall evaluation effect of the section method, and it is more intuitive to show the stress state of three-dimensional vessels. The axisymmetric problem has distinct characteristics: the geometric structure and the boundary constraint condition are symmetrical about the central axis. The axisymmetric problem has different characteristics: the geometric structure and the boundary constraint condition are symmetrical about the central axis.
The thickness d = 25 mm of the wall of the nuclear pressure vessel in this numerical experiment, and Young's modulus of continuous structural damage model isÊ. Similarly, Poisson's ratio of the continuous structural damage model isv. The material density of the pressure vessel changes with the increase in temperature. The specific function expression is shown in Equations (67) and (68), and the variable density function is denoted asρ(t). As for the geometric parameter of the nuclear pressure vessel, the simplified nuclear pressure vessel can be considered a combination of a hollow cylinder and a semi-ellipsoid. The geometric parameters and dimensions are as follows: the thickness t = 0.25 m of the vessel wall, the radius of the bottom of the vessel is R A 1 = D 1 2 = 5m, the radius of the inner wall is R A 2 = D 2 2 = 4.75m, and the dome height formula with a curved radian that satisfies the following relation is: Among them, R k = 0.1D 1 = 1m, R c = 0.9D 2 = 8.55m, the maximum angle between the tangent of the dome and the vertical direction of the container wall is α.
After calculation α = 0.519 rad, the height of the container wall is H A = 10 m. For a pressure vessel, the axisymmetric method is mainly used to solve the problem. The displacement boundary conditions can also be called boundary constraints. u(r, θ, z)| r=r 0 ,θ=θ 0 ,z=z 0 =ū = 0 (89) v(r, θ, z)| r=r 0 ,θ=θ 0 ,z=z 0 =v = 0 (90) w(r, θ, z)| r=r 0 ,θ=θ 0 ,z=z 0 =w = 0 (91) The impact force generated during the reaction of the pressure vessel is assumed to be uniformly acting on the inner wall of the vessel, and this mode of action is generated along the normal direction of the vessel wall surface. The forces received by the three-dimensional nuclear pressure vessel are all face forces, which are different in actual applications. The impact force on the inner wall of the nuclear pressure vessel is different. At a certain moment, the impact force on the inner wall of the nuclear pressure vessel by gas is a constant force P. The surface force is the force on the surface of the object. Motion is an internal force, and only boundary elements may have surface forces. The two components of surface forceP = {P r , P z }; if the surface force on the edge of the element is q, the equivalent nodal load of the element is: In this example, the face force is linearly distributed perpendicular to the surface of the object. This face force is very common, such as dams hitting a flood, the air flow of the aircraft engine hitting the outer wall, and the wheel pressing the ground. The effect placed on the edge of the unit ij is perpendicular to the linearly distributed surface force on the surface of the object, the force at node i is q i , the force at node j is q j , and the surface force at any point on the ij side of the element is decomposed into hoop and axial components.
Then, we determine the linearly distributed surface force perpendicular to the surface of the object, and the equivalent nodal load of node i is: The equivalent nodal load of the linearly distributed surface force node j perpendicular to the surface of the object is: Finally, the axisymmetric stress tensor can be calculated according to the following equation: The axisymmetric finite element method is used to solve the pressure vessel. The solution steps are carried out according to the following points: material parameter definition, geometric region construction, mesh division, and physical field selection. Then, load conditions and fixed constraints are added. Finally, it is transformed into algebraic equations and output stress cloud diagrams. The geometric cross-section is divided into elements on the rz plane, and quadrilateral elements are used for distillation. The discretized area contains 575 domain elements and 250 boundary elements, the number of degrees of freedom is 2541, the total time is t = 10 s, and the symmetric axis is r = 0. As for the fixed constraint condition of the vessel, the outer boundary displacement is 0, and the load boundary load conditions are that the inner side is subjected to uniform outward pressure; the range is fromP = 3 × 10 6 pa toP = 1.68 × 10 7 pa. Figure 7a shows the stress and deformation of the 2D symmetry plane of the pressure vessel. The output result contains two physical quantities: stress and deformation. The deformation is mainly the displacement change, and its amplification factor is α = 128.6. The purpose is to observe the largest displacement change position more clearly. Similarly, the greater value of stress is also concentrated in the middle part. Figure 7b shows the radial strain ε rr cloud diagram of the pressure vessel at F = 13 Mpa, which can also reflect that the 2D section stress results are consistent with the 3D results. The center of the cylinder is larger, and the surrounding forces are relatively large, stable, and uniform. Figure 7c shows the von Mises stress figure of the pressure vessel at F = 3 Mpa. Figure 7d shows the von Mises stress of the pressure vessel at F = 16.8 Mpa. The solution process is also to first calculate the single stiffness matrix, including the application of boundary conditions, and then assemble the total stiffness matrix to form a large sparse matrix and finally solve the linear system.

Three-Dimensional (3D) Multi-Physics Field Model of RPV
The working environment of the real pressure vessel is relatively complex. After nuclear fission, a large amount of heat is released in a short time. Generally speaking, there are four exhaust pipes on the wall of the pressure vessel, two pipes are the inlet and outlet of the coolant, and the other two pipes are the hot steam inlet and outlet. In the initial state, the temperature inside the pressure vessel is very high. With the addition of a cooling system (ECC), the temperature gradually drops from 350 to about 100 • C. The first two examples only consider pressure related to the force condition of the container, the established geometric model does not consider the temperature change around the exhaust pipe and the state of stress distribution, and the stress change of the pipe port is more important for the sealing of the entire container and the pressure distribution. Therefore, it is very important that it is necessary to study the state of stress and temperature changes at the pipe mouth of RPV. Example 3 established an coupled model of three-dimensional temperature field and stress field [58][59][60]. The FDM-FEM numerical method is used to solve the problem, and the output is the corresponding temperature and stress.

Three-Dimensional Transient Elastic Equation
In fact, the three-dimensional nuclear pressure vessel internal force analysis is a rather complicated analysis process. Due to some uncertainty of the impact force inside the vessel, the geometric structure is not completely symmetrical, and there are multiple pipes near the top of the nuclear pressure vessel. These interfaces are mainly responsible for the discharge of exhaust steam, the input of coolant, etc. The stress change analysis near the interface has always been the focus of scholars. Due to the relatively high surface temperature of the interface, the final force analysis is not simply a mechanical problem. It is more scientific to use the knowledge of multi-physics coupling to solve the problem. This example is mainly responsible for the coupled modeling of the thermal field and the force field, and it analyzes the force situation near the pipe mouth of the pressure vessel.
The solution of nuclear pressure vessels in the thermal-mechanical coupling field contains two important equations: namely, the convective-diffusion equation and the equilibrium equation of solid mechanics. When the reaction of the pressure vessel is stable, the overall internal temperature tends to be balanced, but the local temperature changes are quite different due to the action of the coolant and the position very close to the core. Since the nuclear reactor reaction is a continuous process and is closely related to time, the analysis of temperature change should consider the transient 3D heat conduction equation and the 3D transient mechanical equilibrium equation. The following will give the 3D transient heat conduction equation.
First, we establish the force field balance equation of the nuclear pressure vessel under the action of thermal shock. Among them, u = (u, v, w) T is the displacement field function. Since this model belongs to the multi-physics coupling model, the RPV material ρ(t) is a nonlinear function of time t, which will decrease with the increase of temperature, σ is the three-dimensional stress tensor matrix, T is the corresponding thermal shock force. The three-dimensional gradient operator can be expressed as ∇ = ∂ ∂x , ∂ ∂y , ∂ ∂z ; then, the three-dimensional transient mechanical equation is expressed as shown in Equaiton (99): The Cauchy stress tensor introduced by a three-dimensional deformable solid is expressed as: For the 3D stress tensor and strain, the relationship satisfies Equation (101).
where E is the Young's modulus of the elastomer, u is the Poisson ratio, and the Lamé constant formula is: The relationship between the three-dimensional stress tensor and strain is σ = Dε, where D is the elastic matrix, and the final stress can be expressed in the form of displacement: Then, by substituting Equation (103) into Equation (104), we can convert the threedimensional linear elastic stress equation into three displacement classification equations. Then, using the linear weighted form of the displacement basis function in place of the displacement in the equation, the weak form of the Galerkin finite element is obtained. The scalar equations along the x, y, and z directions are: For the boundary conditions of the exhaust pipe of the RPV, there are two cases of displacement boundary conditions and force boundary conditions. The displacement boundary of the outer side of the pipe and the outer wall of the RPV can be assumed to have an initial value of 0, which can be expressed as: The wall of the RPV container is mainly subjected to radial impact force, and there is also a small shear force in the z direction. The inner side of the pipe wall will also be subjected to thermal shock force. If there is cooling liquid inside the pipe, gravity must also be considered, which belongs to fluid-solid coupled heating conduction model. This example only considers the force and temperature changes at the nozzle of the RPV exhaust steam, which satisfies the following mechanical boundary conditions:    σ xx n x + σ xy n y + σ xy n z = f x (t) σ yx n x + σ yy n y + σ yz n z = f y (t) (x, y, z, t) ∈ Γ RPV × [0, T m ] σ zx n x + σ zy n y + σ zz n z = f z (t) As for Equation (106), the time term is discreted by the central difference scheme, the terms u n+1 ∈ V, 1 ≤ n ≤ N − 1 are to be solved such that we can get a semi-discrete variational equation as follows in Equation (107).
The basis function of elastic plate displacement constitutes the finite element solution space u n h = u n 1h , u n 2h , u n 3h T , which satisfies the relationship u 1h , u 2h , u 3h ∈ U h = span ϕ j NF j=1 . NF is expressed as the number of displacement components. The discrete form of the space term of the 3D linear elastic equation is similar to the two-dimensional term, and there are two differences. The first point, The 3D linear elastic equation has one more equation about the displacement w component than the 2D linear elastic equation when the equation is discretized. The second point, When the space term of the 3D linear elastic equation is discretized by FEM, the basis function selected is different from that of the 2D. The 3D lowest-order basis function is a tetrahedral element, while the 2D discretization is a linear CST element. Then, according to the variational principle, the variational form is obtained, and then the basis function is brought in to obtain the Galerkin weak form. This continuous differential equation is transformed into a discrete algebraic equation. Finally, we can obtain a matrix iterative equation, and each iteration needs to solve a linear equation.

Three-Dimensional Stress Analysis Model of RPV
The 3D elastic force analysis model is similar to the 2D model. For spatial dispersion of the 3D pressure vessel, we can use the tetrahedral element to discretize it. When the asymmetric 3D structure is discrete, the axisymmetric method is not feasible, and a Cartesian coordinate system needs to be established, using the general finite element discrete model, the displacement linear function of the tetrahedral element can be expressed as the following form: Then, we combined with the four nodal coordinates of the tetrahedral element, and we use the shape function to represent the element displacement The matrix form of the element displacement can be expressed as {d e } = u v w T = Nq e , and the element stress matrix is The value of the matrix S is mainly related to the tetrahedral node coordinates In Equation (111), the following relational expression is satisfied Similarly, the expression form of the stiffness matrix can be obtained according to the principle of virtual work: The total equivalent nodal load array of the element due to body force, surface force, and concentrated force is: After obtaining the stiffness matrix of the element, it is necessary to synthesize the total stiffness matrix, the stiffness matrix of the N e tetrahedral elements on the three-dimensional solution area, and the node load according to the element node coding rules from the total stiffness matrix. Then, finally, we get a large sparse system of linear equations Kq = F.
The above is the stress solution of 3D RPV from the perspective of a virtual work principle. It is convenient to solve the steady-state problem, for the transient problem, we need to change the force F(t) of the integral term in each calculation. In Section 4.3.2, we will introduce another method to solve the stress at different times according to the transient elastic equation.

Three-Dimensional Transient Heat Conduction Equation
The temperature variable T(x, y, z, t) is a multivariate function of coordinates and time. When T(x, y, z, t), it indicates that Q does not change with time, and it indicates that the temperature of the thermally conductive object does not change with time after heat exchange. This process is called the steady-state temperature field. When ∂T ∂t = 0 is the transient temperature field, the difference between the transient temperature field and the steady-state temperature field is time variable t. According to the Fourier heat transfer law and the energy conservation law, the energy balance differential equation in the rectangular coordinate system satisfies the following relationship: The heat transfer equation [61][62][63] satisfied by the heat transfer inside the RVP material, its equation, and related parameters are as follows: For Equation (116), ρ is the material density, the unit is kg/m 3 , and c is the specific heat capacity at constant pressure of the RPV material J /(kg · k). The internal nuclear reaction of RPV produces enormous heat; Q is the heat generated by the internal heat source. k = k x , k y , k z is the thermal conductivity vector along different directions x, y, z. Furthermore, ∂T ∂x , ∂T ∂y , ∂T ∂z respectively represent the heat that flows in the x, y, z direction, u c = u x , u y , u z is the convection velocity terms. Furthermore, we can get that the equivalent form of Equation (116) is the differential Equation (117).
In addition, the temperature field distribution in the solution domain Ω needs to meet certain boundary conditions.
(1) Class I boundary conditions: The solid surface temperature is a known function of the time t. T 1 (x, y, z, t) =T(x, y, z, t), T 1 (x, y, z, t) ∈ Γ (118) (2) Class II boundary conditions: The thermal flow density of the solid surface is equal to the change value of the temperature T in the direction of each component.
(3) Class III boundary conditions: The difference between the heat flow density of the solid surface is proportional to the surface temperature T and the fluid surface temperature T c .
n x , n y , n z is the direction cosine of the normal line outside the boundary,T(x, y, z, t) is a given temperature, ∇T 3 is the heat flow density vector on the boundary Γ 3 , h is the thermal conductivity coefficient W/ m 2 · K on the boundary, T a is the insulating temperature of the boundary layer under natural convection conditions, and the combination of all boundaries can be expressed as Then, it is assumed that the three-dimensional RPV is a homogeneous material; that is, the thermal conductivity along different directions is the same [64], and there is k x = k y = k z . So, the convection-diffusion equation can also be written as (117), which is equivalent to Equation (121).

∂T ∂t
Perform Galerkin integration, multiply both sides of Equation (120) by the test function ϕ(x, y, z), obtain the corresponding discrete equation according to the variational principle, and consider the right-hand term of Equation (120).
Using the Gauss divergence theorem to further simplify the equation, we can get After considering the left-hand term of Equation (124), the discrete form of Galerkin can be finally obtained:T The tetrahedral element (4-NQ) is used to discrete the RPV, and the discrete equation in the whole region is obtained.
Then, we abbreviate the diffusion matrix of Equation (126) as: Meanwhile, M ij of Equation (127) is called the mass matrix The matrix C ij is the convection matrix: Finally, the partial differential Equation (129) is transformed into a system of ordinary differential equations, and the specific form is as follows: Neumann boundary conditions are imposed on the right-hand term components of Equation (130); its components are of the form: The boundary temperature of the outer surface of the RPV belongs to the Drichlet boundary condition, and the substitution method can be used. If there is a known function , the i − th equation of the overall discrete Equation (129) can be replaced by the Drichlet boundary function T i and the corresponding mass matrix M. The i − th row diagonal element of the convection matrix C and matrix B is i.

h-p Method Error Estimate
(1) p-type adaptive error analysis The finite element cluster is denoted as {e, p e , Σ e }, the continuous function space can be written as Ω, and the adaptive mesh discretization of the region is denoted as T h ; for any element in the region T h , it satisfies ∀e ∈ T h , h e → 0, h e > 0, h e p e ≤ const [65,66]. If π e is a higher-order approximation operator on the element, π h is a higher-order approximation operator on the overall region, then there is a constant C such that the following interpolation error estimation holds: When p = q, H k+1,p (Ω) and H m,p (Ω) are two Sobolev spaces, and we have the following relation established: Among them, the constant C is related to m, n and the reference elementê, σ n is the unit sphere volume in the space R n ; then, the C( T h ) range satisfies the following expression: (2) h-type adaptive error.
h -type finite element, where h represents the maximum size of the element. In the calculation process, the method does not change the type of element but improves the calculation results by continuously reducing the geometric size of the element, that is, refinement mesh. Because the order of elements in this method is generally low, it is also called low-order finite element method.
If v ∈ H k+1 (Ω) V, ∀e ∈ T h , h e → 0, h e > 0, h e p e ≤ const. Then, the following estimation formula is established.
Among them,h is the average side length of the cells of the adaptive mesh, and N A represents the total number of cells in the adaptive mesh region Ω.
When the order of the basis function is p-order, the numerical error of adaptive FEM generally has a relationship with the quality of the adaptive mesh [67]. If the size of the element side length of the adaptive mesh varies from the edge to the center point O where the mesh size is the smallest, and decreases exponentially, the mathematical formula is written as: Then, the newly formed adaptive grid error also has the characteristics of exponential change. The specific error estimation expression is:

Numerical Result Display
The RPV parameters are divided into geometric parameters and basic material property parameters. For the RPV radius R = 5 m we simulated, the wall thickness d = 15 cm, the radius of the pipe mouth r pip = 10 cm, the thickness of the pipe wall d pip = 4 cm, and the height of the RPV calculation area H = 0.6 m. This example uses tetrahedral elements, and we perform discretization. Figure 8 below shows the mesh quality evaluation diagram. This example simulates an RPV with four nozzles. The structure contains 58,910 domain elements, 30,345 boundary elements, and 3985 edge elements. The green mesh in Figure 8a indicates that the mesh quality is relatively good, and the yellow part belongs to the part with poor mesh quality. The overall quality of the mesh divided in this example is relatively good. Due to the symmetry of RPV, we only need to study the stress and temperature changes of a nozzle. Figure 8b is the result of meshing the tetrahedral element of one nozzle. The vicinity of the nozzle is the focus of our research, so the adaptive meshing method is adopted. It has a total of 9735 domain elements, 3297 boundary elements, and 471 edge elements, and the mesh has the characteristics of self-adaptation, which can not only ensure the advantages of fast calculation speed but also ensure the accurate description of the details of stress changes at the interface of the pipe wall. It can smoothly transition the numerical results of positions with large stress gradient changes, which makes the visualization effect better. The mesh quality evaluation chart shows that the mesh effect near the exhaust pipe is relatively poor, and the geometric change of the RPV container wall is gentle, so the mesh quality of the element is better.
Regarding other numerical results, Figure 8c shows the inner temperature variation diagram of the RPV nozzle at the calculation time t = 420 s, which reflects that the temperature of the nozzle is relatively lower than that of other positions. The reason is that the coolant added at the nozzle can reduce the temperature near the nozzle. Figure 8d is the isothermal line diagram of the three-dimensional nozzle. The outer side of the RPV is shown in the figure, and different colors represent different temperatures. Similarly, it can be seen from the figure that the temperature near the nozzle is relatively low, and there is a cooling effect caused by the combined action of air cooling and coolant. Figure 8e is the displacement cloud map of an RPV nozzle at t = 7200 s, which reflects the initial starting time; the force generated by the internal pressure is large, and the temporary change of the cooling system is not very obvious. Figure 8f is the strain ε xx of RPV, and the change range near the nozzle is significantly larger than the change of the internal value of the nozzle. The change of strain and displacement is mainly affected by the thickness of the material, the external load, the elastic modulus of the material, the Poisson ratio, and other basic properties. The results show a comprehensive response; not only the embodiment of the load change but also the change of the displacement is significantly larger than the strain. This can also be reflected from the geometric constitutive equation.
The symmetrical geometric appearance is helpful for element calculation, and the numerical results show that the stress change is relatively large at the intersection of the pipe of the RPV. In the actual application process, the RPV is accompanied by multiple exhaust pipes, which can exhaust gas, add coolant, etc. Therefore, the third numerical model of this paper is mainly to analyze the stress change near the pipe mouth, and we have established a three-dimensional finite element model of thermal-mechanical coupling. The adaptive mesh is used to discretize the pipe mouth area. Finally, the numerical solution cloud map of the stress, strain, and displacement of the RPV pipe mouth is obtained. Then, we compared the radial and hoop stress changes of RPV at different times. The physical model established in this paper and the new numerical method proposed in this paper have important reference values for the stress analysis of RPV, and the method can also be transferred to other coupled physical models. The safe control of nuclear energy production is a meaningful research topic. boundary elements and 471 boundary elements. The grid has the characteristics of self-adaptive. This can not only ensure the advantages of fast calculation speed, but also ensure the accurate description of the stress change details at the interface of the pipe wall. It can smooth the numerical results of the position where the transition stress gradient changes greatly, and make the visualization effect better. The grid quality assessment map shows that the grid effect near the exhaust pipe is relatively poor, and the geometric change of the RPV container wall is gentle, so the grid quality of the unit is better. (d) (c) (b) (a) gradient changes greatly, and make the visualization effect better. The grid quality assessment map shows that the grid effect near the exhaust pipe is relatively poor, and the geometric change of the RPV container wall is gentle, so the grid quality of the unit is better.  Figure 9a is the circumferential stress of the RPV vessel wall. We choose four time points: t = 1200 s, t = 2400 s, t = 3600 s, and t = 7200 s. From the curve trend, it can be seen in the initial t = 1200 s and subsequent t = 7200 s, the hoop stress is larger, mainly because the internal pressure is larger. Of course, the thermal stress formed by temperature will also have a partial influence. It should be noted that the distance from the outer surface of the RPV (d = 150 mm) will produce reverse stress on the vessel wall, which is formed by the interaction between internal stresses. Similarly, Figure 9a,b form a set of stress comparison diagrams. Figure 9b is the axial stress above the pipe mouth. The trend is basically similar to that in Figure 9a, and the stress gradually decreases. The difference is that the stress near the pipe mouth changes obviously, the gradient value is relatively large, and the other positions decrease slowly, which is mainly due to the stress concentration at the pipe mouth. This also shows that the numerical results are in good agreement with the actual situation.

Conclusions
The purpose of this paper is to study the stress variation of nuclear pressure vessels. The mechanical models of pressure vessels established in this paper are from onedimensional to three-dimensional. The theory of each model and three numerical examples are given. For a one-dimensional model, the equilibrium equation is mainly established according to the internal pressure, axial stress, and circumferential stress of RPV. The model belongs to the rough estimation of stress, and the error is too large. The conclusion of the theoretical model is that the circumferential stress is twice the axial stress. In this paper, the continuous damage dynamic model with cross-section finite element method (CMMD-TCFEM) is proposed. The advantage of this model is that it can dynamically describe the change of physical parameters of RPV under the action of the loss factor. The model has the characteristics of fast calculation, layered, and localized display. The model can obtain the stress distribution cloud map of axial and radial 2D sections. The numerical conclusion is that the Poisson's ratio increases with the increase of temperature (û,E t ), and the Young's modulus decreases with the increase of temperature (û,d).
In the second numerical example, axisymmetric theory is mainly given. After giving the geometric parameters and material parameters, the model can output the stress and strain nephogram of RPV under different internal pressures. The other contribution is to compare the initial state of RPV and the deformation of RPV after internal pressure. Compared with the CMMD-TCFEM method, the 3D axisymmetric method can obtain the overall stress and strain cloud map, which is more complete and intuitive. The numerical results show that the stress change at both ends of RPV is significantly greater than that of the middle vessel wall, and the reinforcement method should be adopted at both ends of the RPV. In practical application, RPV is accompanied by multiple exhaust pipes to release steam and add coolant, and stress modeling and solution near the exhaust pipe are also very important. Therefore, in example 3, a three-dimensional finite element model of thermal-mechanical coupling is established, and the adaptive mesh is used to discrete the solution area. Finally, the stress, strain, displacement, and other information of the RPV pipe mouth are obtained. We also compare the changes of radial stress and circumferential stress near the RPV pipe mouth. The numerical results show that the distance d from the inner surface is from inside to outside, and the radial stress and circumferential stress near the RPV pipe mouth gradually decrease.
In short, the new numerical method proposed in this paper corresponds to a number of different models, from one-dimensional to three-dimensional, from a single physical field to a multi-physical field. The continuous damage dynamic model is successfully combined with the finite element method, which can better characterize the elastic modulus and Poisson ratio with the change of the damage factor. The axisymmetric model and the coupling model of the nozzle can provide theoretical and simulation experience for RPV design and stress simulation. Our method reflects the multidisciplinary intersection and can solve the same problem from different angles. Our future work will focus on the intelligent application of nuclear energy. It includes optimizing fuel metering and periodically automatically detecting the performance of RPV materials. Machine learning is used to predict the internal temperature of RPV in real time, and multi-scale theory is used to analyze the defects and crack propagation of RPV materials. The application of these problems will contribute to the effective improvement of nuclear technology and provide a good theoretical basis for the application and design of a new generation of nuclear energy, so that nuclear energy can better benefit mankind and create more energy value for society.