Hygro-Thermo-Mechanical Analysis of a Reactor Vessel

Determining the durability of a reactor vessel requires a hygro-thermo-mechanical analysis of the vessel throughout its service life. Damage, prestress losses, distribution of heat and moisture and some other quantities are needed for a durability assessment. A coupled analysis was performed on a two-level model because of the huge demands on computer hardware. This paper deals with a hygro-thermo-mechanical analysis of a reactor vessel made of prestressed concrete with a steel inner liner. The reactor vessel is located in Temelín, Czech Republic.


Introduction
Many nuclear power plants across Europe are approaching the end of their design durability.This raises issues connected with the high costs of decommissioning and replacing of existing plants.Extending the serviceability of plants is therefore welcome because it is significantly cheaper than replacing them.However, detailed analyses of reactor vessels have to be performed, because nuclear equipment is at the centre of public attention and there are severe safety requirements.
A comprehensive analysis of existing reactor vessels contains mechanical and transport parts which are coupled together.Several analyses are necessary in order to estimate behaviour of a vessel after more than 30 years of service life.During this long period of time, the vessel has undergone through many stages that have to be modelled.
This paper deals with a hygro-thermo-mechanical analysis of a reactor vessel made of prestressed concrete with a steel inner liner.The reactor vessel is located in Temelín, Czech Republic.
Coupled hygro-thermo-mechanical analyses are very demanding, because there are many unknowns that are discretized, and therefore many degrees of freedom are needed.Moreover, the growing number of degrees of freedom in the nodes of a finite element mesh leads to larger bandwidth of the system matrix and direct solvers are significantly inefficient.At the same time, the properties of the system matrix are poor due to the very different orders of the particular entries, and iterative methods require many iterations.Efficient implementation of hygro-thermo-mechanical analysis is described in reference [1].
Due to limited space, only selected models used in the analysis are described.Section 2 is devoted to the orthotropic damage model, and Section 3 describes the Künzel model of coupled heat and moisture transport.Section 4 describes the simulation of the reactor vessel.

Mechanical models -Orthotropic damage model
In the first approach, the containment was computed with the help of a scalar isotropic damage model which showed that tensile strength plays a key role.Several analyses with different tensile strength values for concrete were performed and the results were unrealistic until the tensile strength was greater than 3.7 MPa.Lower tensile strength values led to unrealistic damage of the containment.The tensile strength was not measured experimentally, but a mean value for the compressive strength obtained from earlier experiments was set to 60 MPa.The approximate tensile strength value according to the Czech technical standard, can be set to 3.5 MPa.These results led to the conclusion that it will be necessary to use more advanced damage model for concrete which can better describe the 3D problems.
In reference [2], the authors proposed a general anisotropic model for concrete which contains nine material parameters.Laboratory measurements of the required material parameters have to be performed, but this caused difficulties in certain cases.Additionally, the model required a significant number of internal variables which have to be stored.These difficulties led to the development of a simplified version of the model.This simplified version is based on six material parameters -three for tension and three further parameters for compression.
The model is based on the following stress-strain relation where the subscript α stands for the index of the principle components of the given quantity.εα denotes the principal values of the strain tensor.The model defines two sets of damage parameters, α , for tension and compression, respectively.In the equation ( 1), the symbol H() denotes the Heaviside function, K is the bulk modulus, G is the shear modulus and εv stands for volumetric strain.
There are many evolution laws that can be used for describing D (t) α and D (c) α .In our problems, the two evolution laws for the damage parameters are used similar to the laws used in the scalar isotropic damage model.The first law gives better results for compression but it is more complicated to determine the material parameters.It can be written in the form where the superscript (β) represents indices t or c, which are used for tension and compression.A (β) , B (β) are material parameters controlling the peak value and the slope of the softening branch and ε(β) 0 is the strain threshold.Damage evolves after the strains exceed the limit value of ε0 .
The second law involves correction of the dissipated energy with respect to the size of the elements, and it provides a better description of the tension.It is defined by the non-linear equation (3), which can be solved using the Newton method In the above equation, f β represents the tensile or compressive strength, w (β) cr0 controls the initial slope of the softening branches, E is the Young modulus of elasticity, h is the characteristic size of the finite element and ε(β) α is the principal value of the strain tensor.More details about the implemented models can be found in [3], [4] and [5].

Coupled Heat and Moisture Transfer -the Künzel and Kiessl model
This model introduces two unknowns in a material: relative humidity ϕ [-] and temperature T [K].The model divides the overhygroscopic region into two subranges -the capillary water region and the supersaturated region, where different conditions for water and water vapour transport are considered.For the description of simultaneous water and water vapour transport, relative humidity ϕ is chosen as the only moisture potential for both the hygroscopic and the overhygroscopic range.This model uses certain simplifications.Nevertheless, the proposed model describes all substantial phenomena and the predicted results comply well with experimentally obtained data.This is the main advantage of the model together with easy and quick determination of the material properties measured in a laboratory.

Transport equations
Künzel proposed that the moisture transport mechanisms relevant to numerical analysis in the field of building physics are just water vapour diffusion and liquid transport [6].Vapour diffusion is the most important in large pores, whereas liquid transport takes place on pore surfaces and in small capillaries.
Vapour diffusion in porous media is described in the model by Fick's diffusion and effusion in the form where is the vapour permeability of the porous material, p denotes vapour pressure [Pa], the vapour diffusion resistance number µ is a material property and δ [kg m −1 s −1 Pa −1 ] is the vapour diffusion coefficient in the air.
The liquid transport mechanism includes the liquid flow in the absorbed layer (surface diffusion) and in the water filled capillaries (capillary transport).The driving potential in both cases is capillary pressure (matrix suction) or relative humidity ϕ.The flux of liquid water is described by where the liquid conductivity and the derivative of the water retention function D ϕ = D w • dw/dϕ, where w [kg m −3 ] is the water content of the material.The heat flux is proportional to the thermal conductivity of the moist porous material and the temperature gradient (Fourier's law) where is the thermal conductivity of the moist material.The enthalpy flows through moisture movement, and the phase transition is taken into account in the form of the source terms in the heat balance equation.

Balance equations
The heat and moisture balance equations are closely coupled because, the moisture content depends on the total enthalpy and thermal conductivity while the temperature depends on moisture flow.The resulting set of differential equations for describing the simultaneous heat and moisture transfer, expressed in terms of temperature T and relative humidity ϕ, have the form of partial differential equations defined on a domain Ω dw dϕ where H w [J m −3 ] is the enthalpy of the moisture of the material, h v [J kg −1 ] is the evaporation enthalpy of the water, p sat [Pa] is the water vapour saturation pressure, ρ [kg m −3 ] is the mass density of the porous material, The system of equations ( 7) and ( 8) is accompanied by three types of boundary conditions.The Dirichlet boundary conditions have the form The Neumann boundary conditions have the form The Cauchy/Newton boundary conditions have the form In the previous relationships, T (x, t) is the prescribed temperature, ϕ(x, t) is the prescribed relative humidity, q(x, t) is the prescribed heat flux, J (x, t) is the prescribed moisture flux, β T [W m −2 K −1 ] and β ϕ [kg m −2 s −1 Pa −1 ] are the heat and mass transfer coefficients, T ∞ is the ambient temperature and p ∞ is the ambient water vapour pressure.Besides the boundary conditions, the initial conditions are prescribed, i.e.
where T 0 (x) denotes the initial temperature and ϕ 0 (x) denotes the initial relative humidity.

Discretization of the differential equations
The finite element method is used for spatial discretization of the partial differential equations ( 7) and ( 8).The weighted residual statement [7] is applied to the mass balance equation, assuming δT = 0 on Γ T D and δϕ = 0 on Γ ϕD (17) and is also applied to the energy balance equation In the finite element method, temperature T and relative humidity ϕ are approximated in the form and the gradients of temperature and relative humidity are also needed In the previous equations, N (x) denotes the matrix of approximation functions, B(x) is the matrix of their derivatives, d T denotes the vector of nodal temperatures and d ϕ denotes the vector of nodal relative humidities.A set of first order differential equations is obtained in the matrix form Matrices K ϕϕ , K ϕT , K T ϕ and K T T form the conductivity matrix of the problem, and they have the form where the conductivity matrices of material D ϕϕ , D ϕT , D T ϕ and D T T are diagonal matrices and the diagonal entries are equal to the appropriate conductivities Matrices C ϕϕ , C ϕT , C T ϕ and C T T form the capacity matrix of the problem, and they have the form where the capacity matrices of material H ϕϕ , H ϕT , H T ϕ and H T T are diagonal matrices and the diagonal entries are equal to the appropriate capacities Vectors J ϕ and q T contain prescribed nodal fluxes and have the form where J ϕ denotes the mass boundary fluxes and q T denotes the heat boundary fluxes.

Computer simulation of a reactor vessel
The reliability and the durability of reactor containments depend directly on the prestressing system.The general results from in-situ measurements throughout the operation show the increase in deformations and the increase in prestress losses since The heat transfer analysis runns in parallel with the mechanical analysis, where the effect of temperature on concrete creep is modeled by Bazant's microprestress-solidification theory, described in [9] and [10].The local model is subsequently supplemented by suitable damage models.
The study presented here is a part of the overall reliability and durability model of nuclear power plant containment in Temelín in the Czech Republic.The computation attempts to model and explain the increase in radial deformation and the decrease in tendon forces since the onset of service.There have been many of measurements to explain these phenomena at the Swedish nuclear reactor containment with non-injected (non-bounded) prestress tendons [8] over a time period of 5 years (6.5 years in the Czech Republic).The time evolution of the tendon force in a selected tendon is plotted in logarithmic scale in Figure 1.
Two gradients of tendon force losses were also observed in prestress measurements at the Czech containment.With reference to [8], it can be concluded that an increase in temperature accelerates of creep.Any change in temperature, moisture content, and loading causes changes in the creep rate.There is no doubt that temperature is one of the sources of an increase in prestress losses.The influence of temperature will be the dominant phenomenon, while damage will have a minor effect, and the increase in radial strains will be neglected.

Basic data
The containment of the nuclear power plant in the Czech Republic is a monolithic post-tensioned structure made from reinforced concrete.It consists of two parts -the lower cylindrical part and the upper dome.The cylinder has an internal diameter of 45.00 m and the wall is 1.20 m in thickness.The dome is fixed into a massive girder.The scheme of the structure is shown in Figure 2. The leakproofness of the containment is ensured by the 8 mm thick steel lining placed inside the structure.Unbonded tendons are placed in three parallel layers in the containment wall.

Local model
A local model is presented in Figure 3.The cylindrical segment represents a periodic unit cell (PUC) from the cylindrical part of the containment with channels for prestressing tendons and with vertical, radial and horizontal reinforcement.PUC is 2.12 m in height and it covers the section of an angle of 7.5°.The prestressing tendons are not modeled.Their effect is introduced as the mechanical loading.
The finite element mesh was generated by the T3D automatic mesh generator [11].The thermomechanical coupled algorithm of the SIFEL finite element computer code [12] was used.

Loading
Temperature loading.The impact of temperature is modeled by the Dirichlet boundary conditions.Temperatures from in-situ measurements (inner and outer surface) are applied directly into the computation.The temperature cycle loading was considered in one year intervals.

Mechanical loading.
The mechanical loading of the cylindrical segment is considered as a combination of four types of loading: • Dead weight of the segment.
• Dead weight of the containment over the segment which is considered as a loading on the top surface.
• Vertical loading of the prestress forces which is also considered also as a loading on the top surface.It is computed from the reactions of • Loading prescribed directly in the tendon channels which consist of radial and tangential components.
The first two loadings are instantaneous.The latter two loadings are calculated as a multiple of the prestress forces in the tendons in the place of the anchorage system.The prestress forces values are obtained from in-situ measurements by a magnetoelastic method (MEM) and they are displayed in Figure 4.It should be emphasized that several measurements were averaged and a homogenized prestress force was used in the numerical analysis.The prestress force is not reduced due to the effect of temperature on concrete creep, and it is therefore slightly different from the force depicted in Figure 1.The data was approximated by a logarithmic regression method.In the graph, jumps which simulate in the cycle service time-the planned stop-are obtained from the global model.

Material properties and equations
Coupled heat and moisture transfer was assumed because of the climatic conditions on the exterior side.It was recognized that the relative humidity varies very little, and its influence on the containment response is negligible.The non-stationary heat transport was solved assuming constant material parameters.The mechanical part of the computation considered four types of constitutive material models, namely creep, damage, plasticity and thermal dilatation.The B3 creep model influenced by temperature and moisture changes and a damage model describe the behaviour of the concrete.Several damage models were used in the computer simulation.There were local and nonlocal versions of the scalar isotropic damage model, the anisotropic damage model and the orthotropic The application of damage models is described in detail in [5].The steel reinforcement was modeled using the bar finite elements with a plasticity model, using the Huber-Mises-Hencky condition.The thermal dilatation model was assumed in both materials (concrete and reinforcement).All material parameters were obtained from the nuclear power plant operator.Our analysis was concerned with capturing the trends in containment behaviour, and with identifying important processes in concrete that have to be modelled.The key point for future analyses is the accessibility of the material parameters.From this point of view, several experiments on concrete specimens have to be performed in order to obtain exact material parameters, e.g.tensile strength and fracture energy.Moreover, the model needs to be calibrated.

Computation Results
The relation between the response of the local model and the tensile strength of the concrete in the damage models was observed during the computer simulation.Several calculations were therefore made with different tensile strengths in order to verify the damage evolution.The scalar isotropic damage model gives the upper estimate, because the damage parameter influenced all principal directions.
More realistic anisotropic or orthotropic models should therefore be used for a reliable prognosis of the durability of containment.The distribution of the damage parameter for an orthotropic damage model is captured in Figure 5.
From the point of view of concrete creep, the different levels of the effect of temperature on concrete creep were also studied.For an explanation of the increase in radial deformation, the most monitored graphs are the strains in the radial reinforcement depicted in Figure 6.For the accompanying effect due to cracking strain evolution, the increase in radial deformation and the decrease in tendon forces during service life can be observed.The strains in the radial reinforcement are plotted in Figure 6 and are compared with the average data from in-situ measurements.Noticable increments in strains are caused by temperature increases after planned reactor shutdowns.

Conclusions
The explanation for the increase in radial strains and the decrease in tendon forces since the beginning of service is based on theoretical knowledge about concrete creep influenced by temperature changes, and also partly on measurements of prestress losses mainly performed at Swedish nuclear reactor containments.The influence of an increase in temperature during service has been proved.
The results obtained from the connecting the simplified global model and the local model show relatively good coincidence with in-situ measurements.
For the best coincidence between the computer simulation and the measurements, it is necessary to calibrate all material models and their parameters and to compare them with laboratory and in-situ measurements.In particular, it is necessary to measure the tensile strenght, which is the basic property for monitoring the hypothetical damage to the containment.
is the specific heat capacity and t [s] denotes time.The boundary of domain Ω is split into parts Γ T D , Γ ϕD , Γ T N , Γ ϕN , Γ T C and Γ ϕC .D indicates the Dirichlet boundary conditions (prescribed values), N denotes the Neumann boundary conditions (prescribed fluxes) and C denotes the Cauchy/Newton boundary conditions (flux caused by the transmission).The parts Γ T D , Γ T N and Γ T C are disjoint and their union is the whole boundary Γ.The same is valid for the parts Γ ϕD , Γ ϕN and Γ ϕC .The heat fluxes are prescribed on the part Γ q = Γ T N Γ T C and the moisture fluxes are prescribed on the part Γ J = Γ ϕN Γ ϕC .

Figure 2 :
Figure 2: Geometry -section view of the containment

Figure 3 :
Figure 3: Tendon channels and reinforcement 4.2 Global model A global model is depicted in Figure 2. It serves only for determining the prestress losses caused by temperature fluctuations.The prestress losses are determined analytically on the basis of shell theory.

Figure 4 :
Figure 4: Change in the prestress force in the anchorage system since the end of construction the anchorage system decreased by the prestress losses caused by friction in the tendon channels.

Figure 5 :
Figure 5: Isosurfaces of damage parameter D t in the radial direction at the end of the pre-stressing phasemodel.The results obtained using the orthotropic model showed the best coincidence with the in-situ measurements.The application of damage models is described in detail in[5].The steel reinforcement was modeled using the bar finite elements with a plasticity model, using the Huber-Mises-Hencky condition.The thermal dilatation model was assumed in both materials (concrete and reinforcement).All material parameters were obtained from the nuclear power plant operator.Our analysis was concerned with capturing the trends in containment behaviour, and with identifying important processes in concrete that have to be modelled.The key point for future analyses is the accessibility of the material parameters.From this point of view, several experiments on concrete specimens have to be performed in order to obtain exact material parameters, e.g.tensile strength and fracture energy.Moreover, the model needs to be calibrated.

Figure 6 :
Figure 6: Comparison of the strain in the radial reinforcement obtained from computation (black line) and from in-situ measurements (red line -averaged data)