Open access peer-reviewed chapter - ONLINE FIRST

A Combination of Finite Difference and Finite Element Methods for Temperature and Stress Predictions of Early-Age Concrete Members

By Tu Anh Do

Submitted: June 9th 2020Reviewed: December 2nd 2020Published: December 26th 2020

DOI: 10.5772/intechopen.95331

Downloaded: 45


A combination of finite difference and finite element methods was employed to develop a model for predicting the temperature development and thermally induced stresses in early-age concrete members (such as bridge footings, piers, columns, girders, and slabs). A two-dimensional finite difference (FD) scheme was utilized for heat generation and transfer within a hydrating concrete member. A finite element (FE) plane strain model was then established to compute the thermal stresses in the concrete subjected to the temperature changes. The FD-FE model can be easily created using any programing language, and the methodology can be used to predict the temperatures and stresses as well as assess the possibility of early-age cracking in concrete members.


  • finite difference
  • finite element
  • early-age concrete
  • heat of hydration
  • thermal stress
  • thermal cracking
  • insulation layer

1. Introduction

Thermal cracking is one of the biggest concerns regarding early-age concrete. Hydration of a large amount of cement results in higher peak temperatures as well as larger temperature differences between the concrete surface and the core. Such large temperature differentials can cause substantial tensile stresses that might increase the likelihood of early-age cracking in the concrete [1].

In order to control the temperature in early-age concrete structures, thus mitigating the risk of thermal cracking, temperature and stress analyses should be performed beforehand. Different methods have been used for predicting the temperatures and thermal stresses concrete structures at an early age. Among them, the Schmidt’s method is a simple approach but has been widely used for computing the temperatures for single nodes in the concrete [1]. The finite difference (FD) method was also employed in spreadsheet programs [2, 3] or in computer programs [4, 5] for calculating temperature–time histories in concrete elements. A two-dimensional model for thermal analysis based on the finite volume method (FVM) was introduced by Yikici and Chen [6]. The finite element (FE) method has been commonly utilized for both thermal and stress analyses of early-age concrete structures [7, 8, 9, 10, 11, 12, 13].

This chapter presents a two-dimensional FD scheme for thermal analysis of a concrete element. An FE analysis was then used to calculate the temperature-induced stresses in the concrete. The analysis results were compared with measurements of actual concrete elements. The combined approach can be a simple and useful tool for analyzing temperatures and thermal stresses in early-age concrete elements.

2. FD scheme for solving heat transfer

The heat evolution and temperature in a concrete element can be known by solving the governing differential equation as described in Eq. (1):


where ρ is density; cp is specific heat; T is temperature; t is time; k is thermal conductivity; x, y, and z are coordinates; and is heat evolution rate.

The finite difference formulation for any node in the system can be written as [14]:


where i = rate of heat conduction at time step i; Ėi = rate of heat generation at time step i; Tim and Ti+1m = temperatures of node m at time step i and i + 1, respectively; and Δt = time interval.

During the actual construction stage of concrete structures, the concrete is usually covered by formwork and/or insulation materials. Heat generated from cement hydration is conducted through the formwork and/or insulation layer before being dissipated to the surroundings by surface convection (Figure 1).

Figure 1.

FD mesh for heat conduction of concrete covered with insulation layer.

Considering a formwork/insulation layer covering the concrete, and assuming a unit square mesh for the concrete (Δx = Δy = l) and an insulation thickness of d (Figure 2), the FD formulation for the interior node can be computed using Eq. (3).

Figure 2.

FE plane triangular element.


where Tim,n = temperatures of node (m,n) at time step i; Tm1,ni,Tm+1,ni,Tm,n+1i,Tm,n1i= temperatures at neighboring nodes; and τF = dimensionless Fourier number,


Using Eq. (2), the FD equations for each of the four outer corner nodes of the insulation can be derived. For instance, the quarter size volume element of the insulation layer (d × d × 1) represented by the top left outer corner node (1,N) is subjected to convection on both sides and to conduction from the right and bottom nodes, the energy balance relation above (Eq. (2)) becomes:


where ρs = insulation material density; cps = specific heat of insulation material; ks = thermal conductivity of insulation material; Ta is the ambient temperature; h is the convection coefficient; and


The insulation volume element at the surface node (2,N) adjacent to the top left outer corner node is subjected to convection at the top and conduction at the left, right, and bottom surfaces. An energy balance on this element gives:

T2,Ni+1= T2,Ni12hd+lA2ks2A22dlA22d+lA2d+2A2T1,Ni+dlT3,Ni+d+ldT2,N1i+hd+lksTaE7



Similarly, the insulation volume element of half size at a surface node is subjected to convection at the top and conduction at the left, right, and bottom surfaces. An energy balance on this element gives:




The “mixed” volume element at the concrete’s corner node (2,N-1) is subjected to conduction at the four sides. An energy balance on this element gives:




The “mixed” volume element at a concrete’s top surface node is also subjected to conduction at the four sides. An energy balance on this element gives:




It is noted that the stability criterion of the explicit method requires all primary coefficients to be positive or zero for all nodes:


The maximum time step used to solve the problem must satisfy Eq. (15) above.

If an insulation layer is not used, the new corner temperature Tim,n+1 will be simplified to:


and the next time step temperature of a top surface node will be simplified to:


The maximum time step in this case is as follows:


2.1 Rate of hydration heat

The rate of heat liberated from cement hydration depends on the temperature of the concrete element itself. The heat rate can be experimentally determined using isothermal [10, 15], adiabatic [16, 17], or semi-adiabatic calorimetry [4]. The experimental adiabatic temperature rise (ATR) can be converted into a maturity-based heat rate as presented by Ballim and Graham [18], in which the total heat (Q) liberated at any time (t) is firstly computed from the ATR using the following relationship:


where Tt = sample temperature at time t; T0 = initial sample temperature; ms = mass of concrete sample; and mc = mass of the cementitious materials in the mix. The heat rate in the adiabatic condition is then calculated by differentiating Eq. (19):


The “maturity heat rate” (qte), as shown in Eq. (21), is used in a further thermal analysis of concrete, which considers the maturity of concrete.


where te is equivalent age (or maturity) [19]:


where Ea is apparent activation energy (J/mol); R is the universal gas constant (8.314 J/mol-K); Tc(t) is concrete temperature (K); and Tr is reference temperature (K).

The activation energy (Ea) of a cement blend can be estimated from its chemical compositions using the following relationship derived by Poole [20]:

Ea=41230+1416000pC3A+pC4AFpcempSO3pcem347000pNa2Oeq        19.8Blaine+29600pFApFACaO+16200pslag51600pSFE23

where pFA = % fly ash in the cementing blend; pFA-CaO = % CaO in fly ash; pslag = % slag in the cementing blend; pSF = percentage of silica fume in the cementitious materials; Blaine = cement fineness (m2/kg); pi = percentage of i component in the cement (C3A, C4AF, SO3, cem = cement); and pNa2Oeq = % Na2Oeq in cement (= 0.658 × %K2O + %Na2O).

The actual heat rate, which will be used in a numerical model, can be reconstructed from the maturity heat rate using the following equation:


The maturity-based heat rate curve qte should be built from an isothermal or adiabatic test, before the actual heat rate can be computed at each time step for the analysis [18]. The drawback of this method is that the total time of the constructed maturity-based heat rate is limited by the test duration.

There are several models to mathematically characterize the heat generation from the cement hydration. The 3-parameter exponential degree of hydration model show in Eq. (25) [21] has been widely used for predicting temperature development in concrete since it includes the temperature effect through the equivalent age:


where αu is ultimate degree of hydration; τ and β are hydration parameters.

The total cumulative heat Q(te) is proportional to the degree of hydration α(te) as expressed in Eq. (26). The rate of heat generation with respect to equivalent age and real age can be determined using Eqs. (27) and (28), respectively.


where Qc is the total available heat (J/m3).

The hydration parameters (αu,τ and β) can be determined from the fitted curve, Eq. (25), using the experimental ATR data. These parameters can also be calculated from an experimental isothermal cumulative heat curve without converting the real time into the equivalent age because in the isothermal condition (i.e., at a reference temperature of 23°C), the test time is identical to the equivalent age.

3. FE method for solving thermal stresses

Since a common concrete structure has one dimension larger than the other two, the middle cross section should be analyzed; hence, a FE plane strain problem is selected for the stress computation. A triangular element is chosen with nodes i, j, m numbered in a counterclockwise order as illustrated in Figure 2 [22]. The strain at any point within the element is estimated by Eq. (29):


where ae = element displacement vector, and


in which ai=xj ymxmyj; bi=yjym; ci=xmxjwith the other coefficients obtained by a cycle permutation of the subscripts in the order i, j, m; and Δ is area of the triangle.

The stress vector in the element can be calculated as:




and the thermal strain is derived as [22]:


in which ν = Poisson’s ratio, αc = coefficient of thermal expansion, and θ e = temperature change (from the previous time step to the current time step) subjected to the element. The element stiffness matrix ijm is calculated using the following equation:


where tt = element thickness.

The nodal forces due to thermal strain is computed as follows:


in which E = elastic modulus. The nodal displacement vector U is derived by solving the global system of equations:


Computational Procedure

  • FD thermal analysis:

    1. Define geometry of the structure (including the nodal grid), initial material properties, initial temperature and boundary conditions, and time interval.

    2. Compute the nodal degree of hydration and the rate of heat evolution.

    3. Compute the new temperature at each node.

    4. Iterate (2) & (3) and record the temperatures.

  • FE stress analysis:

    1. Divide the nodal grid into triangular elements (the vertices coincide with the FD grid nodes).

    2. At t = n (n = 1, 2,…), calculate average temperature, equivalent age and degree of hydration of each element.

    3. Let i = 1, compute each element’s effective modulus.

    4. Compute element stiffness matrix, global stiffness matrix, and equivalent nodal forces; solve for nodal displacements and element stresses.

    5. Let i = i + 1 and iterate (7) and (8) till i = n. Sum all the stresses at step (8) to get the total stress.

    6. Let n = n + 1. Iterate (6) through (10) until the final time step is achieved.

4. Temperature analysis of bridge pier footing

A bridge pier footing constructed in Orlando, Florida was monitored for temperature development within 7 days after casting (Figure 3). The concrete footing had dimensions of 6.71-m × 3.05-m × 1.75-m and was insulated with 25.4-mm thick polystyrene foam boards at its bottom, top, and sides during 7 days.

Figure 3.

Bridge footing for monitoring in Orlando, Florida.

The cementitious materials of Mix #1 was experimentally measured for the heat of hydration using an isothermal calorimeter. The hydration parameters and calculated activation energy (Ea) for Mix #1 are presented in Table 1.

Mixτ (h)βαuQc (J/m3)Ea (J/mol)
#116.730.87648.3141.26 × 10835,451
#214.00.940.7031.67 × 10841,800

Table 1.

Hydration parameters and activation energy.

The concrete had a density of 2238 kg/m3, specific heat of 1045 J/kg-K, and thermal conductivity of 1.87 W/m-K. The footing was insulated with Styrofoam that has density, thermal conductivity, and specific heat of 16 kg/m3, 0.04 W/m-K, and 1200 J/kg-K, respectively [14]. The boundary conditions consist of the initial temperature and the external ambient temperature over time. The air convection coefficient at the insulation surfaces was assumed to be 13.9 W/m2-K [5, 23].

Two temperature sensors were installed at the center and at the center top surface of the footing to record temperatures within 7 days after placement. The measured temperatures at the center and top, and the ambient temperature are shown in Figure 4. A peak measured temperature of 74°C occurred in the middle 42 hours after concrete casting. Figure 5 shows the temperature distribution in the footing at 40 h calculated by the FD model. The predicted FD temperatures at the center and the top of the footing are also plotted in Figure 4. It is clear that the temperature histories computed using the FD model show very close agreement with those collected in the field.

Figure 4.

Predicted and measured temperature profiles in the footing.

Figure 5.

FD temperature contour in the footing at 40 h (°C).

5. Temperature and thermal stress predictions of cap beam

A bridge concrete cap beam (pier cap) was analyzed for temperatures and thermal stresses due to the heat of cement hydration. The cross section of the pier cap was 1.6-m by 2.1-m. The concrete used in the cap beam is Mix #2 with the hydration parameters and activation energy listed in Table 1. The concrete coefficient of thermal expansion (CTE) of 8.5 × 10−6/°C, density of 2287 kg/m3, the specific heat of 1028 J/kg-K, and thermal conductivity of 1.87 W/m-K were assumed in the analysis.

Figure 6 shows the calculated temperature profiles at the core, corner and side of the section. The center temperature peaked at 70.7°C approximately 28 hours after casting. The temperature contours are also depicted in Figure 7.

Figure 6.

Temperature profiles at different points in the cap beam.

Figure 7.

Temperature distribution of the section at 30 h.

The thermal analysis was followed by a stress calculation using the FE model with the element mesh shown in Figure 8a. The 1st principal stress and stress component σyy contours are shown in Figure 8b and c, respectively. The figure reveals that the maximum stress is σyy occurring at the mid-sides and having almost the same magnitude as the 1st principal stress.

Figure 8.

FE mesh and stress distributions (MPa) in pier cap at 21 h.

The calculated stresses over time at different locations of the pier cap are plotted in Figure 9. Clearly, the maximum stress is σyy at the mid-sides, while σxx is the maximum stress at the corner (compared to σyy), thus the middle sides have a higher risk of cracking.

Figure 9.

Calculated stresses in the pier cap.

To assess the model’s accuracy, the computed stress-time histories are compared with those obtained from the 3-D ABAQUS FE model developed by Lin and Chen [12]. It is worth noting that the ABAQUS model was validated using measurements on 2 concrete blocks. Figure 9b shows that the 2-D FE results reasonably match with those of the 3-D ABAQUS model.

The 3-D ABAQUS results reveal that the maximum stress is the component σzz at the corner. Nevertheless, the 2-D FE model cannot compute this stress component because it is out-of-plane. The plane element over-predicts σ xx at the corner compared with that of the 3-D ABAQUS. The maximum stress σzz at the corner predicted using the 3-D ABAQUS model is about the same magnitude as σyy at the mid-sides, hence the critical stress magnitude as well as cracking risk can still be forecast by the 2-D analysis.

6. Conclusions

In this study, FD and FE formulations were created for solving the transient heat transfer equation and thermal stresses in a concrete element. The results of this study show that the approach that combines the FD and FE methods can be a useful and effective tool for predicting temperature evolution and thermally induced stresses in early-age concrete members with simple geometries. The FD model can analyze thermal behavior of a concrete placement covered with formwork or an insulation layer, thus it can help engineers/contractors control concrete temperature during construction.


This work is financially supported by the Ministry of Transport of Vietnam. Special thanks are given to Prof. H. L. Chen and Mr. G. Leon at West Virginia University (WVU) for their contributions to this study.

Download for free

chapter PDF

© 2020 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Tu Anh Do (December 26th 2020). A Combination of Finite Difference and Finite Element Methods for Temperature and Stress Predictions of Early-Age Concrete Members [Online First], IntechOpen, DOI: 10.5772/intechopen.95331. Available from:

chapter statistics

45total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us