Solution of Nonlinear Coupled Heat and Moisture Transport Using Finite Element Method

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.


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.Mechanical analysis is used for assessing the limit state, and transport analysis is used for describing the leakage of pollutants 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.Therefore 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 theoretical basis of the most widely used mathematical models describing the coupled heat and moisture transport in deforming porous media, to provide a set of governing equations together with the finite element method.The theory discussed below is based on porous media theories given in [1].

Mass and heat transfer in deforming
porous media -a review of theory

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 S w is a function of capillary pressure p c and temperature T, which is determined experimentally The capillary pressure p c is defined as where p w >0 is the pressure of the liquid phase (water).
The pressure of the moist air, p g >0, in the pore system is usually considered as the pressure in a perfect mixture of two ideal gases -dry air, p ga , and water vapor, p gw , i.e., In this relation r ga , r gw and r 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, M g , 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 The water vapor saturation pressure, p gws , is a function of the temperature only, and can be expressed by the Clausius-Clapeyron equation where T 0 is a reference temperature and h vap is the specific enthalpy of saturation.Materials having heat capacities is the term deliberately used to emphasize the similarity to the description of moisture retention.It is simply expressed as where H is the mass specific enthalpy (J×kg -1 ), T is the temperature (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.

C H T
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 capacity as

Transfer equations
The mass averaged relative velocities, v v as , are expressed by the generalized form of Darcy's law [1] where a = w for the liquid phase and a = g for the gaseous phase.
Dimensionless relative permeabilities k ra Î 0 1 , are functions of the degree of saturation In Equation ( 9), k sat (m 2 ) is the square (3×3) intrinsic permeability matrix and m a is the dynamic viscosity (kg×m -1 ×s -1 ).The intrinsic mass densities r a are related to the volume averaged mass densities r a through the relation r r The relative permeability k rw goes to zero, when water saturation S w approaches S irr , which is the limiting value of S w as the suction stress approaches infinity [2].
The diffusive-dispersive mass flux (kg×m -2 ×s -1 ) of the water vapor (gw) in the gas (g) is the second driving mechanism.

It is governed by
Here, J g ga is the diffusive-dispersive mass flux of the dry air in the gas.Conduction of heat in the normal sense comprises radiation 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 where q is the heat flux (W×m -2 ), and c eff is the effective thermal conductivity matrix (W×m -1 ×K -1 ).
The thermal conductivity increases with increasing temperature due to the non-linear behavior of the microscopic radiation, which depends on the difference of temperatures raised to the 4 th power.The presence of water also increases the thermal conductivity.A suitable formula reflecting this effect can be found in [1].

Concept of effective stress
The stresses in the grains, s s , can be expressed using a standard averaging technique in terms of the stresses in the liquid phase, s w , the stresses in the gas, s g , and the effective stresses between the grains, s ef .The equivalence conditions for the internal stresses and for the total stress s lead to the expression [3].
The assumption that the shear stress in fluids is negligible converts the latter equation into the form where e 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 where K s is the bulk modulus of the solid material (matrix).
When admitting only this effect and combining Eqs. ( 16), ( 19) and (20), we get where and is the bulk modulus of the porous skeleton.For a material without any pores, K K sk s = .For cohesive soils, K sk K s and a =1 .The above formulas are also applicable to long-term deformation of rocks, for which a £ 05 ., and this fact strongly affects Equation (21) [4].
Changes of the effective stress along with temperature and pore pressure changes produce a change in the solid density & r 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 relationship 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 where b 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 where K w is the bulk modulus of water and b w is the thermal expansion coefficient of this phase.

Set of governing equations
The complete set of equations describing the coupled moisture and heat transport in deforming porous media comprises 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.

Initial and boundary conditions
The initial conditions specify the full fields of gas pressure, capillary or water pressure, temperature and displacement and velocities: where n is the unit normal vector of the surface of the porous medium, r ¥ gw and T ¥ are the mass concentration of water vapor and temperature in the undisturbed gas phase far away from the interface, and q ga , q gw , q w 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: where t is the imposed traction.

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 equations are discretized in space by means of the finite element method, yielding a non-symmetric and non-linear system of partial differential equations: The Eqs. (34) can be rewritten in concise form as where , , , , C(X) is "the general capacity matrix", 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.

Method of solution
Coupled mechanical and transport processes after discretization by the finite element method are described by a system of ordinary differential equations which can be written in the form 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 system of ordinary differential Eqs.(36) is based on expressions for unknown values collected in vector r at time n + 1 r r v where the vector v n +a has the form The vector v contains time derivatives of unknown variables (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 This formulation is suitable because an explicit or implicit computational scheme can be set by parameter a.The advantage of the explicit algorithm is based on possible efficient solution of the system of equations, because parameter a can be equal to zero and capacity matrix C can be diagonal.Therefore the solution of the system is extremely easy.The disadvantage 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 complicated for nonlinear problems, where a nonlinear system of algebraic 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 where vectors f int and f ext contain internal values and prescribed values.Due to the nonlinear feature of the material laws used in the analysis, the Eq. ( 40) is not valid after computation of new values from the equations summarized in Table 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 established.This is done with the help of matrices B s and B g where particular partial derivatives are collected.Concise expressions of strains and gradients are written as New stresses and fluxes are obtained from the constitutive relations where D s is the stiffness matrix of the material and D q is the matrix of conductivity coefficients.The real nodal forces and the discrete fluxes on an element are computed from the relations   method is used.If the matrices are updated only once after every time step, the modified Newton-Raphson method is used.

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 compared after 1, 2 and 4 hours.

Conclusion
Hydro-thermo-mechanical problems are extremely complicated due to many nonlinearities, and therefore only numerical 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-dimensional 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 Acta Polytechnica Vol.44 No. 5 -6/2004 Transport processes lead to nonsymmetric indefinite systems.This means that usual methods, such as LDL T decomposition, do not work for such problems, and a more general LU decomposition must be used.It seems to us that there are where pivoting is necessary.

Fick
where D g (m 2 ×s -1 ) is the effective dispersion tensor.It can be shown [

} m 1
of a porous skeleton associated with the grain rearrangement can be expressed using the constitutive equation written in the rate form & ( is the tangential matrix of the porous skeleton and & equilibrium equation (the linear balance equation) must still be introduced to complete a set of governing equations div flux boundary conditions for water 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 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 ( ) matrices C and K are updated after every computation of the new increment from Eq. (44), the full Newton-Raphson Acta Polytechnica Vol.44 No. 5