Structural Modelling at the Micro-, Meso-, and NanoscalesView this Special Issue
Fire Spalling Prevention via Polypropylene Fibres: A Meso- and Macroscale Approach
A deep understanding of concrete at the mesoscale level is essential for a better comprehension of several concrete phenomena, such as creep, damage, and spalling. The latter one specifically corresponds to the separation of pieces of concrete from the surface of a structural element when it is exposed to high and rapidly rising temperatures; for this phenomenon a mesoscopic approach is fundamental since aggregates performance and their thermal properties play a crucial role. To reduce the risk of spalling of a concrete material under fire condition, the inclusion of a low dosage of polypropylene fibres in the mix design of concrete is largely recognized. PP fibres in fact evaporate above certain temperatures, thus increasing the porosity and reducing the internal pressure in the material by an increase of the voids connectivity in the cement paste. In this work, the contribution of polypropylene fibres on concrete behaviour, if subjected to elevated thermal ranges, has been numerically investigated thanks to a coupled hygrothermomechanical finite element formulation. Numerical analyses at the macro- and mesoscale levels have been performed.
In the last three decades, catastrophic fire events in concrete tunnels (such as the Danish Great Belt one, the Channel Tunnel, Mont Blanc, and Tauern tunnels) led to the development of new fire protection systems to increase the safety of the people as well as the strength of concrete structures under high temperatures.
One of the principal problems in concrete under fire exposure is spalling, corresponding to the ablation of concrete segments until collapse of the structure, when exposed to high and rapidly rising temperatures. This phenomenon is influenced by many factors but can be explained by taking into account two main contributions: a thermal stress, generated by a thermal gradient between the heated surface and the internal concrete zone, and the pore pressure increments in concrete that occurs when the internal water evaporates .
In order to prevent spalling, when concrete structures are subjected to elevated temperatures, the role of polypropylene (PP) fibres into concrete mix design is largely recognized. These fibres are anisotropic monofilaments (diameter μm and length mm) that are not able to increase the material stiffness, but under high temperatures the explosive spalling risk decreases if the volume content in the mix design is between 1 and 3 kg/m3 of concrete.
Studies have proved that PP fibres reduce pore pressure into the cement paste, so decreasing the risk of explosive spalling. In fact, the fibres are in a solid state at room temperature but when it increases the state changes; melting phase starts at about 165°C and over 325°C the fibres vaporize and the connected voids channels increase the permeability and diffusivity [1, 2] into concrete. PP evaporation process terminates when temperature is over 475°C.
In this work, the complex mechanism of polypropylene contribution on concrete behaviour under thermal conditions will be numerically investigated through a 3D thermohygromechanical finite element code [3, 4], appropriately updated to take into account the effect of the polypropylene fibres if they are added in the mix design . Innovative concrete systems are investigated in [6, 7].
Numerical analyses at the macro- and mesoscale levels have been performed and validated considering experimental tests by literature. Recent theoretical and computational advances about composite structures also at nanoscale can be found in [8, 9].
To simulate spalling phenomenon, it is not possible to use a linear constitutive law of the material, as concrete has a brittle behaviour; numerically, the softening branch of the material can be described through several theories such as the fracture [10, 11] and damage [12, 13]. In this work, Mazars’ damage law  with nonlocal correction has been adopted.
A mesoscopic approach has a remarkable importance for understanding specific concrete phenomena, such as spalling. Indeed, this representation is able to determine the effects of internal hyperstaticity due to the different mechanical characteristics, triggering stress concentrations that can lead to damage.
2. Concrete as a Multiphase Material
Although traditional engineering studies consider concrete as a homogeneous material, idealized as a continuum medium with average properties (macroscopic approach), concrete is a highly heterogeneous material and its composite behaviour is exceedingly complex.
On a macroscopic approach, most of the works proposed in literature assume phenomenological relationships based on macroscopic observations; even if this approach implies a series of simplifications, using continuum-type constitutive models, a satisfactory description of the basic features of the mechanical behaviour of concrete has been reached. Anyway, to obtain a deeper understanding of the macroscopic constitutive behaviour of concrete, it is necessary to adopt a lower scale of observation, that is, mesoscale.
A mesoscale approach will provide a more realistic description of concrete than the macroscale, influenced by the geometry and the properties of its multiple constituents. This could be expected, since the observed macroscopic behaviour is a direct consequence of the phenomena, which take place at the level of the material heterogeneities.
At this level, concrete becomes a mixture of cement paste with aggregates inclusions of various sizes. Aggregates generally occupy 60–80% of the volume of concrete and greatly influence its properties, mix proportion, and economy. Aggregates can be divided into two distinct categories: fine (often called sand) and coarse aggregates; the latter represents around 40–50% of concrete volume. However concrete is not just a two-phase composite; it has been found that the presence of grains in the cement paste causes a thin layer of matrix material surrounding each inclusion to be more porous than the bulk of the surrounding cement paste matrix. This layer is named interfacial transition zone (ITZ) and has relevant effects on the properties of concrete, being likely to act as the “weak link in the chain” when compared to the bulk cement paste and the aggregate particles [15, 16].
For the numerical simulations at the mesolevel mesoscopic continuum models have been adopted; each single composite constituent itself has been approached as a multiphase material, fully described and characterized via a coupled thermohygromechanical model.
Coarse aggregates have been simplified assuming spherical shapes, in order to eliminate possible stress concentrations generated by the angularities; they are distributed randomly in the concrete sample and have an elastic behaviour (they do not creep and do not damage).
The mortar matrix, comprehensive of the cement paste and of the fine aggregates, and the ITZ, whose thickness is strictly related to the diameter of each aggregate, are homogeneous materials; they can be subjected to creep and damage.
Finally PP fibres, having a size on the order of micrometer, are not explicitly represented in the mesoscale approach; their presence and effect have been taken into consideration updating the concrete formulation as explained in the subsequent section.
3. Theoretical Background
Concrete is considered as a multiphase system where the voids of the skeleton are partly filled with liquid and partly with a gas phase. The liquid phase consists of bound water and capillary water, while the gas phase, that is, moist air, is a mixture of dry air and water vapor and it is assumed to behave as an ideal gas.
When higher than standard temperatures are taken into account, several phenomena are considered within the code, dealing with concrete as a porous medium: heat conduction, vapor diffusion, and liquid water flow in the voids.
As regards the mechanical field, the model couples shrink, creep, and damage within the constitutive law of the material. For details, please refer to [17, 18].
In order to take into account the effect of PP fibres, concrete porosity formulation has been enriched and the microcracking that appears after PP fibres vaporization around the void channels has been considered .
3.1. PP Fibres Effect on Concrete Porosity
Total porosity at different temperatures and different dosages of fibres can be expressed as where the capillary porosity and the micro- and macroporosities and terms depend only on the mix design, while the rest of the equation is related to the PP fibres effect at different temperatures. Specifically, the porosity terms dependent on fibres mean the following: is related to aggregates ITZ connected by fibres; is due to ITZ formation around fibres; is related to voids variation due to PP fibres relaxation; takes into account micro air bubbles connection and consequent porosity increment; considers void channels formation, at high temperatures, after fibres evaporation; and finally is the porosity increment to the microthermal cracks.
Appling this porosity formulation to the hygrothermomechanical FE code, the hydraulic diffusivity will vary accordingly.
Following Bazant’s formulation , is so expressed: where is the diffusivity part, dependent on the physics characteristics related to the ratio.
The PP fibres employed in the concrete material at elevated temperatures, increasing the porosity and the connections between the material voids and not affecting the chemical reactions in the matrix, can be taken into account modifying only the term : where and are constitutive parameters.
3.2. Microcracking Related to PP Fibres
An aspect that has to be considered, when determining the concrete diffusivity containing PP fibres, is the microcracking that appears around the void channels once the fibres evaporate . The cracks, being a mechanical response of the cement paste due to a local stiffness variation, cause a void increment, so reducing the internal concrete pressure after the peak and causing a variation in hydraulic diffusivity.
For the determination of the microthermal cracks, a formulation that considers a damage variable dependent only on temperature has been developed :where and are the minimum and maximum temperatures when microcracks occur and is the initial thermal expansion. The term is the linear variation of the thermal expansion during heating and is the thermal variation when microcrack takes place:where is the tangent modulus of the damaged curve .
The hygrometric diffusivity has been modified to take into account the microcrack contribution as follows:where can be found experimentally.
3.3. Visco-Elasto-Damaged Formulation
The skeleton of the concrete material has been represented as a viscoelastic material coupled with Mazars’ damage formulation .
Following a FE formulation the stress-strain relation can be so expressed:where is the analysis time, is the relaxation time, , are the stress and strain tensors, respectively, is the scalar damage parameter, B is the derivative shape function operator, and is the relaxation function, dual to the compliance or creep function, in accordance with the Maxwell-Chain model:where represents the th elastic modulus of the Maxwell unit and y is the reduction time parameter.
The damage variable d is evaluated following Mazars’ formulation , where the load function is assumed as follows:where k is an internal variable and is the equivalent strain:The damage law is an exponential function, subdivided in a compression and in a traction part:where , (with ), and are the material parameters; are the coupled coefficients.
4. Numerical Models
4.1. 3D Numerical Analyses at the Macro- and Mesolevel without PP Fibres Inclusions
Before investigating the effect of polypropylene fibres, the behaviour of a base cell of concrete, without addition of PP fibres, modelled at the macro- and mesolevel and subjected only to a heating rate has been presented in this subsection.
A representation of the concrete at the mesolevel will allow us to determine the effects of internal hyperstaticities due to the different characteristics of the individual phases, which are due to concentrations of stresses and strains during the heating process of the material and that can lead to a mechanical damage. Such concentrations are not visible by using macroscopic models.
Details of the cell adopted in the analyses are visible in Figure 1; symmetry conditions have been assumed to simplify the numerical model.
A cylindrical cell of radius equal to 10 mm has been considered. As regards the cell at the mesolevel, it contains approximately a 40% of coarse aggregates with an average diameter equal to 1.2 cm, assumed spherical in order to eliminate any possible stress concentration due to angularities, and a thin layer of ITZ whose thickness is closely related to the diameter of the inclusions. The discretization adopted here consists of 1402 nodes and 6708 linear tetrahedral elements.
The hygrothermomechanical properties of the various components have been properly assumed based on literature data as follows.
Material Parameters Matrix [MPa]: ) elastic modulus. : () Poisson’s coefficient. : () material density. °C)] () thermal conductivity. °C)]: () specific heat. : () coeff. Bazant formulation. : () exponent in diffusivity equation (Bazant formulation). : () initial diffusivity. [°C]: () room temperature. Age [day]: () concrete age. : () initial porosity in the concrete. : () aggregate diameter (average). Aggregates : () elastic modulus. : () Poisson’s coefficient. : () material density. °C)] () thermal conductivity. °C)]: () specific heat.In the performed analysis, the sample was subjected to a linear thermal ramp applied on the outer surface, which starts from a room temperature of 20°C, up to 300°C in 1000 sec, and then is kept constant.
Mechanical damage triggering for the macro- and mesoscale models is shown in Figure 2.
As can be noted, while at the macroscale damage interests at the beginning the external ring and then enters in the sample, at the mesolevel its evolution is driven by aggregates distribution and ITZ. Indeed, damage initially triggers within ITZ (being the weakest zone due to its properties), and then it will spread inside this layer and finally will involve the cement paste. Damage concentrations are located where the distance between two different aggregates are minimal and the stress concentrations are maximal.
In the next section, a model including fibres will be presented and damage evolution comparisons will be discussed.
4.2. Numerical Analyses Considering the Effect of PP Fibres
4.2.1. Macroscopic Numerical Models
At first, the numerical model has been validated at the macroscale level referring to the experimental tests of Phan , where a concrete sample of size mm has been heated with a thermal ramp of 5°C/min, reaching a maximum temperature of 600°C. The concrete sample has been coated with insulating material on all sides except the face that will be subjected to heating, to reproduce the behaviour that would have a wall if subjected to a fire. A dosage of 1.5 kg/m3 of PP fibres has been added to the mix design (fibres length equal to 13 mm, with a diameter of 0.1 mm).
Materials parameters for concrete and PP fibres are summarized as follows.
Properties of Concrete and PP Fibres Concrete [MPa]: () elastic modulus. : () Poisson’s coefficient. : () material density. [W/(mm°C)]: () thermal conductivity. [J/(T°C)]: () specific heat. [°C]: () room temperature. Age [day]: () concrete age. PP Fibres : () volume. : () mean diameter. : () mean length. ITZ [mm]: () ITZ thickness. : () relaxation coefficient of PP fibres.The specimen did not show spalling phenomenon. Pore pressure data, measured at 25 and 50 mm from the heated surface, were collected; the temperature was also measured, in addition to the data of the heated surface and the temperature of the furnace.
The results reported in Figure 3 show a good agreement between the experimental and numerical results in terms of pore pressure as a function of time at a distance of 25 mm from the heat source. An abrupt drop of the pressure after the peak can be noticed, due to the fusion of the fibres and to microcracking, which increase the porosity and consequently the diffusivity of water vapor in the specimen.
This set of macrolevel analyses did not capture the real local behaviour related to the presence of aggregates, that is, stress concentrations, but they did correctly assess the mean value of the pore pressure into the sample at different time steps, although the obtained hygral and thermal distributions were uniform.
4.2.2. Mesoscopic Numerical Models
In this example, a cylindrical concrete sample with a diameter D = 100 mm and a height H = 200 mm (Figure 4(a)) has been subjected to a fire condition. A three-dimensional model has been carried out in a mesoscale approach where an aggregate random distribution, in agreement with the curve reported in Figure 4(b), has been taken into account.
The concrete mix design in volume is composed of water for a 16%, inerts for a 72%, and cement for a 12%. The inerts have been subdivided as follows: 53% is sand and 47% is gravel.
Only the coarse aggregates have been explicitly represented in the mesoscale model as spheres having different diameters between 3 and 25 mm. Sand and fine aggregates have been taken into account into the cement matrix. The hygral, thermal, and mechanical characteristic for the cement matrix and the coarse aggregates are reported in material parameters. In this example, the ITZ interface between the matrix and the aggregates is neglected.
Two different concrete samples have been evaluated, with and without PP fibres in the matrix. In the first model, the addition of 0.3 kg/m3 of fibres with a length mm has been considered, while in the second model the same concrete without PP fibres has been represented.
In a mesoscale approach, the porosity formulation has to be modified as coarse aggregates have been explicitly represented.
So, porosity in the cement matrix is evaluated as follows:with an initial porosity .
Coarse aggregates are considered with a constant small porosity .
In order to simplify the numerical analyses, symmetry conditions have been assumed (Figure 5), obtaining a numerical model composed of 17600 nodes and 91300 tetrahedral elements with linear shape functions and five degrees of freedom (displacements along the three directions, humidity, and temperature).
A random aggregates distribution in the sample, according to the mix design grading curve of Figure 4(b), has been realized following [23, 24]. Numerically, a continuous mesh discretization has been realized, where no numerical contact condition has been considered at the interface between matrix and aggregates.
Mechanical boundary conditions have been set on the symmetry surfaces and a slow thermal heating curve has been applied on the external surfaces of the cylindrical sample (see Figure 6). No mechanical loads have been considered in this numerical example.
Considering the same thermal conductivity and specific heat for the cement matrix and the aggregates, the thermal distribution inside the sample during the analysis is not affected by the different components, as visible in Figure 7.
As regards the hygral evolution, different results are obtained (Figure 8). The humidity reduction caused by thermal field is influenced by the aggregates distribution, being considered with a low porosity.
Moisture clog is the humidity increment that occurs in the internal part of the sample due to the fact that the temperature pushes water inward during thermal heating, as visible in Figure 9. This effect causes a consequent water pressure increment which, together with the thermal gradient, is one of the spalling causes.
Considering the presence of PP fibres into the concrete mixing there will be a porosity increase, so reducing the water pressure and consequently the spalling risk.
Around the grains of the concrete sample, a lower humidity reduction can be noticed (Figure 10(a)), as well as higher water pressure results if compared to other zones where the distance between the aggregates is bigger (Figure 10(b)).
The small humidity variations around the aggregates can be associated to the initial conditions assumed for the sample. The analysis starts considering that the hygrothermal conditions of the aggregates and the matrix are in equilibrium with the environment (50% of humidity and 20°C of temperature have been assumed). The hydraulic diffusivity in the inert has been assumed very low, causing a slow water release during the analysis that will maintain humidity of the cement matrix, so increasing also the internal water pressure.
The water pressure has been evaluated in agreement with  where is the water vapor saturation pressure, obtained following the Clausius-Clapeyron relation as a function of the molar mass ; the evaporation enthalpy ; and the gas constant :where is the room temperature and is the initial water pressure (see ).
In Figure 11, the water pressure variation inside the two samples (with and without PP fibres) has been reported, where an internal pressure reduction is visible if polypropylene fibres have been added in the concrete mixing.
In the sequence shown in Figure 11, the water pressure is approximately the same in the sample with and without fibres before PP melting point (Figures 11(a), 11(d), 11(b), and 11(e)). This result can be explained considering that the variation of the initial porosity (porosity at room temperature) in the matrix having PP fibres, due to the increment of micro air bubbles trapped in the paste, in the ITZ around fibres, is not very high and therefore the hygrothermic diffusivity does not change too much. Exceeding the melting temperature for PP fibres (that occurs between 165°C and 200°C), the diffusivity in the paste increases, so reducing the water pressure. Diffusivity increases again when temperature gets over 325°C, that is, the polypropylene vaporization point.
The water pressure reduction in the sample with PP fibres can be seen also in Figure 12, where at different times a decrease of is visible in the sample containing polypropylene; this is related to the porosity increment caused by the presence of fibres.
The analyses presented in this subsection, differently from the macroscopic ones reported in Section 4.2.1, did allow seeing humidity and temperature concentrations due to the barrier effect related to the explicit representation of aggregates; that is, accurate evaluations of the local effects in terms of humidity, temperature, and pore pressure.
In this work, the complex mechanism of polypropylene fibres contribution on concrete behaviour at elevated temperatures has been investigated; an additive decomposition of porosity law has been carried out to evaluate the PP effects that occur during thermal evolution.
Three-dimensional hygrothermomechanical finite element models have been performed and validated considering experimental tests by literature to evaluate the diffusivity variations in concrete, with and without fibres addition. As could be noted, the porosity variation law is able to evaluate the reduction of internal pressure in the concrete material when PP fibres are added to the mix design.
Macro- and mesoscale numerical analyses have been carried out; the latter has provided a more realistic description of concrete than the macroscale analysis, as the mesoscale analysis is influenced by the geometry and the properties of multiple constituents of concrete.
The authors declare that they have no competing interests.
G. A. Khoury and B. Willoughby, “Polypropylene fibres in heated concrete. Part 1: molecular structure and materials behaviour,” Magazine of Concrete Research, vol. 60, no. 2, pp. 125–136, 2008.View at: Publisher Site | Google Scholar
C.-G. Han, Y.-S. Hwang, S.-H. Yang, and N. Gowripalan, “Performance of spalling resistance of high performance concrete with polypropylene fiber contents and lateral confinement,” Cement and Concrete Research, vol. 35, no. 9, pp. 1747–1753, 2005.View at: Publisher Site | Google Scholar
G. Xotta, V. A. Salomoni, and C. E. Majorana, “Thermo-hygro-mechanical meso-scale analysis of concrete as a viscoelastic-damaged material,” Engineering Computations, vol. 30, no. 5, Article ID 17090567, pp. 728–750, 2013.View at: Publisher Site | Google Scholar
C. E. Majorana, V. Salomoni, and B. A. Schrefler, “Hygrothermal and mechanical model of concrete at high temperature,” Materials and Structures, vol. 31, no. 210, pp. 378–386, 1998.View at: Google Scholar
G. Mazzucco, C. E. Majorana, and V. A. Salomoni, “Numerical simulation of polypropylene fibres in concrete materials under fire conditions,” Computers & Structures, vol. 154, pp. 17–28, 2015.View at: Publisher Site | Google Scholar
R. Barretta, R. Luciano, and J. R. Willis, “On torsion of random composite beams,” Composite Structures, vol. 132, pp. 915–922, 2015.View at: Publisher Site | Google Scholar
R. Barretta, R. Luciano, and F. M. de Sciarra, “A fully gradient model for Euler-Bernoulli nanobeams,” Mathematical Problems in Engineering, vol. 2015, Article ID 495095, 8 pages, 2015.View at: Publisher Site | Google Scholar | MathSciNet
R. Barretta, L. Feo, R. Luciano, and F. Marotti de Sciarra, “A gradient Eringen model for functionally graded nanorods,” Composite Structures, vol. 131, pp. 1124–1131, 2015.View at: Publisher Site | Google Scholar
F. Greco and R. Luciano, “A theoretical and numerical stability analysis for composite micro-structures by using homogenization theory,” Composites Part B: Engineering, vol. 42, no. 3, pp. 382–401, 2011.View at: Publisher Site | Google Scholar
Z. P. Bažant, M. R. Tabbara, M. T. Kazemi, and G. Pyaudier-Cabot, “Random particle model for fracture of aggregate or fiber composites,” Journal of Engineering Mechanics, vol. 116, no. 8, pp. 1686–1705, 1990.View at: Publisher Site | Google Scholar
S. Eckardt, S. Hafner, and C. Könke, “Simulation of the fracture behaviour of concrete using continuum damage models at the mesoscale,” in Proceedings of the European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS '04), Jyväskylä, Finland, July 2004.View at: Google Scholar
D. Krajcinovic, “Distributed damage theory of beams in pure bending,” Journal of Applied Mechanics, vol. 46, no. 3, pp. 592–596, 1979.View at: Publisher Site | Google Scholar
L. Lanzoni and A. M. Tarantino, “Damaged hyperelastic membranes,” International Journal of Non-Linear Mechanics, vol. 60, pp. 9–22, 2014.View at: Publisher Site | Google Scholar
J. Mazars and G. Pyaudier-Cabot, “Continuum damage theory—application to concrete,” Journal of Engineering Mechanics, vol. 115, no. 2, pp. 345–365, 1989.View at: Publisher Site | Google Scholar
G. Xotta, G. Mazzucco, V. A. Salomoni, C. E. Majorana, and K. J. Willam, “Composite behavior of concrete materials under high temperatures,” International Journal of Solids and Structures, vol. 64, pp. 86–99, 2015.View at: Publisher Site | Google Scholar
Z. Sun, E. J. Garboczi, and S. P. Shah, “Modeling the elastic properties of concrete composites: experiment, differential effective medium theory, and numerical simulation,” Cement and Concrete Composites, vol. 29, no. 1, pp. 22–38, 2007.View at: Publisher Site | Google Scholar
P. Baggio, C. E. Majorana, and B. A. Schrefler, “Thermo-hygro-mechanical analysis of concrete,” International Journal for Numerical Methods in Fluids, vol. 20, no. 6, pp. 573–595, 1995.View at: Publisher Site | Google Scholar
D. Gawin, C. E. Majorana, and B. A. Schrefler, “Numerical analysis of hygro-thermal behaviour and damage of concrete at high temperature,” Mechanics of Cohesive-Frictional Materials, vol. 4, no. 1, pp. 37–74, 1999.View at: Publisher Site | Google Scholar
Z. P. Bažant and L. N. Najjar, “Nonlinear water diffusion in nonsaturated concrete,” Matériaux et Construction, vol. 5, no. 1, pp. 3–20, 1972.View at: Publisher Site | Google Scholar
C. Rossino, F. L. Monte, S. Cangiano, R. Felicetti, and P. G. Gambarova, “Concrete spalling sensitivity versus microstructure: preliminary results on the effect of polypropylene fibres,” in Proceedings of the MATEC Web of Conferences, vol. 6, EDP Sciences, Paris, France, September 2013.View at: Google Scholar
K. Willam, I. Rhee, and Y. Xi, “Thermal degradation of heterogeneous concrete materials,” Journal of Materials in Civil Engineering, vol. 17, no. 3, pp. 276–285, 2005.View at: Publisher Site | Google Scholar
L. T. Phan, “Pore pressure and explosive spalling in concrete,” Materials and Structures, vol. 41, no. 10, pp. 1623–1632, 2008.View at: Publisher Site | Google Scholar
G. Mazzucco, G. Xotta, B. Pomaro, C. E. Majorana, F. Faleschini, and C. Pellegrino, “Mesoscale modelling of concrete with recycled aggregates,” in Proceedings of the 10th International Conference on Mechanics and Physics of Creep, Shrinkage, and Durability of Concrete and Concrete Structures (CONCREEP '15), pp. 853–863, Vienna, Austria, September 2015.View at: Publisher Site | Google Scholar
S. Torquato, Random Heterogeneous Materials: Microstructure and Macroscopic Properties, vol. 16, Springer Science & Business Media, Berlin, Germany, 2013.
R. W. Lewis and B. A. Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley & Sons, 2nd edition, 1998.