S04 Geomechanics and granular materials [download abstracts]

List of abstracts

ID 11: Viscoelastic damage model for rate-dependent shear behavior of shales (keynote) – M. Gutierrez

ID 52: Numerical analyses of hydraulic fracture growth in a 3D rock specimen due to high pressure water R. Abdi, M. Krzaczek, J. Tejchman

ID 55: An effect of content and particle stiffness on the mechanics of particle mixtures: experimental and numerical study (poster) – J. Wiącek, M. Molenda, J. Horabik, M. Bańda, P. Parafinuik

ID 80: Computational method of predicting the crack propagation in confrontation with the laboratory testsJ. Gontarz, J. Podgórski

ID 111: Numerical modelling of concrete-to-concrete interfaces for slim-floor composite slabs analysisS. Dudziak, P. M. Lewiński

ID 120: Essence and usefulness of a new thermodynamic method for determining porosity of solid materials (poster) – Z. Żmudka, S. Postrzednik

ID 139: A quadratic constitutive law for viscous creep of isotropic ice inferred from experimentsR. Staroszczyk, L. W. Morland

ID 180: On the anisotropy and hydromechanical properties of intact/jointed rocks (keynote) – S. Pietruszczak

ID 190: Influence of undercutting anchor geometry on the extent of the rock failure zone during pull-out (poster) – A. Wójcik, J. Jonak, J. Podgórski

ID 226: On a control of contact interaction in the analysis of friction and wearZ. Mróz, I. Paczelt

ID 275: Numerical modeling and improvement of underground blasting efficiencyM. Kucewicz, P. Baranowski, Ł. Mazurkiewicz, J. Małachowski


ID 11

Viscoelastic damage model for rate-dependent shear behavior of shales

Marte Gutierrez1

1 Civil and Environmental Engineering, Colorado School of Mines, United States


Shales are the most ubiquitous material on the earth's surface. Improved understanding of the mechanical behavior of shales is one of the important challenges in the geomechanics community, as shales are often encountered in many energy and environmental engineering applications. The strain rate-dependent mechanical behavior of shale was extensively characterized using triaxial compression tests carried out at different loading rates. A constitutive model for shale under high confining pressure and different strain rates was proposed based on the experimental results. The model is based on a combination of viscoelasticity and damage mechanics, such as the model developed by Katsuki and Gutierrez for asphalt concrete. The model is formulated to predict the quasi-brittle behavior of shales in the pre-peak stage and the strain softening post-peak behavior. The post-peak response is believed to be caused by the growth of microcracks in shales de-bonding and de-cohesion mechanisms responsible for the microcrack evolution. The combined laboratory testing and modeling showed that higher axial loading strain rates lead to higher elastic modulus and higher peak shear strength, with both following exponential function relationships with axial strain rate. Failure results in a single linear shear fracture when the axial strain rate is less than 10-5/s. In contrast, when the axial strain exceeds 10-5/s, more complex development of multiple micro-cracks occurs, forming a crisscrossing fracture network. Failure in shale can be described by a damage variable D which is slightly larger than zero in the compaction stage and is almost zero in the elastic phase, gradually decreases when the load reaches the yield stress, then reduces to a minimum after the yield stress. Taking the yield stress as the cut-off point for damage and considering the evolution of the shale damage variable, a new measure of micro-mechanical strength F is proposed. Based on Lemaitre's equivalent strain assumption and the new variable F, a statistical strain-rate dependent damage constitutive model for shale is established with physically meaningful model parameters. The model predictions are compared with the results of an extensive series of triaxial tests on shale shear behavior under different shearing rates. The excellent agreement between the model predictions and experimental data demonstrates the validity and reliability of the proposed model formulation.

ID 52

Numerical analyses of hydraulic fracture growth in a 3d rock specimen due to high pressure water

Rezvan Abdi1, Marek Krzaczek1, Jacek Tejchman1

1 Faculty of Civil and Environmental Engineering, Gdansk University of Technology, Poland


Hydraulic fracturing is a stimulation technique to extract petroleum, gas, or heat from rock reservoirs that are fractured by a pressurized fluid of about 70 MPa at the target deep depth (of about 3-5 km). The production of petroleum, gas and heat greatly depends on the effectiveness of hydraulic fracturing stimulations. Thus, numerical investigations are crucial for optimizing a hydraulic fracturing process (e.g., pumping rate, fluid viscosity et ca.) [1] to improve their extraction from deep shale layers. The rock masses are characterized by very low porosities which leads to very high packing of particles in numerical simulations. The complex packed structures cause difficulty in mesh generation in common numerical methods like the finite volume method (FVM).

In this study, the FVM was used [2] to perform 3D fluid flow simulations during the propagation of fracture inside a rock specimen that was modelled by densely packed overlapping spheres which were randomly distributed in one layer. This geometry is suitable for the validation and calibration of other numerical approaches like coupled DEM-CFD models due to the lack of valid experimental data. The 3D FVM model was validated with experimental and numerical data from the literature. The compressible single-phase flow of water was considered for simulations under isothermal conditions. The Baseline (BSL) Reynolds Stress model was chosen for turbulent flow analyses. The fluid moved inside the porous domain caused by a pressure drop between the inlet and outlet. Flow characteristics at the pressure drop of 30 MPa inside the porous domain were visualized by streamlines. In addition, the distribution of pressure, velocity, vorticity, and density were computed. Analyzing the results reveals that flow dynamics considerably altered from one pore to another in porous specimens. In addition, the permeability, tortuosity and inertial term coefficient were estimated and compared with well-known relationships.

[1] Bazant ZP, Salviato M, V. Chau VT, Viswanathan H, Zubelewicz A. Why Fracking Works? Journal of Applied Mechanics, 2014: 81.
[2] Abdi R., Krzaczek, M., Tejchman, J. Comparative study of high-pressure fluid flow in densely packed granules using a 3D CFD model in a continuous medium and a simplified 2D DEM-CFD approach. Granular Matter, 10.1007/s10035-021-01179-2, 2021.

ID 55

An effect of content and particle stiffness on the mechanics of particle mixtures: experimental and numerical study

Joanna Wiącek1, Marek Molenda1, Józef Horabik1, Maciej Bańda1, Piotr Parafinuik1

1 Institute of Agrophysics of Polish Academy of Sciences, Poland


The mixtures of granular materials with different elastic properties and content of components have been widely applied over the last few years. The possibility of application of such binary granular systems in many brunches of industry provides necessity to clarify their mechanical behavior under different load conditions.
The objective of this study was to investigate experimentally and numerically the effect of particle stiffness and content of components on the mechanical properties of binary granular mixtures under uniaxial compression.
The series of numerical tests have been conducted with use of the discrete element method (DEM) [1] for two types of spherical particles of varying sizes and material parameters. The initial DEM input parameters corresponded to the material parameters of glass and polylactic acid PLA. The diameters of spheres were 30 mm and 12 mm, respectively. The Young modulus of larger particles was gradually reduced to decrease the ratio between the Young moduli of particles from 8 to 1. Granular packings with various contents of components were generated into a cylindrical chamber 0.21 m in diameter and 0.16 m high. Number of particles varied from 1795 to 30.600. Samples of the same initial volume were subjected to uniaxial compression test. Granular material was compressed through top platen that moved vertically with a constant velocity of 0.005 m/s. A maximum vertical pressure on the uppermost particles was 100 kPa. The volume fraction of larger particles in the sample increased from 0% to 100% providing granular assemblies with different packing characteristics and mechanical properties.
A study of the mechanics of binary mixtures included distribution of contact forces, stress transmission through the contacts, and stiffness of the sample.
Discrete element method simulations were complemented by experimental tests. These were conducted for glass beads and airsoft pellets made of PLA. Experiments were carried out with use of apparatus made of galvanized steel, equipped with the cylindrical chamber with dimensions the same as in numerical tests. Samples of the same size, with different contents of components were prepared and subjected to uniaxial compression.

[1] P.A. Cundall, O.D. Strack, Géotechnique 29, 47-65, 1979.

ID 80

Computational method of predicting the crack propagation in confrontation with the laboratory tests

Jakub Gontarz1, Jerzy Podgórski1

1 Department of Structural Mechanics, Lublin University of Technology, Poland


The presented paper is a continuation of our work, which focuses on the search for effective and accurate methods of simulating the crack propagation processes in quasi-brittle materials subjected to complex stress states. The result of these works is a method of predicting the direction of propagation, which was implemented as a subroutine attached to the popular FEM system Abaqus, in which the X-FEM technique is used to simulate the crack propagation process. In our previous publications, we showed that the proposed algorithm gives results closer to reality.
This paper concerns the verification of the above method in relation to the laboratory testing of the pull-out of an anchor embedded in a cylindrical concrete sample. The sample is a cylinder with two M10 bolts embedded on both sides in one axis. Two variants of the test were performed - only the screw thread is countersunk and the inverted screw with a countersunk head. The test consists in applying a tensile force to two bolts to break the specimen. The tear-out force and the course of the fracture line were observed.
The above study was modeled in the Abaqus system as an axially symmetrical task. The material adopted in the simulation is concrete C25 with the following parameters: Young's modulus E = 30 GPa, Poisson's ratio v = 0.2, tensile strength ft = 2.6 MPa, Critical energy of failure in the mode I KI = 0.1 N/mm. The simulation results were compared with the results of laboratory tests.
Simulations were performed using the default method built into the Abaqus system, i.e. a simple implementation of the maximum principal stress criterion. It assumes that the direction of the crack is the same as the direction of the line perpendicular to the direction of the maximum principal stress, which for complex states gives results significantly different from reality. It is because the algorithm used to find the principal directions of the stress tensor near the crack tip is too simple. Therefore, a method was developed that gave much more precise values of the stress field by approximating the FEM results at many points located near the crack tip, which allowed for the correct implementation of several criteria.
This paper compares the results obtained using the above criteria and presents the advantages of the developed method. The possibilities of using this method to determine the strength parameters of the tested materials were also proposed.

ID 111

Numerical modelling of concrete-to-concrete interfaces for slim-floor composite slabs analysis

Sławomir Dudziak1, Paweł M. Lewiński2

1 Faculty of Civil Engineering, Warsaw University of Technology, Poland
2 Building Research Institute (ITB), Poland


Concrete-to-concrete interfaces can be found in many modern composite structures, e.g. slim-floor slabs on hybrid beams made of high performance concrete and high strength steel [1] or composite T-beams. They influence load capacity and deformability of such structures a lot, therefore the issue of determining their mechanical characteristics is a current and important research topic.

In the present work, a strategy for calibration of concrete-to-concrete interfaces is presented. It is based on cohesive finite elements and traction-separation constitutive law formulated in the framework of damaged elasticity, which are available by default in the Abaqus code [3,4]. The standard strength envelope from the Abaqus code, which is not suitable for concrete-to-concrete interfaces, is modified with USDFLD user's subroutine. The proposed model takes into account different post-cracking behaviours in tensions and shear traction states as well. However, some limitations of the proposed approach are apparent, e.g. inability to reproduce of shear-friction phenomenon in fully correct manner.

Then, the proposed numerical model is validated using results of popular laboratory tests performed on small specimens, i.e. a splitting, a slant-shear test and a three-point bending of a beam with a notch. For concrete regions, brick elements with the selective integration (C3D8) and the concrete damage plasticity (CDP) constitutive model are used. The agreement between results of the laboratory tests and the performed simulations is satisfactory.

Finally, the applicability of the proposed approach to the simulation of full-scale composite systems is demonstrated by the non-linear analysis of a slim-floor slab on hybrid beams. The results of the analysis are confronted to the experimental outcomes previously presented in the paper [1]. The numerical model is able to capture behaviour of this complex structural system, i.e. load-deflection curve and expected failure modes are correctly reproduced. Proposals for the model refinement are also suggested.

[1] P. M. Lewiński, J. Derysz, S. Dudziak, and P. Więch, Newly designed structural solutions for the 'slim floor' composite system, Proc. fib Symp. 2019 Concr. - Innov. Mater. Des. Struct., no. May, pp. 873-880, 2019.
[2] N. S. Ottosen and M. Ristinmaa, Thermodynamically based fictitious crack/interface model for general normal and shear loading, Int. J. Solids Struct., vol. 50, no. 22-23, pp. 3555-3561, 2013.
[3] S. Dudziak, W. Jackiewicz-Rek, and Z. Kozyra, On the Calibration of a Numerical Model for Concrete-to-Concrete Interface, Materials (Basel)., vol. 14, p. 7204, 2021.
[4] Abaqus Theory Guide, Dassault Systemes 2013.

ID 120

Essence and usefulness of a new thermodynamic method for determining porosity of solid materials

Zbigniew Żmudka1, Stefan Postrzednik1

1 Silesian University of Technology, Poland


Among the various porous materials, two categories should be indicated as the most typical: granular bed and porous structures with a rigid framework (e.g. catalyst substrate). The basic parameters that determine the structure of any solid porous system include the following:
- specific surface area (reaction area) of solid phase,
- porosity of solid material and its structure,
in the case of a granular (polydisperse) bed, the additional information shall be:
- size analysis (mass frequency function of the grains).
Specific surface area determines first of all intensity of surface phenomena e.g. chemical reactions in the bed. Bed porosity influences mainly transport conditions of continuous (gaseous) phase through solid material, ensuring among other things constant inflow of substrates to superficial reaction zone as well as outflow of products.
The research carried out showed that the porosity of solid materials should be determined accurately as far as possible using a suitable measurement method, i.e. a measurement method that is independent of the flow characteristics of the granular bed or solid structure. For this purpose, the own, original method of determination of any structure porosity was developed, which has been called thermodynamic method. The method takes advantage of the thermodynamic properties (mainly compressibility) of the gas filling the intergranular space of the solid structure. An analytical algorithm was developed to calculate the porosity. In fact, with the developed method it is possible to determine so-called total porosity of a granular deposit, which in the case of homogeneous grains (without internal gaps) is equal to the intergranular (external) porosity. Non-homogeneous grains additionally show internal porosity.
In order to verify the practical usefulness of the developed measurement method, a number of validation tests were carried out. Among other things, the porosity of the bed consisting of grains of known geometry was investigated. The total porosity of the ceramic catalyst substrate of an IC engine was also tested. The porosities of these structures were determined in two ways:
- on the basis of the known geometry of the tested structures,
- using the presented thermodynamic method and calculation.
A comparison of the porosity calculated from the geometry and the porosity determined by the presented method showed that the relative error is between 1 and 3%. The article presents the essence of the method, sample measurement results and an assessment of its accuracy.

ID 139

A quadratic constitutive law for viscous creep of isotropic ice inferred from experiments

Ryszard Staroszczyk1, Leslie W. Morland2

1 Institute of Hydro-Engineering of the Polish Academy of Sciences, Poland
2 School of Mathematics, University of East Anglia, United Kingdom


A viscous fluid flow law for isotropic ice is formulated, in which the response of ice to applied stress is described by a general frame-indifferent Rivlin-Ericksen quadratic relation in the strain-rate tensor, with two response coefficient functions assumed to depend on the second principal invariant of the strain-rate tensor and on temperature. The creeping ice is treated as an incompressible and thermo-rheologically simple material, in which the strain-rate at a given stress can be normalized by a temperature-dependent rate factor.
Three different forms of the general quadratic law for the viscous creep of isotropic ice have been examined, for which analytic relations for uni-axial compression and simple shear configurations have been derived, and the criteria necessary for the validity of the predicted ice responses have been determined. Lacking simple shear experimental data for ice, an alternative analysis is made of Steinemann's (1958) torsion experiment configuration, in which measurements have been performed on ice samples in the form of a hollow cylinder that is subjected to a torque applied on its top horizontal surface. The results of correlations of the proposed constitutive law with the combined uni-axial compression and torsion experimental data sets show that the quadratic term in the flow law is significant over the stress and strain-rate ranges investigated in the laboratory measurements, with the maximum ratio of the quadratic term contribution to that of the linear term equal to about 23%. This clearly implies that the conventional coaxial (linear) constitutive relation with the response coefficient depending only on one strain-rate invariant is not valid in those stress and strain-rate ranges. Extrapolation of the laboratory Steinemann’s data to the low stress levels arising in large polar ice-sheet flows indicates that the quadratic term contribution is also not negligible, since it reaches about 12% of the linear term contribution to the total stress in ice.

ID 180

On the anisotropy and hydromechanical properties of intact/jointed rocks

S. Pietruszczak1

1 McMaster University, Canada


This paper deals with description of hydromechanical response of intact rocks as well as discontinuous rock mass. In the former case, the focus is on the notion of anisotropy, which is addressed in the context of argillaceous Cobourg limestone. For a rock mass containing discontinuities, the average hydro-mechanical properties are established in a domain that comprises the intact material and the pre-existing and/or newly developing fracture zones.

The anisotropy of Cobourg limestone stems from heterogeneity of its fabric. In this work, an approach that incorporates a fabric tensor is developed for the description of strength and deformation characteristics of this rock. The quantification of fabric is carried out by employing the basic principles of stereology and the specific fabric descriptor used is the so-called 'mean intercept length'. The results demonstrate that the material may be perceived as transversely isotropic and the mechanical response is affected by the volume fraction of argillaceous partings. The latter are spatially variable, so that the notion of Representative Elementary Volume (RVE) is not, in general, applicable. A numerical study is included involving simulation of triaxial tests conducted on differently oriented samples of Cobourg limestone.

For a rock mass containing discontinuities, the formulation employs the averaging of the field operators within the referential volume adjacent to macrocrack. This leads to an enriched form of Darcy's law, which incorporates the notion of equivalent hydraulic conductivity. The latter is defined as a symmetric second-order tensor whose components are function of hydraulic properties of constituents (viz. intact material and fractured region) as well as the internal length parameter. Such an approach does not require any additional degrees of freedom to account for the presence of discontinuities, which significantly improves the computational efficiency.

The mechanical analysis is based on an enhanced embedded discontinuity approach, which is conceptually similar to that employed for specification of equivalent hydraulic conductivity. It incorporates the same internal length parameter related to geometry of fractures and enables a discrete tracing of the propagation of new cracks.

The formulation is illustrated by a numerical study involving 3D simulation of an axial splitting test carried out on a saturated sample under displacement and fluid pressure-controlled conditions. The finite element analysis incorporates the PPP stabilization technique and a fully implicit time integration scheme.

ID 190

Influence of undercutting anchor geometry on the extent of the rock failure zone during pull-out

Andrzej Wójcik1, Józef Jonak1, Jerzy Podgórski2

1 Department of Machine Design and Mechatronics, Lublin University of Technology, Poland
2 Structural Mechanics Department, Lublin University of Technology, Poland


The paper presents results of numerical analysis of influence of geometric parameters of undercutting anchor head on forming of failure zone range (failure cone) in the initial stage of medium failure development. The issue is important for optimization of rock removal technology with the use of shear bolts. In order to assess the effectiveness of the process, one of the indicators is the volume of the potentially spalling rock under the action of the anchor. Determination of factors influencing the range of rock dislodgement and their most favorable combination for the process is a subject of a broader study, a fragment of which, concerning the influence of geometrical parameters of the anchor head on the propagation and range of the rock failure zone, is the subject of this research.

ID 226

On a control of contact interaction in the analysis of friction and wear

Zenon Mróz1, Istvan Paczelt2

1 Institute of Fundamental Technological Research, Poland
2 University of Miskolc, Hungary


The present work is related to control of local boundary displacement of a loaded elastic structure by introducing additional action of one or two punches. A simplest control is aimed at assuring a fixed normal displacement value u = u0 at a selected point of a loaded boundary. A more advanced control is applied when a body enters in contact and the required force-penetration interaction Pn= Pn(un) is executed by applying punch action. Such problem formulation generates a new type of optimization problems, discussed already in [1] and [2]. The analysis results can be applied to design of robot grippers or tools of mechanical processes (such as riveters) or control and measurement devices. For some interactions the control of normal displacement and its tangential component related to the displacement gradient (slope) is executed by the synchronized action of two punches. In such class of problems the proper design of acting punches, selection of boundary constraints in order to preserve elastic stress states is important. Finally, for periodic action of punches, the attached indenter performs a combined sliding and penetration motion at the contact interface of a substrate element. In practical applications the flexural deformations of beams and plates can be transformed to 2D or 3D indenter tool motion.

[1] Páczelt, I. and Mróz Z. Optimized punch contact action related to control of local structure displacement. Structural and multidisciplinary optimization, 1, 1921-1936, 2019.
[2] Páczelt, I. and Mróz Z. A New Class of Optimization Problems Related to Structural Control by Contact Interaction, In: Krivtsov A.M.; Indeitsev D.A. (eds.), Advanced Problems in Mechanics, Springer International Publishing, 2020.

ID 275

Numerical modeling and improvement of underground blasting efficiency

Michał Kucewicz1, Paweł Baranowski1, Łukasz Mazurkiewicz1, Jerzy Małachowski1

1 Faculty of Mechanical Engineering, Military University of Technology, Poland


Underground blasting is a commonly used method for excavation of orebodies contained in the hard rocks as dolomite. This method induces a high costs of mining, which is demanded to be reduced as much as possible. The rapid development of numerical simulations affords for reliable reproduction of rock behavior subjected to dynamic loadings that are induced in blasting. Simulations could be adopted for optimization of borehole spacing, dimensions, filling etc., to maximize the excavated material to costs ratio. This task, however, requires a careful selection of constitutive material model for rock. Several demands such a faithful reproduction of hardening, softening, strain rate and pressure dependency as well as cracking and fragmentation have to ensured by the model.

In this paper three constitutive models: Karagozian and Case Concrete (KCC) [1], Johnson-Holmquist 2 (JH-2) [2] and Johnson-Holmquist-Concrete (JHC) [3] will be compared in terms of modeling a dynamic experimental test. Three experiments were selected to verify the performance of mentioned models: uniaxial dynamic compression, ball bearing drop test and small scale blasting. The selection resulted from covering by the tests the most critical aspects of the material's response in compression and tension. Lately, the selected model was applied for numerical simulations of cut hole and underground blasting, and various loading scenarios and blasting sequences were considered. The JH-2 model, that was applied for these simulations, allowed to verify how to maximize the crushing of the material in blasting. The numerical aspects, such as reproduction of HE, mapping of quasi-plane strain model based on 3D geometries and its effects for simulations of underground blasting were investigated.

[1] M. Kucewicz, P. Baranowski, J. Małachowski, Determination and validation of Karagozian-Case Concrete constitutive model parameters for numerical modeling of dolomite rock, International Journal of Rock Mechanics and Mining Sciences. 129 (2020).
[2] P. Baranowski, M. Kucewicz, R. Gieleta, M. Stankiewicz, M. Konarzewski, P. Bogusz, M. Pytlik, J. Małachowski, Fracture and fragmentation of dolomite rock using the JH-2 constitutive model: Parameter determination, experiments and simulations, International Journal of Impact Engineering. 140 (2020).
[3] M. Kucewicz, P. Baranowski, J. Małachowski, Dolomite fracture modeling using the Johnson-Holmquist concrete material model: Parameter determination and validation, Journal of Rock Mechanics and Geotechnical Engineering. 13 (2021) 335-350.