D1. Civil engineering materials

D1. Civil engineering materials

The objective of this chapter is to highlight the peculiarities of civil engineering materials regarding FE modeling. Among those, there are the particularities of the civil engineering materials themselves (concrete, timber, steel), and their use within the structure. 

1. Concrete

Concrete is a compound material composed of cement paste, and aggregates of various sizes. The cement paste that contains water and additives acts as a binding agent (glue) between all the aggregates. Even though concrete is a heterogeneous material, it revolutionized the world of construction because of its characteristics such as its versatility, adaptability (offering a variety of aspects: form, color, texture), mechanical resistivity, fire resistivity, acoustic insulation, and diversity (high-performance concrete, precast concrete, etc.). However, concrete is subjected to degradation processes. The latter can be due to environmental factors such as humidity, rain, cold, freeze-thaw, heat, wind, drought, aggressions from deicing rock salt, alkali-reactions, internal or external sulfate-reactions, corrosion, mechanical solicitations due to collisions, and load increase. Depending on the degradation phenomenon, the deterioration of the structural element is progressive and initiates as soon as the concrete tensile stresses are applied. Degradation leads to the emergence of small cracks and localized deformations. The increase of the deformation leads to the decrease of the material’s strength, which can then lead to softening behavior and the collapse of the structure. 

Constitutive law – In the FE calculations, generally, only some aspects of the concrete behavior are considered because of its complexity. To express the most relevant constitutive law, it is necessary to well understand the characteristics of concrete. Thus, it is essential to know the behavior of the cement paste, as well as the behavior of the aggregates, both responsible for the stiffness properties of the material. The paste is the cause of the mechanical strength of the concrete and, because of its water content, is also the cause of the deficiencies such as the increase of porosity, the decrease in mechanical strength, the delayed effects (shrinkage, creep), and the transfer of aggressive agents. 

Concrete is considered as a “quasi-fragile” material because of its complex behavior. Its main characteristics are highlighted in the literature. The most significant one is its fragile behavior in tension and more ductile in compression and when subjected to unilateral forces in a tension/compression cycle. Other properties must be considered such as the shrinkage and creep according to the environment and the loading applied to the concrete. 

In the literature, the two main approaches to model the behavior of concrete are:

• Models relying on a continuous approach, meaning that concrete is considered as a continuous media (Bazant, 1979). Using this approach, the calculations of the cracks are deduced from the displacement-force relationship. There exist different models such as the elasto-plastic (Ottosen, Drucker-Prager), damage (Mazars, Laborderie), smeared crack model, of gradient type, or energy type models (fictitious crack model (Hillerborg, 1984), crack band model). In some cases, elasto-plastic models are used, but the user must be very careful, once the elastic phase is overcome, it is possible to obtain deformations that do not reflect the behavior of concrete at all. 
  
  

• Models based on a discrete approach. Indeed, concrete cracks generate geometrical discontinuities, which are integrated between the boundaries of the connected elements. There exist discrete models of Ngo, (Ngo, 1967), Bazant (Bazant, 1979), Blaauwendraad (Blaauwendraad, 1999), and other lattice-type models.

However, in the literature, the two most used model categories are the ones based on mechanical damage and the ones considering explicit discontinuities. 

Although the elastic component is not exactly linear ([Baron] p.276; [Sargin 1971]), the calculations are generally conducted using the isotropic linear elastic assumption (Hooke’s law) to represent the elastic phase of the concrete behavior. The values of Young’s modulus and Poisson’s ratio are computed by characterization of the concrete. In the eventuality that other mechanical characterizations of the concrete cannot be conducted, some properties such as the compressive strength can be deducted using model laws that can be found in model-codes such as Eurocode 2 for instance. 

Note – Design offices use the model-codes except for study cases when they conduct detailed analyses on a specific structural element on which particular mechanical phenomena occur: cracks, heave, or local heterogeneity. Indeed, there exist finite elements of “cracked element” type, damage models, viscoplastic laws, or even poro-mechanical coupling laws. However, to conduct these studies, data from on-site or sample analyses of the material must be available to orient the engineer towards adequate constitutive laws, for which, coefficients and moduli will have to be determined anyway.  

Elastoplastic model – For the sake of simplicity, concrete is often modeled as an isotropic homogeneous elastoplastic material, which is acceptable until the concrete experiences cracking. 

Regarding an elastoplastic law, concrete cracking is not directly modeled: the cracking zones are characterized by important anelastic deformations (>1-2%) and a fixed constraint state between tensile f_t and compressive strength f_c of concrete. 

As shown in chapter 3, an elastoplastic constitutive law for concrete is composed of an elastic law and a plastic flow rule associated with a “criterion of plasticity”.

However, the post-cracking behavior can be approximated using a stress-strain curve going further than the tensile strength. The latter can cover the effects of tensile softening (representing the necessary work to open the crack) and tensile hardening (concrete contribution in between the cracks, stresses due to reinforcement adherence).

If the decreasing behavior of the material leads to a global decreasing behavior, one must pay attention to the eventual local effects: the size of the finite elements will limit the size of the anelastic zones and the solution will depend on the used mesh. Different numerical techniques enable us to solve or limit this issue. 

Damage – A damage law is a law that considers one of the main macroscopic effects of concrete cracking: the loss of material stiffness. The fundamental idea is to renounce following the eventual cracks (their apparition and propagation) and consider that the concrete of a given structure deteriorates by multiplication of cracks in the damaged areas. This type of law allows us to describe the decrease in material stiffness undergoing small cracks. This stiffness loss is measured by an internal variable called damage, denoted D, that evolves from 0 (undamaged material) to 1 (totally damaged material). This variable is generally a scalar value. 

To best represent the behavior of concrete, the damage laws consider post-peak softening behavior. This enables determining the stresses as follows:

with  

The advantage of this method is to consider the concrete as a “continuous medium”, for which the FE methods are perfectly adapted. 

In the case of concrete behavior, the two main damage modeling families are the anisotropic and isotropic ones. The isotropy characterizes the invariance of the physical properties of concrete regardless of the direction. On the other hand, the anisotropy depends on the direction. An anisotropic law presents different responses due to loading according to its orientation. 

One of the most used damage models in the industrial and research world is the Mazars model [Mazars, 1984]. It is certainly the first damage model for concrete that works accurately. 

The main difficulties raised by the damage models are:

• a dependence a priori of the meshing results: note that one should, in principle, demonstrate that the ruin mechanism obtained with this type of model is independent of the refinement of the mesh, at least after a certain threshold. Besides, this dependence led to the development of the regularization method. 

• the absence of an analytical solution in these simple cases. 

Several regularization classes exist including the non-local regularization and the regularization by cracking energy [Hillerborg, 1976] (that solves the problem only partially). Among those non-local methods, it is possible to cite the integral methods [Pijaudier et al., 1987], [Giry et al., 2011] or gradient methods (deformation gradient or internal variable gradient [De Borst et al., 1992], [Peerlings et al., 1996], [Nedjar, 2005], [Lorentz, 2017]. These methods require the use of relatively thin meshes, which renders the analysis computationally expensive. 

Delayed effects of the stress relocation – When analyzing the behavior of a concrete structure several weeks after pouring, but especially for the long term, it is necessary to consider the delayed effects such as shrinkage and creep. 

These phenomena specific to concrete can, in principle, be modeled adopting a visco-elasto-plastic constitutive law (Bingham law), or sometimes only visco-elastic (“scientific” creep, [Eymard]): this approach is generally applied by research laboratories to analyze tests on materials. However, in the case of a refined model considering the delayed effects, it is necessary to incorporate phenomena such as drying and cracking in addition to shrinkage and creep. Indeed, all these phenomena interact with one another, so it is important to model them numerically. 

The deformation caused by the shrinkage is induced by the drying of the concrete due to the environmental effects. The shrinkage leads to a differential deformation, meaning that the stresses are more significant where the concrete dries faster. This phenomenon causes tensile stresses at the surface, which then leads to cracks and compressive stresses inside the body. 

Concerning creep, its deformation is generally separated into two deformations, one is due to creep itself, and the other one to shrinkage from the drying process. The latter can be explained by the strong dependence between relative humidity and creep. 

In Eurocode 2, it is possible to determine concrete deformations due to the delayed effects (without external loading). For this purpose, one must compute the deformation due to the endogenous shrinkage (caused by the internal humidity) and the desiccation deformation (caused by the drying process and the size of the structural element). According to the §3.1.9 of the EN1992-1-1, the creep deformation under compressive stresses σc can be written as:

with E_c the tangent modulus (equal to 1.05 E_cm) and φ the creep coefficient. 

The Eurocode 2 (EN 1992-1-1) describes a simplified calculation method of φ (∞, t0). Annex B gives a more complete method enabling the estimation of φ(t, t₀) as well as the shrinkage evolution (see the annex B of EN 1992-2). Moreover, it is important to remember that the calculation of the delayed deformation is relative to the type of cement paste. 

2. Structural steel

The behavior of steel is much simpler than the behavior of concrete for several reasons: it is an isotropic material with identical strength and moduli in tension and compression. Moreover, it is subjected to industrial quality control to ensure its homogeneity. 

Elastic models – Even though steel behaves essentially as anisotropic elastoplastic material (Eurocode 3, part 1.1, §5.4.3), the typical models for steel structures, or compound containing steel, adopt a linear elastic behavior. One must then check that the elastic constraint threshold fy was not reached. In the case of beam or column type elements, the codes allow, if the sections have a size that guarantees a sufficient local ductility, to exceed in the analysis the elastic strength and consider reaction moments based on the plastic distribution of the stresses. 

Elastoplasticity and strain hardening – The elastoplastic theory was developed from the study of steel alloys, mainly to predict the rolling and forging stresses ([Hill], [Nadai]). For steel, the usual model is composed of Hooke’s law for the elastic deformation and the plasticity criterion denoted “Von Mises criterion”, as explained in the Eurocode 3 (part 2, §7.3; part 1.5 §10 and part 1.7 §5.2.3.2): for this criterion, only the elastic constraint threshold fy of steel must be provided. 

Considering the FE calculations, another issue is the strain hardening of the material, meaning the hardening in the purely plastic phase. This aspect is notably explained in the annex C6, part 1.5 of the Eurocode 3. Since it is hard to verify the model, it is preferable to conduct a preliminary computational model with no strain hardening: indeed, one can then verify the quality of the obtained results in the plastic zones by observing the eigenvalues of the Von Mises stresses (the plastic zones must be about monochrome). Nevertheless, it is important to note that the use of an elastic perfectly plastic law with no strain hardening can lead to convergence issues of the linear analysis. Indeed, the plastic zones have zero tangent modulus and, therefore, no stiffness. Considering strain hardening can stabilize the numerical resolution. 

In the case of strain hardening, until a stress fy+X is reached, the possible loading/unloading phases are elastic, with a modulus equal to the initial one. Thus, fy+X becomes the new elastic yield stress of the material. Moreover, the behavior of the steel is close to the kinematic strain hardening model. As a first approach, one can consider a tensile strain hardening up to a certain stress value (fy+X) that will decrease the compressive elastic yield strength to fy-X, and vice versa. 

The most common steel model is the isotropic strain hardening, which depends on the accumulated plastic deformation: 

where  and f_y an increasing function.

In this case, there is no distinction between tension and compression, so a strain hardening in one direction leads to an increase of the compressive elastic yield stress. This hypothesis is accurate only if the considered loading is monotone and not cyclic. 

In most numerical computational models, one can retrieve the eigenvalues of the cumulated plastic deformation. Thus, it is possible to verify that the Von Mises constraint is distributed according to this deformation in the plastic zones. 

In the specific case of stainless steel, the hypothesis of a linear behavior up the limit of elasticity fy corresponding to plastic deformation of 0.2% is not respected, so it is recommended by the Eurocode 3-1-4, § 4.2. to consider the effects of the nonlinear behavior on the force-displacement relationship in the deflection calculations. The material law is of Ramberg-Osgood type and recommendations to implement this model can be found in annex C of this section of the Eurocode. 

3. Pre-stressing steel

Regarding their behavior, steel cables (braces, tension cables) and steel for pre-stressing differ from other structural steel only by an elastic limit about three times greater (from 1680 to 2140MPa in France), which is justified by the serviceability solicitations, generally greater than 1000MPa. One can then, firstly, model the pre-stressing reinforcement as linear elastic elements or elastoplastic of Von Mises type. Steel elements for pre-stressing are all hot-rolled, so the linear isotropic strain hardening law in the plastic phase should be taken into account. The strain hardening coefficient h can be calculated using the characteristics guaranteed by the manufacturer:

with R_m the constraint to failure, f_p0,1 the usual elastic limit, E the Young modulus, and A_gt the elongation at the point of failure. 

A direct consequence of this high level of stress is the initiation of the relaxation mechanism. Relaxation of steel is a non-elastic delayed effect that depends essentially on the time spent since the loading was applied. It leads to a progressive decrease of the constraint for bars and cables subjected to a constant elongation. This mechanism initiates at ambient temperature only for elongations exceeding about 60% of the elastic limit (≈1000MPa). Relaxation increases slightly with temperature. 

Alike concrete creep, it is possible to model the relaxation using a linear visco-elastic constitutive law. However, this approach is mostly used in research work. It is not common to model steel for pre-stressing: generally, it is accounted for by introducing a distributed force within the pre-stressed concrete. Nevertheless, for detailed analysis, in which the interaction between the steel cables and the surrounding medium is of interest, relaxation can be considered using an incremental time-dependent calculation with the relaxation losses introduced as initial stress conditions. 

4. Passive steel

Normal mechanical actions of passive steels – For FE calculations, passive steel elements are frequently modeled by linear 1D elements of “bar” type. According to the historical reinforced concrete methods, the reinforcement is considered to work only in tension/compression along its own axis, which is exactly what FE models enable. Connecting those bar elements to the nodes in the solid elements representing the concrete makes it possible to avoid mechanisms specific to “bar-only” assemblies.

In most computational software, the bar elements are considered by default as elements with a linear elastic behavior. Considering the 1D characteristic of these elements, the law is of the type N = E A u/L, with u the axial displacement, and N the normal stress at the nodes. It must be verified afterward that the axial force N/A stays within ±f_y.

Friction – What was said in the above paragraphs is valid only if there is a continuous friction force acting between the steel reinforcement and the concrete. 

5. Timber

Timber is a material that possesses some inherent peculiarities: it is not homogeneous (this is true at various scales), it is not isotropic, it does not show symmetrical behaviors, and it can be subjected to ductile or fragile failures according to the solicitation and its orientation. Timber is sensitive to humidity, which has a direct impact on its dimensional characteristics, stiffness, and strength. The duration of the loading has an important impact on the strength and the deformation of the wooden elements. The variations in humidity can amplify the deformations (mechano-sorption). Relying on the previous observation, one can imagine a method to implement the FE analysis in this particular case. 

Homogenization – The hypothesis must be validated for the Representative Elementary Volume (REV), or at least, the volume of the smallest finite element. Knowing that according to the different timber species the growth ring can exceed 1cm, it becomes difficult to assume this homogenization hypothesis close to the connections (peaks, pins, bolts…) with a diameter of the same order of magnitude, or smaller. 

The presence of nodes is rarely considered when modeling structures or structural elements. 

Orthotropy –Timber has a structure and characteristics that are a function of three directions, the lengthwise directions (the axis of the trunk, wood fibers), the radial directions, and the tangential directions perpendicular to the lengthwise direction. The last two constitute a plan, often named the cross-sectional plane of the beams, on which the growth rings can be visible. The orthotropy coordinate system is therefore quasi-cylindrical. On the contrary, the coordinate system in which the elements are inscribed is cartesian. The representation of this orthotropy, when considered, is in most cases limited by an isotropic transverse hypothesis (radial and tangential axes having identical characteristics) in a cartesian frame. The slope of fibers associated with the presence of nodes is generally not considered in the FE calculations. Rather, it is modeled in the assembly calculation especially for the assemblies by contact.  

The flexibility matrix (inverse of the stiffness matrix) can be defined as follows for an orthotropic plan:

Note: the orthotropy and symmetry hypotheses of Sij reduce the number of independent terms from 36 (most general 3D case) to 9 terms. 

Timber – orthotropy hypothesis

Elastic modeling – One must simply know the isotropic transverse matrix of behavior for a 2D or 3D model. The moduli between the lengthwise/longitudinal, radial, and tangential directions can present ratios of the order of 20. They evolve as a function of the loading period (creep), timber humidity when constructing, and during its lifetime (environment). The representation of the FE method will depend on the pertinence of the parameters considered. 

Plastic modeling, failure criterion – For an isotropic transverse non-symmetric material, the Hill, Tsai… criterion can be used knowing that it will be necessary to describe the perpendicular tensile and shear fragile failures. The significant strength variability makes it harder to fix the criteria parameters. The sequence of growth rings or assembled layers with different mechanical characteristics can render a homogeneous model difficult to fix in terms of strength.  Indeed, the strength and stiffness are highly correlated for timbers and “systemic” effects appear quickly in terms of element strength. Consequently, there are strength limits in axial tension, axial compression, bending….

Finally, timber structures are particularly sensitive to the behavior of their assembly. Thus, the elements do not need to be modeled in detail, only the assembly zones do. Nevertheless, problems can arise because of the contact zones, the multitude of materials, some of them entering plasticity, the fragile failure of others, and the homogenization limits reminded earlier. 

It is clear that the modeling effort is related to the scale of the investigation or the advancement of the project. Provided the above singularities are considered, the modeling of timber inside of a structure can be conducted like any other material. 

Modeling of the delayed effects and the interaction with the hydric/water phenomena – Timber is a hygroscopic material sensitive to the relative variation of humidity of the air. Moreover, it is subjected to creep under loading. If one wants to model these phenomena, a rheological viscoelastic model in a variable environment and respecting the thermodynamic principles can be adopted. The generalized viscoelastic behavior of Kelvin-Voigt type, characterized by rheologic parameters and dependent on the level and history of humidity can be associated with the Ranta-Maunus non-viscoelastic behavior to characterize the shrinkage-swelling and the mechano-sorption.