Modeling of rigid inclusion ground improvements in large-scale geotechnical simulations

An increasing number of studies highlight the significance of conducting large-scale geotechnical simulations (e.g., Lu et al., 2011, Łydżba et al., 2021). Such simulations are essential for the accurate design of large-scale structures, including earth dams, heaps, and tailings storage facilities. These computations frequently encompass areas spanning several square kilometers. To ensure the results are accurate and reliable, an appropriate level of discretization is required. In these types of problems, mesh refinement can lead to a substantial increase in the number of nodes. Such increase significantly extends the computation time, rendering numerical modeling of this scale computationally prohibited for real-world engineering applications. However, certain simplifications assumed in such models can lead to a satisfactory compromise between the number of nodes and the quality of results.

One of the simplest approaches to significantly reduce the number of nodes is to simplify a three-dimensional (3D) task to a two-dimensional (2D) one. Instead of solving a full 3D problem, it is sometimes possible to identify a characteristic cross section and perform calculations under the assumption of plane strain conditions or to find an axis of symmetry and treat the problem as axisymmetric. However, this simplification is not always feasible. If the subsoil conditions are complex, that is, the layers are not parallel, the top and bottom of subsequent layers are highly variable, or there are discontinuous layers, converting a 3D problem to 2D can lead to significant errors. A similar situation occurs when the designed structure varies in all three directions, making it difficult to identify a characteristic cross section. In such cases, it is necessary to use 3D modeling, which, as mentioned earlier, is very time-consuming.

To handle large 3D problems, researchers and engineers continuously seek effective methods and computational approaches. An example of such an approach is the computational method for designing dams of dry flood reservoirs proposed in works by Łydżba et al. (2021) and Sobótka et al. (2022). In these works, the authors introduced a comprehensive computational approach to design a dam according to the European Geotechnical Standard (EN 1997-1:2004 Eurocode 7). It is important to note that this approach is dedicated to linear structures such as flood embankments and flood reservoirs, where the geometry of the construction is relatively constant in one direction; however, significant variability in subsoil can still occur in that direction. The major streamlining of the proposed methodology was to avoid using full poromechanical coupling in 3D calculations (Truty et al., 2020). For the evaluation of pore pressure values, the authors propose 3D analyses, which are necessary to account for the spatial variability of the subsoil. Deformation and stability analyses, however, are performed on 2D models, taking into account pore pressure calculated in 3D model. The factor of safety (FOS) for 2D calculations is usually smaller than for 3D (Fredlund et al., 2010); therefore, this approach should be considered safe.

In large-scale geotechnical projects, as in most other structures, the structure must meet the requirements of both ultimate limit state (ULS) and serviceability limit state (SLS) (Bogusz & Godlewski, 2019; Brinkgreve & Post, 2015). However, this is not always feasible with direct foundations due to the insufficient mechanical and strength properties of the underlying soils. In such situations, subsoil reinforcement becomes necessary. When weak subsoil exists at relatively shallow depths, a highly effective method of reinforcement is using rigid inclusions. This technique is particularly effective for reducing settlements (Jiang et al., 2014) and thus fulfilling SLS requirements. Rigid inclusion ground improvement, sometimes referred to as columns, involves inserting into the soil objects, which are characterized by significantly greater stiffness than the surrounding soil (Polańska & Rainer, 2020). The use of such inclusions allows for the unloading of soil between the columns. Various construction technologies exist for rigid inclusions, including deep soil mixing (Jończyk-Szostek et al., 2023), full displacement columns (Krasiński, 2023), and gravel columns (Boru et al., 2022). The choice of technique depends primarily on soil conditions and specific improvement objectives.

While the inclusions can be sometimes necessary, modeling such objects, whose dimensions are a few hundred times smaller than the overall problem, is problematic for large-scale models. This difficulty arises due to the limited possibility of mesh refinement, as the dimensions of individual mesh elements often exceed the dimensions of the inclusions themselves. This paper aims to present the basic methods for modeling rigid inclusions in large-scale geotechnical simulations. The paper introduces three basic methods of modeling inclusions, that is, modeling the columns by continuum elements, beam elements, or modeling the column and surrounding soils by an equivalent medium with parameters determined by homogenization theory. The modeling principles are demonstrated using a periodic cell containing a single inclusion. The application of these methods to real projects is illustrated through the analysis of the foundation of the conduit pipe for the Szalejów Górny antiflood reservoir in Poland.

The remainder of this work is structured as follows: the next section presents numerical methods for modeling rigid inclusions, concluding with a comparative analysis of the selected methodologies. Following this, an example of a dry reservoir dam, where the foundation must be placed in reinforced soil, is discussed. The paper concludes with a discussion of the obtained results.

A substrate with rigid columns can be modeled using finite elements according to three basic methodologies: –

Modeling both soil and columns as continuum elements (denoted in this paper as Approach I).

To compare these methods, a reference task was modeled according to the different approaches, serving as a benchmark for further analysis. The comparative analysis focused on the displacement and stress distribution maps.

The benchmark consisted of a single periodic cell with dimensions 2.8 m × 2.8 m, containing a concrete column with a diameter of 0.88 m and a length of 10 m. A 1.0-m-thick concrete slab was placed on the improved ground and loaded with a uniformly distributed load of 500 kPa. Materials were modeled using an elastic model: soil with a Young’s modulus E = 100 kPa and Poisson’s ratio ν = 0.2; and the inclusions and slab with parameters typical for concrete, E = 30 GPa and ν = 0.2. Displacements in the horizontal directions were constrained on the vertical faces of the model, while both vertical and horizontal displacements were constrained on the bottom surface. The scheme of the benchmark is shown in Fig. 1. All calculations were performed using ZSoil software (Truty et al., 2020).

Figure 1:

Scheme of the reference problem.

2.1

Inclusion as a continuum element (Approach I)

The first method of modeling the column–soil system involves representing all components as continuum elements (ASIRI National Project, 2017). A 3D model of such a system offers the most precise depiction of the actual geometry and behavior of the soil–column system, which is the primary advantage of this approach. However, this method’s drawbacks include the complexity and labor-intensiveness of the modeling process, as well as the considerable computation time required. The Mohr–Coulomb interface on the surface separating soil from the column can be included in the model. This allows for the description of displacement discontinuities along the column shaft and potential shearing at the soil–column interface (Van Langen & Vermeer, 1991; Jalali et al., 2012). The benchmark model used for the comparative analysis is illustrated in Fig. 2.

Figure 2:

Numerical model for Approach I: (a) model subdivision by type of continuum material; (b) finite element mesh-axonometric view; and (c) finite element mesh-top view.

In Fig. 2a, the model is schematically presented with division into areas with different mechanical strength parameters. The elastic parameters assigned to the group of elements modeling the column and soil are identical to the ones presented in Fig. 1. In the depicted model, only the continuum element can be distinguished. Since the columns typically have a circular cross section, the finite element mesh needs to be significantly refined to accurately correspond to such geometry. Fig. 2b and c illustrates the finite element mesh in axonometric and top views, respectively.

2.2

Inclusion as a beam element (Approach II)

The second method of modeling rigid inclusion involves embedding the beam elements within the volume finite elements modeling the subsoil. The beam element does not need to satisfy the compatibility condition between discretization of the soil and the column. This lack of compatibility is realized through a special interface, which ensures displacement compatibility between a selected beam node and the continuum element in which the node is embedded. It is possible to model the interface with the Mohr–Coulomb criterion, allowing for the simulation of potential shearing mechanisms at the soil–column interface, similar to Approach I. Detailed theoretical aspects of this approach are discussed in the work of Truty & Podleś (2009).

The benchmark model for Approach II, with the column represented as a beam element, is shown in Fig. 3. Notably, in Fig. 3a, the column is not visible in the cross section because it is not modeled using continuum elements. Consequently, since the column is not described volumetrically, there is no need for significant mesh densification. Comparing Fig. 3c with Fig. 2c, it is evident how much simpler the mesh geometry can be. It is worth emphasizing that in this particular case, the column is located at the intersection of the symmetry axes of the problem, and therefore at a position where nodes are present (Fig. 3c). However, this is not a necessary condition in the general case.

Figure 3:

Numerical model for Approach II: (a) model subdivision by the type of continuum material; (b) finite element mesh-axonometric view; and (c) finite element mesh-top view.

In this approach, all the volumetric elements representing the soil and the elastic parameters assigned to them are identical as those assumed for soil in Fig. 1. For the beam, an elastic modulus E corresponding to the column, namely 30 GPa, was assumed, and the cross section was assumed to be circular with a diameter of 0.88 m.

2.3

Inclusion and surrounding soil replaced by a homogeneous equivalent medium (Approach III)

The last method presented to model the column–soil system uses homogenization theory. To introduce the main assumptions of homogenization theory, it is necessary to introduce two different scales of observation, namely the microscale and the macroscale. The microscale refers to a situation when the parameters of a given material vary with location, and thus, microscales can be expressed as heterogeneity scales (Różański et al., 2013). On the macroscale, the values of the parameters do not depend on the location, and therefore, the given material is homogeneous (Łydżba, 2011). It should be emphasized here that the macroscale is the typical scale of engineering applications. Applying homogenization theory allows us to replace the real heterogeneous medium with an equivalent homogeneous medium, which, on the macrolevel, exhibits the same response as the heterogeneous material. In the case of soil with rigid inclusions, two main components can be distinguished on the microscale, that is, the column and the surrounding soil.

There are three main steps in the procedure for determining the parameters of continuum homogenous material equivalent to the column–soil system: 1)

Separate a representative elementary volume (RVE) from the heterogeneous medium

In general, determining RVE for various media, including random ones, is not straightforward (Rózański et al. 2013). However, rigid inclusion ground improvements are typically implemented as regular square or triangular grids. For such periodic composite, RVE typically corresponds to a periodic cell. Therefore, RVE dimension and geometry are identical to the cell modeled in Approach I.

2)

Solve the boundary problem for RVE

A sequence of boundary problems must be solved for a single column–soil periodic cell. This appropriately chosen sequence of boundary problems aims to estimate the relationship on a microscale between different values of load/stress and the corresponding state of deformation.

3)

Determine the equivalent parameter of a homogenous medium providing the same macroscopic response as RVE using the homogenization technique, for example, volume averaging methods

The goal of the homogenization technique is to determine the effective stiffness tensor Cijkhhom C_{ijkh}^{hom } based on the microscopic fields of displacements u_i, stresses σ_ij, and strains ε_kh, which were identified in the previous step. This tensor is used to describe the macroscopic constitutive relations according to the equation: (1) <σij>=Cijkhhom<εij> < {\sigma _{ij}} > = C_{ijkh}^{hom } < {\varepsilon _{ij}} >

where <σ_ij> and <ε_ij> denote the macro stress and strain, which can be determined based on microscale results using the volume averaging methods, according to to Eqs (5) and (6) (2) <σij>=1VRVE∫VRVEσijdV < {\sigma _{ij}} > = {1 \over {{V_{RVE}}}}\int\limits_{{V_{RVE}}} {{\sigma _{ij}}dV} (3) <εij>=1VRVE∫VRVEεijdV < {\varepsilon _{ij}} > = {1 \over {{V_{RVE}}}}\int\limits_{{V_{RVE}}} {{\varepsilon _{ij}}dV}

The procedure to determine the effective parameter is described in detail in Łydżba (2011). Please note that since the current analysis is focused on elastic parameters, modeling of the nonelastic interface is not performed on a microlevel, and thus is not described on the macrolevel. While in a general case, such an interface could be taken into account (given that both elastic and inelastic behaviors are homogenized); it significantly complicates the analysis and is not considered here.

Fig. 4 illustrates the reference task for Approach III, where the soil and inclusion are replaced by an equivalent medium. The inclusion–soil system exhibits obvious anisotropic behavior. In the vertical direction (along the axis of the single column), the stiffness of this system is significantly higher than in the perpendicular directions. The equivalent medium exhibits the same characteristics: its stiffness is considerably greater in vertical direction than in the horizontal directions, where it is very close to that without the column. However, for simplicity (including limitations of numerical programs adopted for large-scale geotechnical simulations), it was assumed that the equivalent medium would be treated as isotropic and the assigned parameters would correspond to those obtained for the vertical direction. Thus, for the equivalent medium in the example shown in Fig. 4, the parameters were estimated according to the above-mentioned procedure as E_eff = 2431 MPa and ν_eff = 0.01.

Figure 4:

Numerical model for Approach III: (a) model subdivision by the type of continuum material; (b) finite element mesh-axonometric view; and (c) finite element mesh-top view.

Similar to Approach I, only continuum elements can be distinguished in this analysis. However, on comparing Figure 4c with Figure 2c, it is evident that the grid used in this approach is significantly simplified (similar to Approach II).

In this section, the result obtained with approaches proposed in the previous section are compared. To allow the comparison of methods from Approaches I and II with Approach III, perfect elastic contact is assumed in Approaches I and II on the inclusion’s side surface, which prevents shearing on the side surface. The results of simulations obtained directly from all considered approaches will be discussed first.

Fig. 5 presents the displacement maps according to different modeling methods. Fig. 5a–c shows the displacement maps on the boundary surface, and Fig. 5d–f corresponds to displacement maps in the plane intersecting the column axis of symmetry. Fig. 5a and d corresponds to Approach I, Fig. 5b and e to Approach II, and Fig. 5c and f to Approach III.

Figure 5:

Comparison of displacement maps according to selected design approaches.

The analysis of Fig. 5 shows that the displacement maps on the surface of the model are consistent across all methods (Fig. 5a–c). In addition, the calculations indicate that the final value of the plate displacement is very similar across all approaches. For an assumed load of 500 kPa, displacements are approximately equal to 5.0 cm at the top of the plate. This demonstrates that, for the estimation of displacement, all the aforementioned approaches can be used interchangeably. The only noticeable differences in the displacement maps appear in the cross section through the vertical column symmetry axis. In Approaches I and II (Fig. 5d and e, respectively), it is evident that the inclusion is stiffer and settles more at the column base level than the surrounding soil at the same depth. In Fig. 5f, representing Approach III with a homogeneous medium, the cross-sectional maps are identical to those on the boundary surface shown in Fig. 5c.

The observable differences between the selected methodologies are more visible in the stress state. The maps of the total vertical normal stress are presented in Fig. 6. The results are presented with the same division into rows and columns as shown in Fig. 5. The color scale of the bar in Fig. 6 was selected based on Approach III with an equivalent medium (Fig. 6c). In this approach, the vertical stresses increase linearly with depth. The minimum stress (in absolute terms) occurs at the top of the slab and is equal to 500 kPa, which equals the value of the surface load acting on the top surface of the slab. With depth, the stresses increase, and their increment is related solely to the self-weight of the medium. Since the color scales for Fig. 6a and b were assumed to be identical to that for Fig. 6c, the red color in these two part figures indicates values of 500 kPa and below, while the dark blue color indicates values exceeding 1000 kPa. As seen in Fig. 6a for Approach I, the stresses in the column are greater than 1000 kPa, whereas in the surrounding soil, they are less than 500 kPa. This confirms that, due to its stiffness, the column bears a significantly higher stress, numerically validating the effectiveness of ground improvements. In Approach II (Fig. 6b), similar to Approach I (Fig. 6a), the soil at the depths of the column (up to 10 m) bears much smaller loads. In both Approaches I and II (Fig. 6a and b), a stress concentration zone can be observed in the soil directly below the base of the inclusion. These concentrations gradually spread to the surrounding soil, so that at a depth of about 12 m, the stress maps according to all approaches are similar (Fig. 6a–c).

Figure 6:

Comparison of total vertical stress maps according to selected design approaches.

From the presented comparison of all approaches, it can be concluded that all methods lead to the same results for surface displacement. However, significant differences can be observed in the stress state, particularly vertical stress. The stress maps effectively illustrate approaches where the column is modeled as a separate element, whether as a beam (Fig. 6b) or as a continuum (Fig. 6a). This is especially important if a structure to be dimensioned is located directly over the columns. In the homogenization-based approach, there are no stress concentrations in the slab (Figure 6c). Therefore, dimensioning the slab under this condition may be unconservative, as it does not account for local stress concentrations. The most distinctive of these concentrations are visible in Fig. 6a (Approach I).

The purpose of this work is to evaluate the approaches for applications to large-scale numerical simulations. In this context, additional factors such as computation time, ease of model modification, and the number of elements (nodes) required for modeling should be analyzed.

To accurately represent the shape of the column in Approach I, the circle describing the cross section of the column had to be divided into 32 equal segments. This ensured that the difference in the areas of the circle and the inscribed regular polygon with 32 edges was less than 1%. For Approaches II and III, such a dense grid was not necessary, so a 4 × 4 square grid in the cross section was used. The effect of finite element mesh density on the results was not analyzed in this paper; however, for these meshes, the difference in the number of elements exceeds a factor of 10 (Approach I used about 12,000 nodes, while Approaches II and III used about 1,200 nodes). Although further detailed study and analysis on the mesh used might slightly reduce this difference, it will likely remain significant.

The number of elements directly impacts the computation time. For the prepared benchmark, Approach I required approximately 210 s, while the other two approaches took only 10 s each. It is important to note that there is no clear difference in computation time between Approach II and Approach III. In addition, the use of computer resources such as RAM and CPU is significantly higher for Approach I.

From a practical standpoint, the number of finite elements, computational time, and use of computer resources almost disqualify Approach I with the column modeled as volume elements in modeling of large-scale geotechnical problems. Furthermore, during the design stage of a structure, it is not always initially apparent that ground improvement is needed. The necessity of the application of rigid inclusions often becomes evident after multiple simulations. In Approach I, this necessitates remaking the model and remeshing it, which leads to significantly more working hours for the designer. In Approach II, the beam elements are embedded in the continuum and do not need to completely conform to the existing mesh. In Approach III, adding rigid inclusions only requires changing the mechanical parameters in selected elements of the model. Moreover, the design process often involves parametric analysis, where the length and diameter of the inclusions are determined based on subsequent simulations with varying range, spacing, and diameter of the improvement. In Approaches II and III, the changes in column geometry can be applied easily and quickly, whereas in Approach I, such changes require remeshing, making this method impractical for iterative and parametric design processes. According to the authors, only simplified Approaches II and III are applicable in large-area geotechnical models.

The effectiveness and operational mechanisms of the aforementioned methods in large-scale simulations were demonstrated using a computational model of the earth dam at the Szalejów Górny dry flood control reservoir, located on the Bystrzyca Dusznicka River in Poland. This reservoir was constructed as part of the Odra-Vistula Flood Management Project, co-financed by the World Bank (The World Bank, 2021). The project, situated in the Klodzka Valley in southern Poland, encompassed the construction of four reservoirs: Krosnowice, Szalejów Górny, Roztoki Bystrzyckie, and Boboszów, with Szalejów Górny being the largest. The key characteristics of the Szalejów Górny reservoir are as follows: a dam height of 19.3 m, a dam length of 735 m, a capacity of 9.9 million m³, and a flooding area of 118.7 ha (Łydżba et al., 2021).

The primary function of a dry anti-flood reservoir is to intercept and temporarily retain the passing flood wave, releasing it downstream when it can be done safely. These reservoirs are typically empty, only filling with water when there is a risk of flooding. They are not designed for long-term water retention. Reservoir dams generally consist of an impermeable core, a body, and an anti-filtration barrier that extends the filtration path. In addition, these dams are equipped with a special concrete conduit pipe through which the river normally flows in the absence of flood risk. During periods of high water inflow, such as floods, when the inflow exceeds the conduit pipe’s capacity, water retention occurs within the reservoir. An example of a cross section through the dam at the location of the conduit pipe is shown in Figure 8.

Figure 8:

Cross section through the dam.

The conduit pipe in this particular dam is composed of 11 reinforced concrete segments, each measuring 10 m by 7.5 m by 8.0 m. Due to their robust and massive dimensions, conduit pipes are generally not sensitive to the overall magnitude of settlements. However, the greatest threat to their structural safety is nonuniform settlements. This is primarily because special strips are placed in the dilation joints between successive segments. These strips can be damaged if the differential settlement exceeds 1.0 cm. Consequently, from the perspective of earth structure design, controlling differential settlements is crucial.

The dam subsoil consists primarily of Quaternary and Tertiary formations. The near-surface layers directly beneath the dam are composed of Quaternary silty clay and medium sand with varying consistency and density indices. The physical and mechanical parameters of these materials were determined based on laboratory tests. Beneath the Quaternary formations, there are clays and marls with highly variable properties. The subsoil contains both continuous and fractured layers, leading to a wide dispersion of deformation parameters.

All calculations were performed using ZSoil software (Truty et al., 2018). Simulations assumed full poromechanical coupling and employed the Van Genuchten model (Van Genuchten, 1980) to describe the water flow process. The deformation process of the subsoil was modeled using the elastic–plastic Mohr–Coulomb model. For the reinforced concrete elements, such as the conduit pipe, an elastic model was assumed. A numerical model of the entire dam, including a zoomed-in view of the conduit pipe, is shown in Fig. 9.

Figure 9:

Numerical model of the area under consideration.

The layer arrangement in the subsoil (see Fig. 9) was reconstructed based on available geological documentation, including geological boreholes and the built-in tools in the ZSoil software that utilize kriging techniques. The parameters adopted for the calculations are presented in Tab. 1 for soils and in Tab. 2 for rocks.

Calculations accounted for time-varying loads as well as the sequential nature of the construction process. This method reflects how the structure is built over time and takes into account the execution of construction activities. Fig. 10 illustrates the example stages of the dam construction.

Figure 10:

Visualization of construction stages of the dam.

The preliminary numerical simulations demonstrated that FOS required by law and standards is maintained, and the proposed dam construction system complies with all ULS requirements as specified in EN 1997-1:2004. These simulations also indicated that the direct foundation is insufficient to ensure a displacement difference of less than 1.0 cm between successive segments of the conduit. Fig. 11 presents a contour plot of settlements at the foundation level of the dam, showing the settlement values caused by the dam construction across the entire width of the 3D numerical model.

Figure 11:

Contour plot of vertical displacement at the foundation level of the conduit pipe.

The maximum settlement of the dam subsoil within the conduit pipe, as determined by the computations, was 6.2 cm. Relative to the overall size of the structure, this value is not significant. However, as previously mentioned, the differential settlement between successive segments of the conduit pipe is critical from a structural standpoint. Fig. 12 presents a contour map of the displacement of the conduit pipe, showing both side and perspective views. The results are displayed on the deformed mesh.

Figure 12:

Displacement Y of conduit on deformed mesh: (a) perspective view (b) side view.

The maximum displacement difference between two successive segments was approximately 2.2 cm, occurring between segments 5 and 6. As a consequence, it was necessary to improve the ground to reduce the differential settlement between these segments.

The designed soil reinforcement consisted of 150 columns constructed using VDW technology, with a diameter of 880 mm (Sweco Consulting Sp. z o.o., 2020). The column lengths ranged from 8 to 12 m, with an average spacing of 2.8 m in plan. The length of the columns was determined based on the height of the overburden soil layers above the conduit and the depth of the marl layers with relatively uniform mechanical strength parameters (layers B and Cb2 from Tab. 2).

To verify the design solution for ground improvement, the columns were modeled within the dam model (Fig. 9) according to Approaches II and III. Approach I was not applied to this model for reasons discussed in Chapter 4.

Applying Approach II involved adding 150 beam elements embedded within the continuum. The beams were assigned parameters corresponding to a circular cross section with a diameter of 0.88 m, made of concrete. The lengths of the added beam elements ranged from 8 to 12 m, consistent with the actual dimensions of the columns. Approach III required the development of computational procedures described in Chapter 2.3. The first step involved determining the representative volume element (RVE). In the plan view, the RVE dimension corresponded to the column spacing, which is 2.8 m × 2.8 m, with a centrally located column of 0.88 m in diameter. Due to the variability of the subsoil layers and the changing thickness of these layers depending on the location, it was not possible to identify a single characteristic cross section to model the entire soil profile and then determine a single effective parameter for such a medium. Since the layer thickness was variable, it was decided to model a 1-m-high section of the column for each soil layer within the reinforcement area. Subsequently, in the reinforced subsoil area, the effective parameters obtained using homogenization theory tools were implemented in place of the original parameters presented in Tab. 1 and Tab. 2. The calculations for estimating effective parameters were performed using FlexPDE software. An example view of the model is shown in Fig. 13.

Figure 13:

An example view of the RVE model (yellow color – column, blue color – soil).

Following step III from Chapter 2.3, through appropriate integral operations, it was possible to estimate the effective parameters for each specific layer. Due to the numerous geotechnical layers in the subsoil, it was necessary to determine new parameters for each individual layer. Subsequently, the initial soil parameters were replaced with values calculated according to homogenization techniques. As in Section II, only the parameters for vertical stiffness were sought and assigned to the material as isotropic properties. Tab. 3 presents the resulting effective elastic parameters for the layers through which the column passes, along with a comparison to the baseline values from Tabs 1 and 2.

Please note that due to homogenization applied only for elastic parameters, only elastic parameters were modified within Approach III. Other strength parameters remained unchanged relative to the soil parameters without columns. According to the authors, this approach is conservative because introducing concrete columns with much higher strength than the surrounding soil should improve the effective strength parameters.

The results of the displacement contour maps are shown in Fig. 14 for Approach II and in Fig. 15 for Approach III.

Figure 14:

Displacement Y of the conduit on the deformed mesh for Approach II: (a) perspective view (b) side view.

Figure 15:

Displacement Y of the conduit on the deformed mesh for Approach III: (a) perspective view and (b) side view.

The results of the analyses confirmed the effectiveness of the design solutions adopted for soil reinforcement in both Approaches II and III. The maximum displacement differences between consecutive segments do not exceed 1.0 cm, which, as mentioned in Chapter 4, eliminated the risk of damage to the dilation joints between individual sections. For Approach II, the displacement difference was 0.9 cm, and for Approach III, it was 0.8 cm.

On analyzing the contour plots of displacement for two approaches utilized (Figs 14 and 15), some differences in both shape and value can be observed. The most important results obtained during the calculations are summarized in Tab. 1.

The results obtained show a significant reduction in the maximum settlements after the application of the columns (Tab. 4). The maximum vertical displacement of the conduit without the columns was 6.2 cm, while after the installation of the columns, it decreased to 3.7 cm for Approach II and 2.6 cm for Approach III. Regarding the maximum difference in displacements between successive segments, it can also be observed that the displacements decrease to an acceptable level of less than 1.0 cm after the installation of inclusions. The difference between Approaches II and III in this regard was negligible, at only 0.1 cm.

The gap in maximum displacement between approaches was mainly caused by differences in the stress state, which further affects the potential plasticity of the soil. For the elastic analysis and computational examples discussed in Chapters 2 and 3, the displacement results obtained by the different approaches were the same. Stress concentrations in Approach II (visible, e.g., in Fig. 6b) may cause the soil at the inclusion base to exhibit plastic flow, leading to increased settlements. In Approach III, the stress state is more uniform, which likely prevents such effects. In the authors’ opinion, modifying the soil parameters at the column bases may allow Approach III to account for these effects. Further research on such solutions is the subject of future work of the authors.

In Approach II, it is possible to reproduce local stress concentrations at the pile head, enabling the calculation of reinforcement in conduit pipe, for example, for the shear forces associated with the column. The internal forces in the column can also likely be obtained in Approach III by transferring the displacement field obtained on the face of the periodic cell to a small separate problem with a rigid inclusion modeled and then analyzing the column response.

Adding columns to the model according to both approaches does not cause a significant increase in computation time. Compared to the model without columns, the calculation time increased by no more than 10%.

Due to the built-in functions of ZSoil software, the process of modeling the columns according to Approach II is much easier than according to Approach III. Knowing the location of columns and their lengths allows for quick implementation. Approach III requires the use of more mathematical tools; however, it should be emphasized that once the tools are developed and defined parametrically, they can be reused for other types of calculations. The material substitution itself in the software is not problematic and can be even quicker than adding columns as beam elements.

The paper provides guidance for modeling the improvement of rigid inclusion ground in large-area geotechnical simulations. The following numerical modeling methods are discussed in this paper: modeling both column and soil as continuum elements, modeling the column as beam elements with soil as continuum elements, and replacing both soil and column with one equivalent medium modeled with continuum elements. The following conclusions can be formulated from the analyses: –

By modeling the stiff columns using the selected approaches, the displacements can be accurately represented, which is confirmed by results from reference tasks allowing for comparison of all approaches. Notably, in the case of elastic analyses, the results obtained with all the considered approaches are identical.