Acta Polytechnica Vol. 52 No. 6/2012 Hygro-Thermo-Mechanical Analysis of a Reactor Vessel Jaroslav Kruis, Tomáš Koudelka, Tomáš Krejčí Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29 Prague Corresponding author: jk@cml.fsv.cvut.cz Abstract 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. Keywords: durability of a reactor vessel, damage, coupled heat and moisture transport. 1 Introduction Many nuclear power plants across Europe are ap- proaching the end of their design durability. This raises issues connected with the high costs of decom- missioning and replacing of existing plants. Extending the serviceability of plants is therefore welcome be- cause 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 con- crete 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. 2 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 un- realistic until the tensile strength was greater than 3.7 MPa. Lower tensile strength values led to unreal- istic 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 techni- cal 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. Addi- tionally, the model required a significant number of internal variables which have to be stored. These dif- ficulties 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 67 Acta Polytechnica Vol. 52 No. 6/2012 relation σα = ( 1 −H(ε̄α)D(t)α −H(−ε̄α)D (c) α ) × [( K − 2 3 G ) ε̄v + 2Gε̄α ] , (1) where the subscript α stands for the index of the principle components of the given quantity. ε̄α de- notes the principal values of the strain tensor. The model defines two sets of damage parameters, D(t)α and D(c)α , for tension and compression, respectively. In the equation (1), the symbol H() denotes the Heav- iside 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 D(β)α = A(β) ( |ε̄(β)α |− ε̄ (β) 0 )B(β) 1 + A(β) ( |ε̄(β)α |− ε̄ (β) 0 )B(β) , (2) 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 ( 1 −D(β)α ) E|ε̄(β)α | = fβ exp ( − D (β) α h|ε̄ (β) α | w (β) cr0 ) . (3) 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]. 3 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 super- saturated 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 simpli- fications. 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. 3.1 Transport equations Künzel proposed that the moisture transport mech- anisms relevant to numerical analysis in the field of building physics are just water vapour diffusion and liquid transport [6]. Vapour diffusion is the most im- portant 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 Jv = −δp∇p = − δ µ ∇p, (4) where Jv is the water vapour flux [kg m−2 s−1], δp [kg m−1 s−1 Pa−1] 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 diffu- sion 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 Jw = −Dϕ∇ϕ, (5) where the liquid conductivity Dϕ [kg m−1 s−1] is the product of the liquid diffusivity Dw [m2 s−1] and the derivative of the water retention function Dϕ = Dw · 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) q = −λ∇T, (6) where λ [W m−1 K−1] 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. 68 Acta Polytechnica Vol. 52 No. 6/2012 3.2 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 re- sulting 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ϕ ∂ϕ ∂t = ∇T ( Dϕ∇ϕ + δp∇(ϕpsat) ) , (7)( ρC + dHw dT )∂T ∂t = = ∇T ( λ∇T ) + hv∇T ( δp∇(ϕpsat) ) , (8) where Hw [J m−3] is the enthalpy of the moisture of the material, hv [J kg−1] is the evaporation enthalpy of the water, psat [Pa] is the water vapour satura- tion pressure, ρ [kg m−3] is the mass density of the porous material, C [J kg−1 K−1] is the specific heat capacity and t [s] denotes time. The boundary of domain Ω is split into parts ΓTD, ΓϕD, ΓTN, ΓϕN, ΓTC and ΓϕC. D indicates the Dirichlet boundary conditions (prescribed values), N denotes the Neu- mann boundary conditions (prescribed fluxes) and C denotes the Cauchy/Newton boundary conditions (flux caused by the transmission). The parts ΓTD, ΓTN and ΓTC 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 = ΓTN ⋃ ΓTC and the moisture fluxes are prescribed on the part ΓJ = ΓϕN ⋃ ΓϕC. The system of equations (7) and (8) is accompanied by three types of boundary conditions. The Dirichlet boundary conditions have the form T(x, t) = T(x, t), x∈ ΓTD, (9) ϕ(x, t) = ϕ(x, t), x∈ ΓϕD. (10) The Neumann boundary conditions have the form q(x, t) = q(x, t), x∈ ΓTN, (11) J(x, t) = J(x, t), x∈ ΓϕN. (12) The Cauchy/Newton boundary conditions have the form q(x, t) = βT (T(x, t) −T∞(x, t)), x∈ ΓTC, (13) J(x, t) = βϕ(p(x, t) −p∞(x, t)), x∈ ΓϕC. (14) In the previous relationships, T(x, t) is the prescribed temperature, ϕ(x, t) is the prescribed relative humid- ity, 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 coef- ficients, T∞ is the ambient temperature and p∞ is the ambient water vapour pressure. Besides the boundary conditions, the initial conditions are prescribed, i.e. T(x, 0) = T0(x), x∈ Ω, (15) ϕ(x, 0) = ϕ0(x), x∈ Ω, (16) where T0(x) denotes the initial temperature and ϕ0(x) denotes the initial relative humidity. 3.3 Discretization of the differential equations The finite element method is used for spatial dis- cretization of the partial differential equations (7) and (8). The weighted residual statement [7] is ap- plied to the mass balance equation, assuming δT = 0 on ΓTD and δϕ = 0 on ΓϕD∫ Ω δϕ (dw dϕ ∂ϕ ∂t −∇T ( Dϕ∇ϕ + δp∇(ϕpsat) )) dΩ = 0 (17) and is also applied to the energy balance equation ∫ Ω δT (( ρC + dHw dT )∂T ∂t −∇T ( λ∇T ) −hv∇T ( δp∇(ϕpsat) )) dΩ = 0. (18) In the finite element method, temperature T and relative humidity ϕ are approximated in the form T = N(x)dT , ϕ = N(x)dϕ (19) and the gradients of temperature and relative humid- ity are also needed ∇T = B(x)dT , ∇ϕ = B(x)dϕ. (20) In the previous equations, N(x) denotes the matrix of approximation functions, B(x) is the matrix of their derivatives, dT denotes the vector of nodal tem- peratures and dϕ denotes the vector of nodal relative humidities. A set of first order differential equations is obtained in the matrix form ( Kϕϕ KϕT KTϕ KTT )( dϕ dT ) + ( Cϕϕ CϕT CTϕ CTT )( ḋϕ ḋT ) = ( Jϕ qT ) . (21) Matrices Kϕϕ, KϕT , KTϕ and KTT form the con- ductivity matrix of the problem, and they have the 69 Acta Polytechnica Vol. 52 No. 6/2012 form Kϕϕ = ∫ Ω BTDϕϕBdΩ, (22) KϕT = ∫ Ω BTDϕTBdΩ, (23) KTϕ = ∫ Ω BTDTϕBdΩ, (24) KTT = ∫ Ω BTDTTBdΩ, (25) where the conductivity matrices of material Dϕϕ, DϕT , DTϕ and DTT are diagonal matrices and the diagonal entries are equal to the appropriate conduc- tivities kϕϕ = Dw dw dϕ + δppsat, kϕT = δpϕ dpsat dT , (26) kTϕ = hvδppsat, kTT = λ + hvδpϕ dpsat dT . (27) Matrices Cϕϕ, CϕT , CTϕ and CTT form the capacity matrix of the problem, and they have the form Cϕϕ = ∫ Ω NTHϕϕNdΩ, (28) CϕT = ∫ Ω NTHϕTNdΩ, (29) CTϕ = ∫ Ω NTHTϕNdΩ, (30) CTT = ∫ Ω NTHTTNdΩ, (31) where the capacity matrices of material Hϕϕ, HϕT , HTϕ and HTT are diagonal matrices and the diagonal entries are equal to the appropriate capacities cϕϕ = dw dϕ , cϕT = 0, (32) cTϕ = 0, cTT = ρC + dHw dT . (33) Vectors Jϕ and qT contain prescribed nodal fluxes and have the form Jϕ = ∫ ΓJ NTJϕdΓ, qT = ∫ Γq NTqT dΓ, (34) where Jϕ denotes the mass boundary fluxes and qT denotes the heat boundary fluxes. 4 Computer simulation of a reactor vessel The reliability and the durability of reactor con- tainments depend directly on the prestressing sys- tem. The general results from in-situ measurements throughout the operation show the increase in de- formations and the increase in prestress losses since Figure 1: Change in tendon force gradient during service time [8] the onset of service. Most measurements also indi- cate that the temperature has a major influence on the prestress losses. These conclusions have been ob- tained, e.g., from thirty years of measured prestress at Swedish nuclear reactor containments [8]. This section presents a computer simulation of a nuclear power plant containment under cyclic temperature loading during service, when the stages of service and planned stops are changed. It is well known that increase in temperature influences the rate of con- crete creep. This fact can cause significant prestress losses of the structure. Moreover, increasing deforma- tions are observed and additional cracks may occur. An advanced two-level model is used for predicting the prestress losses and the structure response. It is a combination of a global macro-level model and a local model. The aim of the global model is to model the evolution of the prestress forces changed by the temperature and climatic loading. The local model is loaded by the mechanical and thermal load- ing from the global model. The staggered coupled thermo-mechanical analysis is the main part of the local model, which has to explain the time dependent processes in the containment wall. The heat transfer analysis runns in parallel with the mechanical analy- sis, 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 mod- els. The study presented here is a part of the over- all 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 phe- nomena at the Swedish nuclear reactor containment 70 Acta Polytechnica Vol. 52 No. 6/2012 Figure 2: Geometry — section view of the contain- ment 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 ob- served in prestress measurements at the Czech con- tainment. 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 tempera- ture will be the dominant phenomenon, while damage will have a minor effect, and the increase in radial strains will be neglected. 4.1 Basic data The containment of the nuclear power plant in the Czech Republic is a monolithic post-tensioned struc- ture made from reinforced concrete. It consists of two parts — the lower cylindrical part and the up- per 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. 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. 4.3 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 thermo- mechanical coupled algorithm of the SIFEL finite element computer code [12] was used. 4.4 Loading Temperature loading. The impact of tempera- ture is modeled by the Dirichlet boundary conditions. Temperatures from in-situ measurements (inner and outer surface) are applied directly into the computa- tion. 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 seg- ment 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 71 Acta Polytechnica Vol. 52 No. 6/2012 Figure 4: Change in the prestress force in the an- chorage system since the end of construction the anchorage system decreased by the prestress losses caused by friction in the tendon channels. • Loading prescribed directly in the tendon chan- nels which consist of radial and tangential com- ponents. The first two loadings are instantaneous. The lat- ter 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 magneto- elastic method (MEM) and they are displayed in Figure 4. It should be emphasized that several mea- surements 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. 4.4.1 Material properties and equations Coupled heat and moisture transfer was assumed be- cause 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 non- local versions of the scalar isotropic damage model, the anisotropic damage model and the orthotropic Figure 5: Isosurfaces of damage parameter Dt in the radial direction at the end of the pre-stressing phase model. 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 plastic- ity 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 con- cerned with capturing the trends in containment be- haviour, 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 experi- ments 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. 4.5 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 differ- ent 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 dif- ferent 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 de- picted in Figure 6. For the accompanying effect due 72 Acta Polytechnica Vol. 52 No. 6/2012 Figure 6: Comparison of the strain in the radial reinforcement obtained from computation (black line) and from in-situ measurements (red line – averaged data) to cracking strain evolution, the increase in radial deformation and the decrease in tendon forces dur- ing 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 measure- ments. Noticable increments in strains are caused by temperature increases after planned reactor shut- downs. 5 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 simpli- fied 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 mea- sure the tensile strenght, which is the basic property for monitoring the hypothetical damage to the con- tainment. Acknowledgements Financial support for this work was provided by project number VZ 03 CEZ MSM 6840770003 of the Ministry of Education of the Czech Republic. The financial support is gratefully acknowledged. References [1] Kruis, J., Koudelka, T., Krejčí, T.: Efficient Com- puter Implementation of Coupled Hygro-Thermo- Mechanical Analysis, Math. Comput. Simulat, 80, 2010, 1578–1588. [2] Papa, E., Taliercio, A.: Anisotropic Damage Model for the Multiaxial Static and Fatigue Be- haviour of Plain Concrete, Eng. Fract. Mech, 55, No. 2, 1996, 163–179. [3] Koudelka, T., Krejčí, T.: An Anisotropic Damage Model for Concrete in Coupled Problems. In: Pro- ceedings of the Ninth International Conference on Computational Structures Technology, (Ed- itors: B. H. V. Topping and M. Papadrakakis), Stirlingshire, UK, Civil-Comp Press, 2008, paper 157. [4] Krejčí, T., Koudelka, T., Šejnoha, J., Zeman, J.: Computer Simulation of Concrete Structures sub- ject to Cyclic Temperature Loading. In: Proceed- ings of the Twelfth International Conference on Civil, Structural and Environmental Engineer- ing Computing, (Editors: B. H. V. Topping, L. F. Costa Neves and R. C. Barros), Stirling- shire, UK, Civil-Comp Press, 2009, paper 131. [5] Koudelka, T., Krejčí, T., Šejnoha, J.: Analysis of a Nuclear Power Plant Containment. In: Pro- ceedings of the Twelfth International Conference on Civil, Structural and Environmental Engi- neering Computing, (Editors: B. H. V. Topping, L. F. Costa Neves and R. C. Barros), Stirling- shire, UK, Civil-Comp Press, 2009, paper 132. [6] Künzel, H.M., Kiessl, K.: Calculation of heat and moisture transfer in exposed building com- ponents, Int. J. Heat Mass Tran, 40, 1997, 159– 167. [7] Bittnar, Z., Šejnoha, J.: Numerical Methods in Structural Mechanics. New York: ASCE Press, 1996. [8] Anderson, P.: Thirty years of measured prestress of Swedish nuclear reactor containment, Nucl. Eng. Des, 235, 2005, 2323–2336. [9] Bažant, Z.P., Prasannan, S.: Solidification theory for concrete creep: I. Formulation, J. Eng. Mech- ASCE, 115, No. 8, 1989, 1691–1703. [10] Bažant, Z.P., Prasannan, S.: Solidification the- ory for concrete creep: II. Verification and appli- cation, J. Eng. Mech-ASCE, 115, No. 8, 1989, 1704–1725. [11] T3D — automatic mesh generator, http://mech. fsv.cvut.cz/~dr/t3d.html, [12] SIFEL — Simple Finite Elements, http://mech. fsv.cvut.cz/~sifel/index.html, 73