Computational Model for Internal Relative Humidity Distributions in Concrete

A computational model is developed for predicting nonuniform internal relative humidity distribution in concrete. Internal relative humidity distribution is known to have a direct effect on the nonuniform drying shrinkage strains. These nonuniform drying shrinkage strains result in the buildup of internal stresses, which may lead to cracking of concrete. This may be particularly true at early ages of concrete since the concrete is relatively weak while the difference in internal relative humidity is probably high. The results obtained from this model can be used by structural and construction engineers to predict critical drying shrinkage stresses induced due to differential internal humidity distribution. The model uses finite elment-finite difference numerical methods. The finite element is used to space discretization while the finite difference is used to obtain transient solutions of the model. The numerical formulations are then programmed in Matlab.The numerical results were compared with experimental results found in the literature and demonstrated very good agreement.


Introduction
Cracking is detrimental to the serviceability, durability, and the aesthetic quality of concrete structures.A major driving force behind cracking is the nonuniform early-age drying shrinkage [1].Drying shrinkage can be defined as the volume reduction that concrete suffers as a consequence of the moisture migration when exposed to a lower relative humidity environment than the initial one in its own pore system.During production of concrete, more water is added to the concrete mix than necessary for hydration for the sake of workability.This leads to having two types of moisture in concrete.The first type is structural water, chemically bound within the cement paste.As the concrete hydrates, some shrinkage takes place in absence of additional water as hydration takes up some free water.This shrinkage is referred to as autogenous shrinkage and is typically about 50 to 100 microstrains [2].The second type is the excess water that does not take part in the hydration product and as a consequence it will not be chemically bound to the solid phase.This water is contained in the pore structures.
Drying shrinkage in concrete is due to drying of the water contained in the pore structures and the associated decrease in moisture content.Therefore, it is necessary to estimate the moisture loss as accurately as possible in order to study drying shrinkage in concrete members.Accurate prediction of drying shrinkage in concrete requires knowledge of the concrete moisture content and its variation.Extensive research has been conducted in literature to identify the drying mechanisms of water contained in the pore structures and different mechanisms have been proposed in the literature [3][4][5].
The most prominent drying mechanisms proposed in the literature include surface free energy, capillary tension, movement of interlayer water, and disjoining pressure [6].These mechanisms are complex and are often interrelated.This is mainly due to the wide range of pore sizes induced in concrete, which determines the different transport mechanisms during drying.The pore sizes also keep on changing with time because of the continuing hydration process.Drying shrinkage is better expressed in in terms of internal relative humidity than moisture concentrations since the former is continuous while the latter is discontinuous across the surface and interior of the concrete.However, it is worth noting that the water loss is approximately proportional to the drop in the average relative humidity in the pores, and the drying shrinkage is approximately proportional to the water loss; the average humidity also evolves in time similarly as shrinkage [4].
When concrete is exposed to a lower environmental relative humidity than its internal one, the vapor pressure of the water remaining in the region losing water decreases progressively with the moisture content.It is this tendency of the vapor pressure to equalize that causes the slow drift of moisture toward the drying surface.Therefore, the drying process starts at the surface that is exposed to drying and gradually penetrates into the concrete [7].During the drying process, whether it is due to internal or external restrictions, self-equilibrated stresses are usually generated within a concrete member cross-section.
The relative humidity gradients are responsible for a differential drying of the concrete member, causing tensile stresses near the exposed surface and compressive stresses in the inner layers [8].When the induced-tensile stresses exceed the tensile strength of concrete, which is an age-dependent property and usually low at early ages when differential drying is high, cracking will occur [9].This phenomenon affects the durability of many civil engineering structures.Hence, accurate prediction of the drying process will enable structural designers to estimate the drying shrinkage strains a concrete member experiences and to take appropriate steps in design to accommodate these strains.
Moisture loss from concrete is influenced by external factors and internal factors.The external factors affecting loss of moisture from concrete are ambient conditions and size and shape of the concrete member under consideration.Ambient conditions that affect the loss of moisture from the concrete surface consist of relative humidity, wind velocity, and air temperature.Higher drying shrinkage is expected with decrease in relative humidity, with rise in ambient temperature, with increase in air movement around the concrete, and with increase in the length of time for which concrete is subjected to drying conditions.The internal factors affecting drying shrinkage of concrete are those related to its constituents: cements, aggregates, and admixtures; concrete mix design; water-cement ratio and water content; and aggregate properties and volume fraction and those related to the construction of the concrete such as placing, compaction, and curing.
When concrete structures are directly exposed to the environment, the top surface of the structures are usually more susceptible to drying shrinkage than the bottom surface.Hence, the bottom part of these structures show higher moisture content than the top part.This typically creates a drying shrinkage differential through the slab thickness.The existence of a moisture gradient can greatly increase the tensile stresses near the drying surface, potentially exacerbating crack formation [10].
The determination of drying shrinkage stress gradients involves knowledge of internal RH values as a tool to quantify the state of internal stress in concrete.Kelvin-Laplace equation can be used to transform internal relative humidity distribution to stress based on the relationship between the negative pore fluid pressure and the internal relative humidity in the pore structure [11].Hence, determination of the internal relative humidity is of a paramount importance so as to determine the drying shrinkage and/or ultimately the cracking susceptibility of concrete.This is especially important for mass concrete structures like bridge piers, dams, thick concrete pavement, and so forth.

Mathematical Model for Moisture Diffusion
Developing a mathematical model that captures all the aforementioned mechanisms of drying process, which may occur simultaneously at the same region or different regions of the concrete mix, is very complex and beyond the scope of this paper.However, a simplified but reasonably accurate model is developed based on factors affecting rate of moisture transmission.Moisture moves through concrete in a partially adsorbed or condensed state by diffusion.The rate of moisture transmission depends on the degree of saturation, which is a function of the relative humidity on each side of the concrete.Therefore, the major driving force for water vapor movement through a concrete member is the relative humidity differential through the member's depth [12].
In this paper, a formulation based on Fick's law will be used to model the internal humidity distribution.In this formulation, the very important material parameters are the rate of moisture diffusion through the concrete and the diffusion coefficient ().The diffusion coefficient has been developed from experimentally determined relative humidity (RH) profiles in various concrete slabs during drying with different water to cement ratio, the ambient temperature and humidity [13]. is dependent on numerous external and internal factors, such as ambient temperature and humidity; however, the main factors affecting  are the internal relative humidity itself and water to cement ratio.
In this paper, surface factor () represents the cumulative effect of ambient conditions, namely, ambient temperature and wind speed, the / ratio of the mix, and the method of finishing, or closing, of the exposed surface.Surface factor will be used to model the moisture loss through the concrete boundaries.It is simply an equivalent convection coefficient for moisture losses at the concrete boundaries.
In concrete, moisture migrates from regions of high concentration to regions of low concentration.This suggests that the moisture flux, , is proportional to the gradient of the pore relative humidity and is expressed by [7] where ℎ is the pore relative humidity and  is the permeability.The specific water content, , is the function of pore relative humidity, ℎ, in the desorption isotherm; that is,  = (ℎ), so that the mass balance equation can be expressed by [14] where  is the moisture diffusion coefficient and is defined by This indicates that the moisture diffusion coefficient itself is dependent on the relative humidity.
The moisture diffusion coefficient is expressed by (4) as a function of the pore relative humidity, which has been developed from experimentally determined RH profiles in various concrete slabs during drying for 0 < ℎ < 1 [14]: where  1 represents (ℎ) when ℎ = 1.0,  is parameter that represents the ratio  0 / 1 , where  0 is the minimum of (ℎ) for ℎ = 0.0, ℎ  is the pore relative humidity at (ℎ) = 0.5 1 , and  is an experimental parameter. 1 can be estimated from  1 =  1,0 /(  / 0 ), where  1,0 = 3.6 × 10 −6 m 2 /h,  0 = 10 MPa, and the characteristic compressive strength   can be estimated by the mean compressive strength.Alternatively some authors have suggested that  lies between 10 −11 and 10 −9 m 2 /s for normal dense concretes with / ratios between 0.4 and 0.6 [13], which of course reduces the mathematical complexity in solving the differential equations.
The nonlinear differential equation governing the moisture diffusion in two dimensions is expressed by [13] where ℎ(, , ) is the moisture content at (, , ), (ℎ) is the moisture diffusivity as explained below, and  is the time from the start of diffusion process.
To solve the governing differential equation, one needs initial and boundary conditions.Initial boundary condition is given by ℎ (, , ) = ℎ 0 (, , ) .
Prescribed moisture content at boundary is given by The boundary condition on sealed or impermeable boundary is given by The boundary condition on surface of the concrete is given by where ℎ/ is moisture gradient at the drying surface identified by a unit normal, , ℎ en is the equilibrium humidity that an element would reach given particular environmental condition (environmental humidity), and ℎ  is the humidity at the drying surface. is the value of the surface factor or convection coefficient, which depends mainly on ambient temperature and wind speed, the / ratio of the mix, and the method of finishing, or closing, of the exposed surface.
As can be seen in ( 9), the moisture loss through the surface is expressed by the convection coefficient (surface factor), which is driven by the difference between the boundary (surface) and the environmental humidities.Typical values suggested range between 10 −8 and 10 −7 m/s [13].
Figure 1 shows surface factors used by West and Holmes in the FE analysis for the slabs in a laboratory and in a controlled room.This was found by a parametric study to best suit the experimental conditions and results for calibrating models.This helps to have the insight of what value to use while validating the test results as there are only few literatures which quantified the surface factors.The initial rate of evaporation from both slabs is quite high but, over time, it decreases as the RH differential between the concrete and the ambient air slowly converges while, eventually, drying will cease, practically.

Numerical Solution
The numerical solution scheme used in this study is based on the Taylor-Galerkin finite element method in combination with finite difference method.A linear element is used for the one dimensional elements at the boundaries of the concrete and a bilinear element is used for the two dimensional elements inside the boundaries.Upon applying the numerical methods at the governing moisture diffusion equation, (10) is derived as follows: where [] is moisture content stiffness, [] is the capacitance matrix, {} is the total moisture load vector due to moisture fluxes, {ℎ} is the nodal relative internal humidity vector, Journal of Computational Engineering and {ℎ/} is the rate of change of nodal internal relative humidity vector.
For time dependence problems, where temporal discretization is needed, it is effective to employ numerical solutions in the time domain to get the solution of the differential equation [15].From the mean value theorem for the differentiation, ( 11) is derived as follows: where Δ is the time step and {ℎ} +Δ and {ℎ}  are the vectors containing nodal internal humidity values at times  + Δ and , respectively.
The vector {ℎ} at time  =  * within the time step Δ is given by where  = ( * −)/Δ and the vector {} is equal at  and +Δ.
Rewriting (12) in terms of {ℎ}  and {ℎ} +Δ results in Equation ( 13) gives nodal values {ℎ} +Δ in terms of a set of known values, {ℎ}  and the ratio, .There are four different types of finite difference schemes depending on the value of .In this study, the backward finite difference method, for which  equals 1, is used because a Fourier stability analysis showed that the numerical solutions of ( 13) are unconditionally stable for  ≥ 1/2 [16].Hence, (13) further reduces to A Matlab program is developed to obtain numerical solution for (14).

Relation between Relative Humidity and Shrinkage Stress
Determination of drying shrinkage stress is not the purpose of this paper; however, the relationship between internal relative humidity and drying shrinkage is depicted in this section.Once the internal humidity distribution is known, the drying shrinkage strains at any location within the boundary on a concrete can be calculated using [7]  sh =  sh   (ℎ) =  sh (1 − ℎ) , where  sh is the shrinkage coefficient, which can be determined from drying shrinkage on complete drying and the ratio of elastic modulus with time,   (), and   is function of pore relative humidity, ℎ.
The stresses induced due to drying shrinkage under fully constrained boundary condition at a time can be computed using Hook's law as expressed in  sh () =   ()  sh =   ()  sh   (ℎ) . (16)

Results
The main output of the internal RH prediction program is the internal RH values of hardening concrete at each node of the finite element mesh at each time step.In this paper a regular rectangular in shape concrete structure is considered.Comparison from previously done laboratory tests, which were found from the literature, is also conducted with the corresponding program outputs.
Figure 2 shows part of the mesh and the selected 21 nodes along the depth of the concrete section to evaluate temperature development with depth versus time.Table 1 shows the general geometrical, finite element mesh, environmental input parameters used for evaluating the numerical model results based on theoretical scenarios.The Cauchy's boundary coefficients are obtained by dividing the surface factor (the moisture loss coefficient) on the diffusivity of concrete.The mesh associated with the geometric conditions has 400 elements and 441 nodes.Figure 3 shows the plot of the internal RH values with time for selected nodes on the left edge of the concrete when the boundary conditions are as given in Table 1. Figure 4 shows the plot of internal RH values with time under the same conditions and same nodes as Figure 3 when all edges except the top are insulated.As intuitively expected, the internal RH values in Figure 3 are generally lower than the values in Figure 4.This is due to the change in boundary conditions; insulating the left and right edges reduces the rate of moisture loss significantly.
Figures 5 and 6 show the internal RH distributions at different location at the same instant of time, also known as iso-time plots.These plots may be utilized to identify whether the short term or long term internal RH distribution that is inducing critical differential drying shrinkage.
Figure 5 displays the iso-time plots for the 21 nodes on the top surface of the concrete under parameters mentioned in Table 1. Figure 6 displays the iso-time plots for the 21 nodes on the left edge of the concrete under parameters mentioned in Table 1.While about a 7% maximum internal RH difference is observed on the second day in both the left edge and the top surface, about a 5% maximum internal RH difference is observed on the fifth day in both the left edge and the top surface.

Verification
In this paper, the internal RH values measured by Yuan and Wan [17] from two laboratory specimen are used for verification of the developed internal RH prediction model.
The cement used in this test is ordinary Portland cement (ASTM Type I), and the fineaggregate is river sand.The crushed granite gravel passing the 19 mm sieve was used as thecoarse aggregate.Table 2 shows the detailed mix proportions of the specimen.The moisture diffusion constants,  1 , for the two mixes are determined as 1.81 × 10 −6 and 1.0 × 10 −6 m 2 /h, respectively [14]. Figure 7 shows the geometry and size of specimen and the depths of the RH gauges placed.
The actual surface factor for the environmental condition the specimen was kept in is unknown.Hence, in this paper, the first specimen (mix I) is used to determine the surface factor by applying iterative trial and error.The surface factor is obtained as 3.5 × 10 −4 m 2 /hr.The specimen is exposed only in the top surface; all the other surfaces are insulated.
Figure 8 shows internal RH values predicted from the model and measured from laboratory specimen (mix I) under the surface factor mentioned.
Figure 9 shows internal RH values predicted from the model and measured from laboratory specimen (mix II) at the four gauge locations.

Conclusions
A finite element-finite difference model is developed for predicting nonuniform internal relative humidity distribution in concrete.The results from the model show that the internal RH significantly varies according to the depth from drying surface; this directly results in differential drying shrinkage and may cause surface cracks.Environmental humidity, the diffusivity of concrete (mix proportions), the surface factor, size, and the boundary condition of concrete are the main factors that affect the internal RH distribution in concrete.The model developed in this study will enable structural    designers to predict the internal RH and estimate the drying shrinkage strains a concrete member experiences and to take appropriate steps in design to accommodate these strains.
The following can be stated in conclusion of the moisture diffusion mode.It has been observed that the internal RH significantly varies according to the depth from drying surface; this directly results in differential drying shrinkage which in turn may cause surface cracks.Thus, knowledge of predicting the internal RH would help contractors and design engineers to avoid drying shrinkage cracks in thick concrete structures.The ambient conditions, the diffusivity of concrete (mix proportions and material types), the surface factor, size of the concrete, and the boundary condition are some of the factors affecting the internal RH distribution in concrete.The model developed in this study is important in predicting the internal RH and is helpful for design engineers and contractors to avoid drying shrinkage cracks, especially in thick concrete structures.

Figure 2 :
Figure 2: Part of the mesh and the 21 nodes used for evaluating the internal RH values.

Figure 3 :
Figure 3: Internal RH at selected nodes on the left edge when only the bottom surface is insulated.

Figure 4 :
Figure 4: Internal RH at selected nodes on the left edge when the bottom, left, and right surfaces are insulated.

Figure 5 :
Figure 5: Iso-time internal RH distribution on top surface of concrete (node 1 is the far left node).

Figure 6 :Figure 7 :
Figure 6: Iso-time internal RH distribution on left edge (node 1 the bottom most node).

2 Figure 8 :
Figure 8: Internal RH values predicted from the model and measured from laboratory specimen (mix I).

Figure 9 :
Figure 9: Internal RH values predicted from the model and measured from laboratory specimen (mix II).

Table 1 :
General input parameters used.