1 Introduction and scope of simulation

This document aims at providing recommendations for thermo-chemo-mechanical modelling of massive concrete structures, with particular consideration to the assessment of the risk of cracking. It intends to serve designers, practitioners, analysts and researchers with a set of consensus of good practices to apply on several important aspects that can significantly affect the quality and realism of stress and crack risk predictions in practical applications. The scope of application needs to be clearly defined, as all the recommendations will be directly aligned. This document applies to massive concrete structures, understood as concrete structures with significant volume or thickness, which can cause significant temperature variations at early ages, owing to the heat of hydration release of cement (for definitions and identification of a “massive” concrete structure the reader is referred to Ref. [1]). The duration scope of analysis is limited to the period when internal temperatures of concrete are still influenced by the temperature rise inherent to cement hydration. Examples could be thick walls, thick foundations (such as pile caps and rafts), thick slabs, large columns, dams, etc. The numerical analysis under investigation pertains to concurrently conducting time-dependent thermal analysis (heat-transfer), chemical analysis (cement hydration) and mechanical analysis (stress analysis), which for simplicity will be referred to as thermo-mechanical analysis thereafter. Also, although the term thermo-chemo-mechanical analysis can be considered as more descriptive, the inherent chemical modelling approach is phenomenological and contains no detailed information of microstructure evolution. It should be noted that the recommendations given herein are based on the experts group’s theoretical and application experiences for a context of simulation which is affected by a profound number of factors.

On the perspective of simulation, the following specific matters are presumed: (i) a staggered thermo-mechanical simulation is made, with one-directional coupling being considered; (ii) the simulation is made with the widely-used finite element method.

This document is aimed at recommending practices to those that are knowledgeable about the basic principles of thermo-mechanical simulation at early ages, and hence does not pinpoint the basics of all aspects covered. This strategy has been adopted both for the sake of brevity and also for the sake of facilitating the evidencing of actual recommendations in the document. For basic knowledge on phenomenology and simulation, the reader should refer to up-to-date literature, such as the State-of-the-art Report of RILEM TC 254-CMS on Thermal Cracking of Massive Concrete Structures [2].

For the sake of clarification of governing equations and terminology, the classical thermal energy balance equation implied in this document is forwarded (further relevant fundamental theory can be accessed in [3]):

$$\frac{\partial }{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)+\dot{Q}=\rho c\dot{T}$$
(1)

where \(k\), \(\rho\) and \(c\) is the thermal conductivity, specific weight and specific heat capacity of concrete, respectively, \(T\) and \(\dot{T}\) is the temperature and its derivative with respect to time, respectively, \(\dot{Q}\) is the heat generation rate and \(x\), \(y\), \(z\) are the spatial coordinates.

The mechanical simulation relies on simplified generally-used assumptions, which are already implemented by a significant number of commercial software in the market, and reflect the typical intents of those using these recommendations:

  1. 1.

    Modelling disregards cracking, as the intent is to assess cracking risk, through analysis of stress levels or cracking indexes (e.g., ratio of tensile stress to tensile strength). Indeed, if the cracking risk is high, then some cracking will probably be occurring; however, the simulations under investigation herein focus on pre-cracking stages and disregard local cracking in stress singular points around notches, openings or local forces. It is, therefore, not necessary to account for concrete cracking in this type of analysis;

  2. 2.

    The evolving mechanical properties of concrete are following a well-defined evolution (e.g., age/hydration dependent tensile strength gain), normally taking into account maturity (e.g. Equivalent age concept);

  3. 3.

    Viscoelastic behaviour, e.g., creep, is directly considered in the modelling;

  4. 4.

    Reinforcement effects may be omitted but included in exceptional cases (e.g. very heavily reinforced elements).

The concept of degree of hydration has been forwarded to describe the advancement of the chemical reactions involved in cement hydration. It is normally assumed that the degree of hydration may be defined as the ratio between the amount of cement that has reacted at a given instant, divided by the original amount of cement. In the scope of thermo-mechanical simulations for concrete structures, the degree of hydration is usually expressed in a simplified manner as the degree of heat development, which is the ratio between the total heat that has been released up to a given instant in time, divided by a reference cumulative heat (in some cases taken as the ultimately attainable accumulated heat release by the mixture in real conditions, and in other cases by the ultimately attainable heat release if all cement particles would be able to react totally). It has been acknowledged that the ‘degree of heat development’ concept has important merits in the scope of the early ages of concrete (normally not more than a few days/weeks), but tends to be less accurate when the scope of simulation extends to longer term analyses (months or even years). For longer term, the degree of hydration concept should be used only if an appropriate phenomenological law (not based on heat development) is used to predict the continuation of hydration.

The equivalent age approach (the one based on the Arrhenius law, with binder-dependent “apparent” activation energy) has demonstrated similar merits to the degree of heat development approach in thermo-mechanical simulations, allowing proper estimation of properties of concrete within time spans outside the scope of the early ages. While other traditional methods exist for maturity approaches (of which the equivalent age concept is part of), such as that of Nurse and Saul, they are not recommended as they are relatively obsolete and are not thought to be representing the relevant behaviour of the concrete adequately.

Another important aspect worthy of discussion, and that has important ramifications in the simulation methodology and results, concerns the choice of the instant at which the numerical simulation starts, here named as t0,analysis. Figure 1 below, which schematically depicts the evolution of E-modulus of concrete since the contact of water and cement in the mixer (tcontactwc), brings support to the discussion of this matter. The figure schematically depicts a period of a few days, with the time of the end of the mixing (tend-mixing) seeming overlapped with tcontactwc because of the small period of time spent in the initial mixing. After tend-mixing, there is a series of operations comprising transport from the mixing place to the casting place, and casting/vibration operations, which culminate with the end of casting at tend-casting. A few hours later, the final setting time is reached (tsetting), which normally corresponds to a few hours after tend-casting. For the thermo-mechanical analysis envisaged in this document, it is recommended that typically the simulation should start at t0,analysis = tend-casting, as to be able to simulate thermal interactions happening before the setting time and be more accurate in such concern. If the choice is such, care should be taken in adapting the necessary properties evolutions (hydration heat, stiffness, strength, creep, etc.) to match the conventional zero time. Care should also be taken in regard to the stiffness to consider for the concrete between tend-casting and tsetting. As the supporting formwork is not modelled in order to contain the ~ 0 stiffness concrete, very high deformations in the model would occur if one would actually consider the true initial stiffness of the concrete. It is, however, recommended in general to adopt a reference base stiffness (e.g., a value of 100 MPa or other appropriate value – see the end of Sect. 4 regarding self-weight) between tend-casting and tsetting. This initial value of stiffness is likely to be large enough to keep deformation of concrete at a low-enough level (comparable to that which would result from the explicit consideration of the bounding formwork). It is, however, also small enough for the initial stresses in concrete to keep negligible enough by comparison with the reality of a material that is not yet able to bear structural stresses before tsetting.

Fig. 1
figure 1

Key time instants at very early ages in the evolution of E-modulus of concrete

2 Geometrical considerations and modelling complexities

As a principle, one should start with a simple model, aiming modelling and computational time savings. The first step is the identification of symmetries in the problem, considering the geometry, boundary conditions and loading, both in view of the thermal and mechanical simulation fields. Good sense simplifications on the geometry/boundaries/loads are relatively frequent in thermo-mechanical simulations of massive concrete structures (e.g. disregarding the lack of symmetry of solar radiation as a load and modelling only half of the length of the structure). Furthermore, it is frequently advantageous to use symmetry simplifications as they tend to make the mechanical models more statically indeterminate and hence less prone to numerical instabilities upon the solution of the inherent equations.

The choice for 2D plane simulations in opposition to 3D modelling can bring significant computational time savings and yield very realistic results in terms of thermal field simulations (particularly in linear structures such as thick walls). It is a common mistake to choose 2D plane simulation for structures with important out-of-plane stresses, such as analysing only cross-section of a wall-on-foundation ignoring the effect of restraint change over the length of the wall and the effect of slip on the stress field in the wall [4, 5].

The consideration of formwork in thermal simulation should normally be made through an equivalent convective boundary coefficient, thus allowing the simplification of the simulation model without hampering realism (as the thermal capacity of the materials of the formwork can normally be considered negligible in regard to that of massive concrete). Examples of equivalent convective boundaries are given in chapter 3 of the TC 254-CMS state-of-the-art report [6] and discussed further in Sect. 4 of this document.

For heavily insulated sections, e.g., with insulating media exceeding thicknesses found in applications that are considered as normal for massive concrete (e.g., 17 mm plywood, 3 mm steel formwork or equivalent in addition to an indicative 25 mm expanded polystyrene sheet), it would be advisable to explicitly model them, as in such cases the results are more prone to modelling inaccuracies if an equivalent convective boundary is adopted. Nevertheless, as mentioned earlier, for massive concrete structures, the explicit modelling of the formwork and insulation becomes less important and can be modelled by an equivalent convective boundary, owing to the significant differences in thickness between the concrete and the formwork/insulation. Formwork removal or placement of additional insulation should, however, be captured by changing the heat transfer coefficient due to the high impact on the stress field.

Elements adjacent to the massive concrete structure under study (e.g. underlying terrain, connected structures or similar/different material) must be modelled explicitly whenever they are likely to influence either the thermal or mechanical field, or both. Examples of recommendations in this concern can be seen in Fig. 2.

Fig. 2
figure 2

Example of a wall-type massive concrete structure

The terrain underlying the analysed structural element should be explicitly modelled as a soil block in thermal analysis of the element and should not be simplified with a thermal boundary condition with an equivalent convective boundary coefficient. If such soil block is neglected, any replacing ‘equivalent convective boundary’ will inherently fail to model the heat storage effect of the terrain and hence cause inaccuracies in the numerical predictions of temperature fields of the actual concrete element. This can be considered as a common mistake in view of simplification of the model and corresponding analysis. The dimensions of the soil block should be big enough to fully capture the thermal field generated in the underlying soil over the whole analysed time history. For the case of a wall-on-foundation massive concrete structure, as a rule of thumb, the dimensions of the soil block can be initially taken after Fig. 2 as Ladd = BF and D = 2BF, and adjusted accordingly with subsequent simulations until it is ensured that the temperature at the outer boundaries of the block remains unaffected by the thermal field over the whole simulation. Analogically, it is advised to explicitly model other structural elements adjacent to the analysed one to properly simulate heat exchange between the analysed and neighbouring elements. For other configurations of massive concrete structures, analogous approaches accounting for the underlying terrain can be applied (underlying thicker concrete, previous structural part etc.).

In stress analysis, the underlying terrain can be modelled either explicitly or implicitly. There should be a possibility of a loss of contact between the structure and the soil without transmission of normal tensile stresses at the joint plane. In explicit modelling this can be performed by modification of some of the characteristics of the material model for the soil in the contact elements/interface plane (see Sect. 4 for further details); the material should have no tensile strength and limited possibility to transfer shear stresses. Alternatively, mechanical boundary conditions (supports) can be introduced with specific characteristics, e.g., springs with certain stiffness. It is a common mistake to assume full fixation of the structural element to model external restraint which disables both the possibility of elongation/shortening of the element and its rotation. This has essential effect on the distribution of restraint and resulting stresses in the element [4, 5]. Again, if the analysed element is in direct contact with other neighbouring elements which can influence the stress state by inducing additional external restraint, they should be accounted for by either explicit modelling or introduction of appropriate mechanical boundary conditions.

When cooling measures are taken in massive concrete, those related to pre-cooling of concrete constituents or concrete mixture can be modelled through the initial temperature of concrete. If cooling of concrete is made through embedded cooling pipes, accurate estimation of thermal stress states demands for the explicit consideration of the embedded pipes in the simulation model, which account for the heating effect on the pumped water/fluid along its path. This will significantly increase computational cost but is essential for the estimation of cracking risk due to potentially excessive intensity of cooling capacity, and to account for the thermal gradients between the inner and near-surface parts of concrete. There are, nevertheless, several simplified approaches to the modelling of the effect of cooling pipes in massive concrete, which inherently are less accurate on the regions neighbouring the actual cooling pipes. Even though these detailed matters fall outside the scope of this document, the reader is referred to [7] for a recent review on existing methodologies.

Deep-soil temperature is known to be almost equivalent to the average annual environmental temperature of a given location if one considers depths deeper than, e.g., 4 m [8]. It is noted that such depth limits may be location-specific. On the top part of the terrain (i.e. the 4 m mentioned above), the temperature of the terrain is normally affected by diurnal/seasonal environmental variations. For the above reasons, the following is recommended:

  1. (1)

    If excavation for construction of the massive concrete structure is of 4 m depth (or deeper), and a short period spans between excavation and casting (e.g. a couple of weeks), one could consider the initial temperature of the underlying soil to be constant and equal to the average annual environmental temperature.

  2. (2)

    If the conditions mentioned above are not met, one should carefully entertain the possibility of: (a) explicitly modelling the terrain temperature field before the actual casting for a period of more than 3–4 months; (b) resort to references in the literature with simplified proposed equations for consideration of the underground temperature in shallow regions [9] or use in-situ measurements for the same purpose.

Information regarding the thermal properties of different types of soils and rocks can be obtained from relevant literature [10,11,12].

In massive concrete structures, and for the purpose of cracking risk assessment, it is generally considered unnecessary to actually model the reinforcement whenever: (i) the thermal dilation coefficient of the hardened concrete and reinforcement are similar (e.g. within ± 15% of each other); (ii) autogenous deformations are small (e.g. < 100 μm/m developed up to 7 days after concrete casting). It is noted that the former statement is valid when the reinforcement is placed parallel to the surface or when a minor part is placed perpendicularly to the surface and protrudes outwards. The thermal interactions inherent to this situation have been demonstrated to be low and hence render the simplification of ignoring reinforcement in this concern as valid [13]. In particularly cold-weather concreting situations where also a massive concrete element is densely reinforced at a particular location, the thermal interactions mentioned above may be more pronounced. Furthermore, with particular reference to autogenous deformation, in heavily reinforced sections, even a relatively small autogenous deformation of concrete might cause important stresses due to the restraint to free deformation caused by reinforcement. In such cases, the simulation of the effects of autogenous deformations may become quite relevant.

3 Concrete properties and assumptions

Appropriate selection of concrete properties and constitutive laws that can adequately describe the age/hydration- and time-dependent concrete behaviour are of paramount importance in simulating the thermo-mechanical behaviour of concrete satisfactorily.

Thermal properties of concrete relevant to thermo-mechanical simulation are the (volumetric) specific heat capacity and thermal conductivity. Consideration of the age/hydration and even moisture dependence of these input parameters is not mandatory, and most frequently these are treated as constant with no important penalty on simulation results. Nevertheless, it is recognised that the age/hydration dependence of specific heat can have some effect on the simulation result, as temperatures can be overestimated, and age/hydration dependence may be considered where possible. The maximum temperature overestimation due to the adoption of a constant specific heat, rather than using an age/hydration dependent on may reach 10% whilst the temperature differential overestimation can reach up to 15% [14]. For consideration, a 5% variation in specific heat can lead to a variation of 2–3% in the maximum simulated temperature [15].

The values of thermal conductivity and specific heat capacity of concrete to be used in the numerical simulation normally vary between 1.6–3.0 W/(m·Κ) and 1920–2880 kJ/(m3·Κ), respectively, whilst the latter generally affects the simulation results to a greater extent than the former. A value of 2.4 W/(m·Κ) and 2400 kJ/(m3·Κ) for concrete thermal conductivity and specific heat capacity, respectively, is deemed appropriate for most of the normal-weight concretes and the scope of the simulation. It is recognised that appropriateness of the above-recommended values may be affected by the composition of concrete and the location of collection of the raw materials. The modeller should note that supplementary cementitious materials (SCMs) and aggregate types incorporated, as well as the resultant density (e.g., lightweight concrete) influence the thermal properties of concrete, and their effect should be taken into account where possible/important. Wherever concrete mix-specific thermal properties are required, these may be obtained through relevant testing or analytical/empirical prediction formulas based on concrete composition, as explained in Chapter 3 of the state-of-the-art report of TC 254-CMS [6].

The reinforcement ratio is also expected to have an effect in the combined thermal conductivity of reinforced concrete; however, it is frequently the case that in massive concrete structures reinforcement is situated predominantly in the edges of the elements which has a near-negligible effect on the overall thermal conductivity. Nonetheless, very heavily reinforced elements might bring global changes to the thermal properties to consider for concrete when reinforcement is not explicitly modelled. Indeed, simple homogenization equations such as the ones mentioned in the previous paragraph can be used for such concern.

Consideration of the heat release caused by hydration reactions is considered mandatory in modelling hydration-induced temperature stresses. Different models attempt to calculate the evolution of heat release from cement hydration, and each model has its peculiarities. The most commonly used are outlined in Chapter 2 of the TC 254-CMS state-of-the-art report [16]. In any case, the heat release model should consider the influence of temperature on the kinetics of hydration heat release. This is also known as “thermal activation” of a certain mix’s heat of hydration and is important particularly for mixes containing supplementary cementitious materials, which are known to exhibit different hydration kinetics and temperature sensitivity than Portland cements. This is usually accounted for by using the Arrhenius equation which encompasses a mix-specific “apparent” activation energy and the modeller should use an appropriate value (and compatible with the adopted hydration model) from the literature or from experimental results, which reflects the cementitious binder considered in the simulation. It is a common mistake to disregard thermal activation in the simulation, which can lead to modelling inaccuracies [19]. The modeller should also have confidence that the hydration modelling parameters used originate from proper mix-specific model calibration. A word of caution is given in regard to the collection of activation energy values from the literature. Indeed, different modelling approaches for the hydration process (e.g., the method of Reinhardt [17] vs. Jonasson [18]) can require different values for the property named “activation energy” in a given mix, for the resulting thermal output to be equivalent (e.g. see [19, 20]).

Where possible, it is desirable to collect experimental data specific for the concrete mixes under analysis. If one opts for semi-adiabatic or adiabatic calorimetry, the capacity to infer the temperature sensitivity is limited, particularly when compared to that of isothermal calorimetry. However, when selecting isothermal conduction calorimetry, the sample size needs to be very small (few grams) and care needs to be taken in order to ensure representativity (e.g., same proportions of all constituents in the mix, but always at the expense of the shortcoming of having distinct surface-to-volume ratios from the real concrete), as well as to upscale data from paste to concrete scale. It is normally recommended to resort to isothermal conduction calorimetry with a minimum of 3–4 testing temperatures, as to allow understanding the temperature dependence and activation energy of the binder. It is noted that there could be an occasional, relatively slight, underestimation of heat because of omitting the shear interaction between cement particles and coarse aggregate during mixing when an isothermal calorimeter is used. Hence, the modeller should calibrate hydration models with both isothermal and non-isothermal (semi-adiabatic or adiabatic) results where possible, although the former should be adequate for most cases. Whenever testing is not feasible, one may rely on the growing information existing in the literature for several contexts in several countries (e.g., the same materials from the same origin might have been studied before, and extrapolations/interpolations might be viable).

It is also important to account for the fact that the hydration of cement has already endured some development by t0,analysis as defined in Sect. 1. For this reason, in any modelling effort, the due account needs to be done by considering an initial value of the extent of hydration (e.g. through an initial degree of hydration for the cases where hydration is computed through this variable).

When no cooling/heating measures are taken, the initial concrete temperature, i.e. the temperature of the concrete at the time of casting, is recommended to be taken 3–5 °C above the ambient temperature during casting. This recommendation should be applied with caution and/or appropriate adaptations in situations where ambient temperatures are under 5 °C and above 30 °C.

The evolution of mechanical properties in simulating the hardening behaviour of concrete is of paramount importance and should be adequately considered (information on evolution equations is given in Chapter 4 of the TC 254-CMS state-of-the-art report [21]). To begin with, a constitutive law which takes into account the age/hydration dependence of static elastic modulus should be adopted in the numerical simulation as it directly affects the stress calculation. Such formulations can be found in major design codes and are indeed present in most of the commercially available software for thermo-mechanical simulation of concrete. The user should be aware that these formulations are not usually derived for concretes with SCMs, which exhibit somewhat different evolution of elastic modulus at early ages compared to that of neat Portland cement concrete. Similarly to the evolution of the heat of hydration, the modeller should account for the temperature sensitivity of the cementitious binder used, accordingly. In all cases, if experimental data are available, they should either be used directly in the model or used to calibrate the analytical formulation in the software. In the absence of experimental data, the user should use indicative values from major codes, e.g., 28-day values depending on strength class. Avowedly, the static elastic modulus is greatly influenced by the aggregate type, and its effect should be taken into consideration, where possible.

Coefficient of thermal dilation (CTD) of concrete, also known as the coefficient of thermal expansion, will affect the calculation of thermal strains and is affected in turn predominantly by the aggregate type used. Recommended values usually vary from 8 to 12 μm/m/°C and are location-specific. Research has shown that supplementary cementitious materials do not significantly affect the value of CTD of concrete. The age/hydration-dependence of CTD after concrete setting is very limited (whilst a slight increase in CTD may be occasionally observed) and hence the recommendation is to keep CTD constant in the simulations.

Autogenous (basic) strain is generally significantly lower than the relevant thermal strain in massive concrete structures. The modeller should be aware that the autogenous strain is increasing with the concrete strength, and that for some SCM-concretes the temperature dependence may become considerable. It is also important to note that most test standards and codified models neglect the autogenous deformation occurring between tsetting and 1 day, which in some cases may be significant. Since autogenous strain is relatively simple to take into account in the analysis, it should be considered where possible.

The contraction mechanism of drying shrinkage can normally be disregarded from the simulation due to its limited influence and to its relatively slow progression owing to the geometrical significance of the structures under consideration. As a structure becomes more massive, any effects of drying are diminishing and are only limited to the superficial layers of the concrete. Drying may induce some superficial hairline-like cracking during hardening, however, to even model that, moisture gradients would need to be considered which do not normally constitute a built-in function in commercially available software.

Ageing viscoelasticity, however, should be considered in all cases. It is an important property significantly affecting the resulting stresses throughout the entire time of analysis, e.g., [19, 22]. It is a common mistake to neglect creep in analysis of early-age stresses, usually because of complexity this introduces to the analysis, while it is an important property significantly affecting the resulting stresses throughout the entire time of analysis [19]. It is recommended to use a linear constitutive law which represents at least the effects of age-dependent viscoelastic properties. Further influencing factors on the viscoelastic behaviour arising from increased temperatures during hydration may be regarded in very thick members with comparably long periods of concrete temperature higher than 50 °C for more than 3 days. Furthermore, the most relevant mechanism of viscoelasticity in massive concrete structures is found to be that of basic creep, owing to the absence or minor effects of drying and consequently, drying creep. One must be, however, aware that the relevance of drying outside the scope of the first weeks of age (normally those of importance in many massive concrete structures for thermal cracking) may be large enough for the drying to play a significant role in cracking near the surface regions. Nevertheless, this is considered outside of the scope of the present recommendations.

In a macroscopic simulation of hydration-induced stresses, viscoelastic effects can usually be considered linearly-dependent on stresses (a simplification once high stress/strength ratios are attained), whereas the occurring viscoelastic strains in the course of time shall be explicitly determined for each time step and element, and must also be included as such in the simulation. There are two main approaches for such determination. One option is the formulation on basis of rheological models representing the whole range of creep functions for any loading age implicitly. The other option is an explicit incremental analysis of the entire range of creep functions in combination with the stress history occurring in the analysis time. In general, both approaches give similar results when the calibration of the rheological model represents perfectly the entire range of creep functions in the analysis time. The application of rheological models is widely established in science and practice as it requires a relatively low computational effort and it is simple to apply. However, it shall be noted that rheological models are also subject to certain limitations when applied for the simulation of hydration-induced stress histories with inversion of stresses, e.g., the general assumption of the validity of superposition in phases of unloading and equality of creep in compression and tension. More relevant information can be found in [19, 23, 24]. The application of approaches based on the modification of stiffness, i.e. age-adjusted effective elastic modulus, are usually not adequate for thermo-mechanical timestep-based simulations. Creep is, however, still a matter of significant debate and need for future research, particularly in view of complexities such as: (i) creep of new SCM-concretes; (ii) differences between tensile and compressive creep behaviour; (iii) creep behaviour at very early ages; (iv) behaviour upon sign reversals; (v) thermal activation; (vi) coupling between tensile creep and high stresses with microcracking.

Finally, the age/hydration- (and temperature-) dependent tensile strength of concrete should be considered when part of the desired output is the calculation of cracking indexes, frequently expressed as the ratio of simulated tensile stress to the tensile strength at a particular time instant. The temperature sensitivity (particularly influenced by the binder) should be taken into account in the calculation, if possible. Alternatively, codified approaches may be used with an inherent degree of inaccuracy, particularly if mix proportioning deviates from the conditions considered in the standards, e.g., SCMs. In the case where a particular software does not accommodate calculation of the cracking index, the modeller can calculate the age/hydration dependent tensile strength independently and compare it to the simulation results. The selection of appropriate tensile strength value, e.g., to account for scale/size effects, type of test, etc. and the calculation of a probability of cracking has been discussed in Chapter 8 of the TC 254-CMS state-of-the-art report [25].

4 Boundary conditions and loads

Appropriate definitions of thermal and structural boundary and loading conditions are of paramount importance in obtaining reliable simulation results. In general terms, the modeller should not over-simplify the definitions of boundary and loading conditions as this may result in modelling inaccuracies which may lead to unrealistic evaluation of the likelihood of thermal cracking in massive concrete structures. Conversely, the modeller should not by default over-complicate the numerical model with very meticulous definitions of boundary and loading conditions as this often comes at a computational cost, proneness to modelling errors and potentially marginal benefit to the accuracy of the cracking risk assessment in massive concrete structures.

As the early thermal behaviour of concrete (particularly near the surface) is influenced by ambient temperature, it is important to set this thermal load in a proper manner. Modellers should be knowledgeable of the concept of 'dry-bulb' temperature when looking for data to support predictions. Basic data could be obtained from a weather station nearby the construction site (should be as near as possible, and with similar altitude). Looking for data in the same season of the construction to be performed is recommended. In many countries, data from weather stations are readily available for download and use, making this a very easy pathway to reliable information. Of course, if preferred (or if weather station data are not available), one can seek information in national weather data, particularly in “climatological normals” that can assist in the creation of sinusoidal daily temperature variations representative of the season and location (for more relevant information see [26]). Still in this context, it is important to consider the density of temperature information data to use in the process of modelling. When the entire daily cycle is taken into consideration, it is recommended to consider temperature values at a minimum periodicity of 1 h.

In some cases, where the accurate surface temperature/stresses are not a major demand/concern, it can be acceptable for the ambient temperature to be considered as constant. In that case, it is recommended to use the expected average daily temperature. Even though this simplifies the processes of analysis and interpretation, it has smaller accuracy and does not allow studying the effect of the time of casting (e.g. early morning or late afternoon) in the final outcome in terms of developed stresses (and hence cracking risk). Generally, a decision for accounting for variable or constant temperature would be made based on practice/experience. For non-experienced modellers, the consideration of variable daily temperature is recommended, as it can bring significant impacts in some cases (e.g. countries with significant daily temperature fluctuations, and casting taking place at specific times of day as a countermeasure).

Convective boundaries occur whenever an element involved in the analysis is in direct contact with a fluid. For these recommendations, the scope is limited to air as a fluid for convection. As stated in Sect. 2, convection is taken into account through a convection coefficient that implicitly covers the natural convection phenomena, as well as the forced convection, namely through wind, which is a transient phenomenon by nature, and extremely complex to describe, let alone predict. From an engineering perspective, it is usual to consider a constant convection coefficient for concrete (or other materials) in direct contact with the air. Methods to calculate forced convection are included in Chapter 3 of the state-of-the-art report of TC 254-CMS [6]. Generally, the convection coefficient may vary from 6.0 W/(m2·Κ) for stagnant air conditions to even 25.0 W/(m2·Κ) in very windy conditions. In the absence of wind speed data, the convection coefficient including forced convection can be calculated using any of the aforementioned methods with an expected wind speed on site derived from the anticipated wind intensity based on the Beaufort scale after meteorological observations.

Heat transfer due to radiation is related to energy emission of a body because of its temperature in the form of electromagnetic waves. There are two types of radiation, namely shortwave and longwave, with definitions given in [3]. Heat flux from solar radiation can normally be neglected when analysing thick concrete members, as massive concrete structures. However, there are potential implications that the modeller should be aware of in the case where solar radiant heat flux is ignored, that being wrongly assessing the surface cracking risk of a concrete element. If solar radiant heat flux is to be explicitly accounted for, then models which impose a flux on the surface can be considered, or models which convert the flux to ambient temperature increase [27,28,29,30].

Issues related to structural interfaces in modelling are normally raised when two distinct materials are in contact with each other (e.g. concrete and underlying terrain; concrete and steel part embedded into it), or even when concrete elements are cast adjacently to each other at different times/ages (e.g. wall cast on underlying foundation slab). There is a wide variety of possibilities for the interfaces above (e.g. with or without continuity reinforcement or even structural connectors), or even other situations not contemplated above. Nonetheless, some particular recommendations are stated in the following paragraphs.

When massive concrete is cast on the ground (either rock, soil or treated with a regularisation mortar), a structural interface of complexity can arise, evidenced both because of potential slippage (in the tangential direction, caused by differential thermal volumetric changes), or also by the uplift of concrete elements in extremity regions due to global element curling. For this situation, it is recommended that a structural interface is explicitly modelled between the FE of concrete and the FE of the underlying terrain. Stiffness of the interface can be considered as high as the stiffer material involved in the interface. In the tangential direction, if stresses are expected to induce slippage, a bond-slip law can be considered. Accurate definitions would require specific experimentation. In the absence of better information, one can adopt a limit state analysis approach, e.g., through comparing results from a calculation which considers full bond, and those resulting from a calculation that disregards bond. As such, if the discrepancies are significant, the modeller could look into more specialised bond-slip models. In what concerns the behaviour of the interface in the normal direction, it is recommended to consider a perfect linear elastic behaviour in compression, whereas tensile capacity should be neglected.

When a massive concrete structure/element is cast in a sequential manner, causing what is frequently labelled as 'cold joints', the corresponding interfaces are not normally explicitly modelled as such, and continuity of the material is assumed in between the concrete elements. This is generally considered as an acceptable simplification because of the usual continuity role of perpendicular reinforcement through these interfaces, the normally induced roughness at these surfaces (as a good construction practice), and also by the fact that these are not commonly regions of concentration of cracks in massive concrete structures/elements.

Mechanical support conditions can have a decisive influence on the realism of stress/cracking predictions in massive concrete structures. When the terrain underlying the massive concrete element is explicitly modelled, mechanical supports are only normally considered in the bottom surface of the terrain, with all degrees of freedom prescribed. If the underlying terrain is not explicitly modelled, support conditions are applied directly to the massive concrete. In this case, the choice of degrees of freedom involved, stiffness of supports, or even the potential non-linear behaviour (e.g. support uplift) need to be considered judiciously. Naturally, when symmetry is considered in the model, the necessary degrees of freedom need to be prescribed accordingly.

Structural loading is inherently involved in any thermo-mechanical simulation, at least with regards to the self-weight of concrete. One must, however, bear in mind that modelling of massive concrete structures/elements with applied self-weight since the beginning of the analysis can bring about excessive deformations, as the very fresh concrete is kept in place by the forms, which are not normally explicitly modelled in the mechanical simulation. Therefore, the activation of self-weight in the mechanical simulation is only generally made when the stiffness of concrete has reached a value high enough to mitigate erroneous strain estimations. Alternatively, one might also activate self-weight right since the beginning of simulations, provided that a large enough initial stiffness has been considered [31].

Other structural loads such as dead loads (other than self-weight) or live loads, or others (e.g. wind, snow), can be considered with no specific limitations in the context of thermo-mechanical simulation of massive concrete structures.

5 Calculations/procedures

The choice of FE type and size in thermo-mechanical modelling follows the same principles used in general scope applications of FEM. As gradients of temperature/strain are steeper in the vicinity of external boundaries, it is not unusual to see densification of FE meshes in such regions. When cooling pipes are explicitly modelled, the FE mesh of concrete around them needs to be much denser than at any other region because of very significant temperature gradients normally occurring around the cooling pipes. The principles of selection of mesh sizes rely on the intent to reach mesh convergence (i.e. using a mesh size that is dense enough as not to cause errors induced by insufficient mesh density).

Choices of FE mesh sizes for massive concrete end up depending significantly on the size of the model itself (with strong differences between a 300 m tall whole arch dam model and a mere 20 m long spillway wall with 2.5 m thickness). Very relevant hurdles can be brought about by limited computational capacity and limited storage capacity. When relatively small models are considered, desired mesh sizes can have typical edges with 20–50 cm (for the case of linearly interpolated elements in thermal analysis), whereas very large models are likely to include FE edges as large as 1–2 m (particularly in core regions). It is also usual to adopt denser meshes in the near-surface regions of massive concrete because of the expected steeper gradients of temperature (and hence stresses) in these regions at early ages, when compared to core areas (which can tend to have coarser meshes). Despite this, and taking into account the increasing computational power available for the task, it is not infrequent to see relatively uniformly sized meshes matching the element sizes one would expect in the regions of higher gradients.

The time steps in a thermo-mechanical analysis are not likely to be uniform, and they depend on several factors that are described next, with specific recommendations. Regardless of the type of hydration heat model used, there are important volumetric heat releases at early ages in the process of modelling massive concrete structures. For proper accuracy in the modelling (particularly because of the non-linear processes caused by the explicit consideration of thermal activation), it is recommended to consider time steps of 15 min and up to 60 min, with a preference for the former in favour of accuracy. These relatively small step sizes should be kept for as long as the exothermic heat release is relevant. In thick walls, this may be as long as 7–10 days, whereas in dams, this can take several months. Another important aspect to take into consideration when defining the time step sizes is the degree of detail in consideration of thermal loads such as solar radiation of environmental temperature throughout the day.

In order to discuss convergence algorithms and criteria, one has to discuss separately the thermal and the mechanical calculations, with reference to the sources of nonlinearity. The nonlinearity induced by the thermally activated nature of the heat of hydration is normally mild. Standard Newton–Raphson strategies are generally enough to handle this matter. Energy-based convergence checks with a relative tolerance of 0.001 or smaller are recommended. More complex and sophisticated algorithms can be necessary for the application of cooling pipes, as mentioned above.

It is noted that stress averaging in nodes may potentially bring relatively high inaccuracies. This has been traditionally encountered in models involving coarsely-meshed concrete blocks on soil. After cooling, tensile stresses develop on the concrete surface, which are averaged with low soil stresses. Another example where this issue is found is when a block is cast on a concrete layer. Ignoring layer’s compliance (fixed vertical displacement) or debonding between block and the layer may lead to computing different results.

As the present document does not include the modelling of cracking, the only source here considered for stiffness matrix nonlinearity is related to the interface elements. When fracture behaviour is considered, the recommended practice follows the typical procedures already applied in standard mechanical analyses. More theoretical information on numerical modelling and relevant constitutive laws can be found in Chapter 7 of the TC 254-CMS state-of-the-art report [22].