AP04-Bittnar1.vp 1 Introduction Assessment of the durability and serviceability of nuclear power plants is a topical problem in many countries. The most critical part of a power plant is the reactor containment, which is made from concrete with or without a steel lining. There are several requirements on the reactor containment. Basically, it must protect the reactor from external effects, and also the external environment from possible accidents. Me- chanical analysis is used for assessing the limit state, and transport analysis is used for describing the leakage of pollut- ants to the external environment. Possible accidents can lead to great pressure inside the reactor containment, which can cause damage to the concrete, and therefore there is an impact on the transport processes of the pollutants. There- fore coupled hydro-thermo-mechanical analysis is required for a correct assessment of reactor containment properties. Concrete is a heterogeneous and porous material, which leads to relatively complicated material models. The aim of the present study is to show in a condensed form the theo- retical basis of the most widely used mathematical models describing the coupled heat and moisture transport in de- forming porous media, to provide a set of governing equa- tions together with the finite element method. The theory discussed below is based on porous media theories given in [1]. 2 Mass and heat transfer in deforming porous media – a review of theory 2.1 Constitutive relations Moisture in materials can be present as moist air, water and ice or in some intermediate state as an adsorbed phase on the pore walls. Since it is in general not possible to distinguish the different aggregate states, water content w is defined as the ratio of the total moisture weight (kg/kg) to the dry weight of the material. The degree of saturation Sw is a function of capillary pressure pc and temperature T, which is determined experimentally S S p Tw w c� ( , ) . (1) The capillary pressure pc is defined as p p pc g w� � , (2) where pw> 0 is the pressure of the liquid phase (water). The pressure of the moist air, pg> 0, in the pore system is usually considered as the pressure in a perfect mixture of two ideal gases - dry air, pga, and water vapor, pgw, i.e., p p p M M TR M TRg ga gw ga a gw w g g � � � � � � � � � � � � � . (3) In this relation �ga , �gw and �g stand for the respective intrinsic phase densities, T is the absolute temperature, and R is the universal gas constant. Identity (3) defining the molar mass of the moist air, Mg, in terms of the molar masses of the individual constituents is known as Dalton’s law. The greater the capillary pressure, the smaller is the capillary radius. It is shown thermodynamically that the capillary pressure can be expressed unambiguously by the relative humidity RH using the Kelvin-Laplace law RH p p p M TR � � � � � � � gw gws c w w exp � . (4) The water vapor saturation pressure, pgws, is a function of the temperature only, and can be expressed by the Clausius-Clapeyron equation p T p T M R T T gws gws w( ) ( ) exp� � � � � � � � � � � � �0 0 1 1 , (5) where T0 is a reference temperature and hvap is the specific enthalpy of saturation. Materials having heat capacities is the term deliberately used to emphasize the similarity to the description of mois- ture retention. It is simply expressed as H H T� ( ) , (6) where H is the mass specific enthalpy (J�kg�1), T is the tem- perature (K). It is not common to write the enthalpy in an absolute way as is done here. Instead, changes of enthalpy are described in a differential way, which leads to the definition of the specific heat capacity as the slope of the H - T curve, i.e. 68 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004 Solution of Nonlinear Coupled Heat and Moisture Transport Using Finite Element Method T. Krejčí This paper deals with a numerical solution of coupled of heat and moisture transfer using the finite element method. The mathematical model consists of balance equations of mass, energy and linear momentum and of the appropriate constitutive equations. The chosen macroscopic field variables are temperature, capillary pressures, gas pressure and displacement. In contrast with pure mechanical problems, there are several difficulties which require special attention. Systems of algebraic equations arising from coupled problems are generally nonlinear, and the matrices of such systems are nonsymmetric and indefinite. The first experiences of solving complicated coupled problems are mentioned in this paper. Keywords: Coupled hydro-thermo-mechanical analysis, concrete, non-linear system of equations. C H T p p � � � � � � � � � � const. . (7) The heat capacity varies insignificantly with temperature. It is customary, however, to correct this term for the presence of the fluid phases and to introduce the effective heat capac- ity as ( )� � � �C C C Cp p p peff s s w w g g� . (8) 2.2 Transfer equations The mass averaged relative velocities, v v� s , are ex- pressed by the generalized form of Darcy’s law [1] nS k p� � � � � � � �( ) ( )v v � s r sat grad k g , (9) where � � w for the liquid phase and � � g for the gaseous phase. Dimensionless relative permeabilities kr� � 0 1, are func- tions of the degree of saturation k k Sr r w � �� ( ) (m�s�1). (10) In Equation (9), ksat (m 2) is the square (3×3) intrinsic per- meability matrix and �� is the dynamic viscosity (kg�m�1�s�1). The intrinsic mass densities �� are related to the volume averaged mass densities � � through the relation � �� � �� n S . (11) The relative permeability krw goes to zero, when water sat- uration Sw approaches Sirr, which is the limiting value of Sw as the suction stress approaches infinity [2]. The diffusive-dispersive mass flux (kg�m�2�s�1) of the wa- ter vapor (gw) in the gas (g) is the second driving mechanism. It is governed by Fick’s law J v vg gw g gw gw g g g gw g grad� � � � � � � � � � nS � � � � ( ) D , (12) where Dg (m 2 �s�1) is the effective dispersion tensor. It can be shown [1] that Jg gw g a w g g gw g g a w g g grad gra � � � � � � � � � � � � � M M M p p M M M 2 2 D D d ga g g gap p � � � � � � � � � J . (13) Here, Jg ga is the diffusive-dispersive mass flux of the dry air in the gas. Conduction of heat in the normal sense comprises radia- tion as well as convective heat transfer on a microscopic level. The generalized version of Fourier’s law is used to describe the conduction heat transfer q � � eff grad T , (14) where q is the heat flux (W�m�2), and � eff is the effective ther- mal conductivity matrix (W�m�1�K�1). The thermal conductivity increases with increasing tem- perature due to the non-linear behavior of the microscopic radiation, which depends on the difference of temperatures raised to the 4th power. The presence of water also increases the thermal conductivity. A suitable formula reflecting this ef- fect can be found in [1]. 3 Deformation of a solid skeleton 3.1 Concept of effective stress The stresses in the grains, �s, can be expressed using a standard averaging technique in terms of the stresses in the liquid phase, �w, the stresses in the gas, �g, and the effective stresses between the grains, �ef. The equivalence conditions for the internal stresses and for the total stress � lead to the expression [3]. � � � � �� nS nS nw w g g s( )1 � . (15) The assumption that the shear stress in fluids is negligible converts the latter equation into the form � �� ef sp m, (16) where � � � �� � � x y z yz zx xy T T, , , , , , { , , , , , }m 1 1 1 0 0 0 (17) and p S p S ps w w g g� . (18) Deformation of a porous skeleton associated with the grain rearrangement can be expressed using the constitutive equation written in the rate form � ( � � )� � �ef sk� D 0 . (19) The dots denote differentiation with respect to time, D Dsk sk ef� ( �, , )� � T is the tangential matrix of the porous skeleton and ��0 represents the strains that are not directly associated with stress changes (e.g., temperature effects, shrinkage, swelling, creep). It also comprises the strains of the bulk material due to changes of the pore pressure � � �p p � � � � � � � � � m s s3K , (20) where Ks is the bulk modulus of the solid material (matrix). When admitting only this effect and combining Eqs. (16), (19) and (20), we get � � � � � � � � � � �� � � �� ef s sk s sp p pm D m m� � , (21) where � � � � � � � � � � � � 1 3 3 1m I D mT m sk sk sK K K , (22) and K sk sk� m D m T 9 is the bulk modulus of the porous skeleton. For a material without any pores, K Ksk s� . For cohesive soils, K sk� K s and � �1 . The above formulas are also applicable to long-term deformation of rocks, for which � � 0 5. , and this fact strongly affects Equation (21) [4]. © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 69 Acta Polytechnica Vol. 44 No. 5 – 6/2004 Changes of the effective stress along with temperature and pore pressure changes produce a change in the solid density ��s . To derive the respective material relation for this quantity, we start from the mass conservation equation for the solid phase. In the second step we introduce the constitutive rela- tionship for the mean effective stress expressed in terms of quantities describing the deformation of the porous skeleton. After some manipulations we arrive at the searched formula ( ) � ( ) � � ( )1 1� � � � � � � � � � �n n p T � � � � � s s s s s sdiv K v , (23) where �s is the thermal expansion coefficient of the solid phase. A similar approach applied to the mass conservation equation of the liquid phase leads to the following constitutive equation � � � � � � w w w w w� � p T K , (24) where Kw is the bulk modulus of water and �w is the thermal expansion coefficient of this phase. 3.2 Set of governing equations The complete set of equations describing the coupled moisture and heat transport in deforming porous media com- prises the linear balance (equilibrium) equation formulated for a multi-phase body, the energy balance equation and the continuity equations for the liquid water and gas. Continuity equation for dry air � � � � � � � t S S k ( ) ( ) �1 1�� � � � � � � � w ga w ga ga rg sat div div u k � � g g g a w g 2 g gw g grad div grad p M M M p p � � � � � � � � � � � � � D � � � � � � 0, (25) where �u ( �u v� s ) is the velocity of a solid. Continuity equation for the water species � � � � � � � t S S k ( ) ( ) �1 1�� � � � � � � � w gw w gw gw rg sat div div u k � � g g g a w g 2 g gw g grad div grad p M M M p p � � � � � � � � � � � � � D � � � � � � � � � � � � � � � � � t S S k ( ) �w w w w w rw sat w div div (gra u k d gradg c wp p� � � � � � � � � g) . (26) Energy balance equation ( ) ( )� � � � � � C T t T C k p eff eff pw w rw sat w div grad div (gra � � � k d grad grad grad g c w pg gw rg sat g g p p C k p � � � � � � � � � � � g k ) T h t S S k � � � � � = div div (g vap w w w w w rw sat w � � � � � � � � ( ) �u k rad gradg c wp p� � � � � � � � � � � � � g) . (27) The equilibrium equation (the linear balance equation) must still be introduced to complete a set of governing equations div g w c � � �� � � � � �m g( )p S p � 0 (28) with the density of the multi-phase medium defined as � � � � � � �� � � ( )1 n nS nSs w w g g s w g. (29) Initial and boundary conditions The initial conditions specify the full fields of gas pressure, capillary or water pressure, temperature and displacement and velocities: p p p p T Tg 0 g c 0 c� � � �, , ,0 0u u , and � � ,u u� �0 0at t . (30) The boundary conditions can be imposed values on � 1 or fluxes on � 2, where the boundary � � � 1 2 , p pg g� on g 1, p pc c� on c 1, T T� on T 1 and u u� on u 1. (31) The volume averaged flux boundary conditions for wa- ter species and dry air conservation equations and energy equation to be imposed at the interface between the porous medium and the surrounding fluid are as follows ( ) ( ) � � � � � � ga ga g gw ga g gw ga w w g gw c onJ J n q J J J n � � � � � � 2 ( ) ( ) ( � � � � � gw gw gw w c w w vap eff c on grad � � � � � � q q J n 2 h T T T T� �) ,on 2 (32) where n is the unit normal vector of the surface of the porous medium, �� gw and T � are the mass concentration of water va- por and temperature in the undisturbed gas phase far away from the interface, and qga , qgw , qw and q T are the imposed air flux, the imposed vapor flux, the imposed liquid flux and the imposed heat flux, respectively. The traction boundary conditions for the displacement field are: � � �n t on u 2 , (33) where t is the imposed traction. 70 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004 4 Discretization of governing equations A weak formulation of the governing equations (25) to (28) is obtained by applying Galerkin’s method of weighted residuals. For the numerical solution, the governing equa- tions are discretized in space by means of the finite element method, yielding a non-symmetric and non-linear system of partial differential equations: K K K K C C C C uu uc c ug g u u cu cc c cg g c u p p T F u p p T � � � � � � � t t , � � � � � � � � � � � � � K K K K C C C C cu cc c cg g c c gu gc c gg g g u p p T F u p p t t , � � � � , � � � T u p p T F u p p � � � � � � � � � K K K K C C C C gu gc c gg g g g u c c g g t t t t tt t t t tt t � . T u p p T F � � � � � �K K K Ku c c g g (34) The Eqs. (34) can be rewritten in concise form as K X X C X X F X( ) ( ) � ( )� � , (35) where � �X TT � p p ug c, , , , C(X) is “the general capacity ma- trix”, K(X) is “the general conductivity matrix” and they are obtained together with F(X) by assembling the sub-matrices indicated in Eqs. (34). The dot denotes the time derivative. 5 Method of solution Coupled mechanical and transport processes after discret- ization by the finite element method are described by a system of ordinary differential equations which can be written in the form K Cr r F� � d dt (36) where K is the stiffness-conductivity matrix, C is the capacity matrix, r is the vector of nodal values and F is the vector of prescribed forces and fluxes. Numerical solution of the sys- tem of ordinary differential Eqs. (36) is based on expressions for unknown values collected in vector r at time n � 1 r r vn n nt� �� �1 � � (37) where the vector vn �� has the form v v vn n n� �� � �� � �( )1 1. (38) The vector v contains time derivatives of unknown vari- ables (time derivatives of the vector r). Eq. (36) is expressed in time n � 1 and with the help of the previously defined vectors we can find ( ) ( ( ) )C K K� � � � �� �� �� �t tn n n nv F r v1 1 1 . (39) This formulation is suitable because an explicit or implicit computational scheme can be set by parameter �. The advan- tage of the explicit algorithm is based on possible efficient so- lution of the system of equations, because parameter � can be equal to zero and capacity matrix C can be diagonal. There- fore the solution of the system is extremely easy. The disad- vantage of such a method is its conditional stability. This means that the time step must satisfy the stability condition, which usually leads to a very short time step. The previously described algorithm is valid for linear problems, and one system of linear algebraic equations must be solved in every time step. The situation is more compli- cated for nonlinear problems, where a nonlinear system of al- gebraic equations must be solved in every time step. The high complexity of the problems leads to the application of the Newton-Raphson method as the most popular method for such cases. There are several ways to apply and implement the solver of nonlinear algebraic equations. We prefer equilibrium of forces and fluxes (computed and prescribed) in the nodes. This strategy is based on the equation f fint � �ext 0 (40) where vectors f int and f ext contain internal values and pre- scribed values. Due to the nonlinear feature of the material laws used in the analysis, the Eq. (40) is not valid after compu- tation of new values from the equations summarized in Ta- ble 1. There are nonequilibriated forces and fluxes which must be suppressed. When new values in the nodes are computed, the strains and gradients of the approximated functions can be estab- lished. This is done with the help of matrices Bs and Bg where particular partial derivatives are collected. Concise expres- sions of strains and gradients are written as g B r B r� �g , � � . (41) New stresses and fluxes are obtained from the constitutive relations � �� �D q D g� , q , (42) where D� is the stiffness matrix of the material and Dq is the matrix of conductivity coefficients. The real nodal forces and the discrete fluxes on an element are computed from the relations f f qe T e e T e int int,� �� �B B� � d dg� � � � . (43) These are compared with the prescribed nodal forces and discrete fluxes and their differences create the vector of residuals R. Increments of vector v are computed from the equation ( )C K v R� ���� t n 1 (44) If matrices C and K are updated after every computation of the new increment from Eq. (44), the full Newton-Raphson © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 71 Acta Polytechnica Vol. 44 No. 5 – 6/2004 Initial vectors (defined by initialconditions) r0, v0 do until i n� (n is number of time steps) predictor ~ ( )r r vi i it� � � �1 1 � � right hand side vector y f ri i i� � �� �1 1 1K ~ matrix of the system A C K� � �� t solution of the system v yi i� � ��1 1 1A new approximation r r vi i it� �� � ~ 1 1�� Table 1 72 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004 Finite elements: 50 quadrilateral with linear approximation functions Initial conditions: T0 � 298.15 K, pc0 � 9.0e7 Pa(50 % RH), pg0 � 101325 Pa Boundary conditions: Side1 Side2 T0 – heat flux: qT � 1 K�min �1 (up to 475 K), �c � 20 W�m�2�K�1 (both sides) pc – Cauchy’s RH � 10 %, �c � 0.02 m�s�1; RH � 5 %, �c � 0.01 m�s�1 pg – Dirichlet’s pg � 101325 Pa; pg � 101325 Pa Fig. 1: Problem description 292.90 297.90 302.90 307.90 312.90 317.90 322.90 327.90 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 x (cm) T( K ) Fig. 2: Temperature profile 0.0 1000.0 2000.0 3000.0 4000.0 5000.0 6000.0 7000.0 0.0 1.0 2.0 3.0 4.0 5.0 x (cm) Fig. 3: Water vapour pressure profile method is used. If the matrices are updated only once after every time step, the modified Newton-Raphson method is used. 6 Numerical example A simple benchmark test is used to show the differences between a linear and a nonlinear solution. We have a concrete sealed specimen under high temperature conditions, Fig. 1. For this example, concrete is considered as a non-deforming medium (without a displacement field). The material model for ordinary concrete presented by Baroghel-Bouny et al. [8] extended to high temperatures is used. The results are com- pared after 1, 2 and 4 hours. 7 Conclusion Hydro-thermo-mechanical problems are extremely com- plicated due to many nonlinearities, and therefore only nu- merical methods are used for two and three-dimensional problems. A suitable method is the finite element method, where each node has 6 degrees of freedom in three-dimen- sional problems (3 displacements, temperature and 2 pore pressures). The system of linear algebraic equations (after discretization and linearization) contains many unknowns, and an appropriate solver must be used. Sequential computer codes have severe difficulties with memory and CPU time in connection with such systems. Therefore parallel computers are more often used in complicated analysis. Probably the © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 73 Acta Polytechnica Vol. 44 No. 5 – 6/2004 8.0E+07 8.5E+07 9.0E+07 9.5E+07 1.0E+08 1.1E+08 1.1E+08 0.0 1.0 2.0 3.0 4.0 5.0 x (cm) Fig. 4: Capillary pressure profile 0.42 0.44 0.46 0.48 0.50 0.52 0.0 1.0 2.0 3.0 4.0 5.0 x (cm) Fig. 5: Degree of saturation profile most powerful tool for solving of large systems of equations with the help of parallel computers is the family of domain decomposition methods. Transport processes lead to nonsymmetric indefinite sys- tems. This means that usual methods, such as LDLT decom- position, do not work for such problems, and a more general LU decomposition must be used. It seems to us that there are problems where pivoting is necessary. Acknowledgment Financial support for this work was provided by GAČR, contract No. 103/03/D145, and by the EU “Maecenas” project. References [1] Lewis R. W., Schrefler B. A.: “The finite element method in static and dynamic deformation and consolidation of porous media.” John Wiley & Sons, Chichester-Toronto, 1998, (492). [2] Fatt I., Klikoff W. A.: “Effect of fractional wettability on multiphase flow through porous media.” Note No. 2043, AIME Trans., 216, 246, 1959. [3] Bittnar Z., Šejnoha J.: “Numerical methods in structural mechanics.” ASCE Press and Thomas Telford, NY, Lon- don, 1996 (422). [4] Zienkiewicz O. C.: “Basic formulas of static and dynamic behaviour of soils and other porous media.” Institute of Numerical Methods in Engineering. University College of Swansea, 1983. [5] Krejčí T., Nový T., Sehnoutek L., Šejnoha J.: “Structure - Subsoil Interaction in view of Transport Processes in Po- rous Media.” CTU Reports, Vol. 1 (2001), No. 5. [6] Wang X., Gawin D., Schrefler B. A.: ”A parallel algo- rithm for thermo-hydro-mechanical analysis of de- forming porous media.” Computational Mechanics 19, Springer-Verlag, 1996, p. 94–104. [7] Kruis J., Krejčí T., Bittnar Z.: “Numerical Solution of Coupled Problems”, in Proceedings of the Ninth Inter- national Conference on Civil and Structural Engineer- ing Computing, B. H. V. Topping, (Editor), Civil-Comp Press, Stirling, United Kingdom, paper 24, 2003. [8] Baroghel-Bouny V., Mainguy M., Lassabatere T., Coussy O.: “Characterization and identification of equi- librium and transfer moisture properties for ordinary and high-performance cementitious materials.” Cem. And Conc. Res., Vol. 29 (1999), p.1225–1238. Ing. Tomáš Krejčí, Ph.D. phone: +420 224 354 309 e-mail: krejci@cml.fsv.cvut.cz Department of Structural Mechanics Czech Technical University in Prague Faculty of Civil Engineering Thákurova 7 166 29 Prague 6, Czech Republic 74 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No. 5 – 6/2004