Particle Swarm Optimization of Interface Constitutive Model Parameters for Embedded Beam Formulations

Introduction

Background

Embedded Beam Formulations

Over a period of more than four decades, the development of EBFs has been motivated by the increasing demand for computationally efficient slender structure modeling techniques. Early work has been executed in the field of structural concrete, where EBFs have been found eligible for the modeling of reinforcing agents (Phillips and Zienkiewicz 1976; Chang et al. 1987). With regard to the field of computational geotechnics, EBFs have been first employed by Sadek and Shahrour (2004) for the analysis of micropiles. In comparison with the traditional standard FE approach (SFEA) shown in Fig. 2(a), where piles are explicitly discretized by means of 3D solid FEs and interface elements to account for SSI effects, EBFs offer several desirable features of practical relevance, listed as follows:

•

A more flexible mesh generation procedure with independent discretization of the 3D solid domain ( ${\mathrm{\Omega }}^s$) and the 1D beam domain ( ${\mathrm{\Omega }}^b$) representing the slender structure (Sadek and Shahrour 2004) (Fig. 1).

•

Simplified evaluation of stress resultants (e.g., bending moment and normal force) without posterior integration of stress field quantities (Granitzer and Felic 2022).

•

A reduced number of unknowns in the discrete system of equations, resulting in a significantly improved simulation runtime efficiency (Granitzer et al. 2024).

These relative merits render EBFs particularly suitable for iterative design procedures, such as described in Di Prisco et al. (2020) and Jürgens et al. (2022), which greatly benefit from the arbitrary positioning of beam nodes inside ${\mathrm{\Omega }}^s$, as opposed to the SFEA where mesh conformity constraints persist (Steinbrecher et al. 2020).

In view of the available methods for embedding ${\mathrm{\Omega }}^b$ in ${\mathrm{\Omega }}^s$, EBFs can be classified as embedded beam FEs with interaction line (EB-Ls) and embedded beam FEs with interaction surface (EB-Is) (Granitzer et al. 2024). In the pertinent literature, these approaches are also referred to as local and nonlocal embedding methods (Truty 2023), both of which are addressed in this work. Fig. 2(b) visually explains the key difference between both embedding methods based on a single-pile problem: EBFs classified as EB-Ls are equipped with a multipoint line-to-volume (1D-to-3D) coupling scheme, i.e., contact conditions between the beam nodes ( $B_k\in {\mathrm{\Omega }}^b$) and the surrounding target solid FE nodes ( $S_k\in {\mathrm{\Omega }}^s$) are formally described at multiple coupling points ${\mathbf{x}}_i$ placed along the centerline of ${\mathrm{\Omega }}^b$ (Tschuchnigg and Schweiger 2015). Assuming full adherence between ${\mathrm{\Omega }}^b$ and ${\mathrm{\Omega }}^s$ (Fig. 1) means that the motion of $B_k$ is fully controlled by the displacements of the central target 3D solid FE (TSFE) in which it is embedded [Fig. 2(b)]. This can be formally written as (Sadek and Shahrour 2004)

$${\mathbf{u}}_b(B_k)={\mathbf{u}}_s({\xi }_s^c)=\sum _{i=1}^{n_{S_k}}{N}_i^{\mathrm{TSFE}}({\xi }_s^c)\cdot {\mathbf{u}}_i^{\mathrm{TSFE}}$$

(1)

where ${\mathbf{u}}_b$ and ${\mathbf{u}}_s$ are the displacement vectors with colocal origin ${\mathbf{x}}_i$. $N_i^{\mathrm{TSFE}}$ and ${\mathbf{u}}_i^{\mathrm{TSFE}}$ are the element function and displacement vector, respectively, related to the $i$th node of the central TSFE that is equipped with $n_{S_k}$ solid nodes. ${\xi }_s^c$ is the local position vector of ${\mathbf{x}}_i$ described in the local coordinate system of the reference element associated with the central TSFE. To determine ${\xi }_s^c$, it is required to solve the following system of equations for all ${\mathbf{x}}_i$’s:

$${\mathbf{x}}_i-\sum _{i=1}^{n_{S_k}}{N}_i^{\mathrm{TSFE}}({\xi }_s^c){\mathbf{x}}_i^{\mathrm{TSFE}}=\mathbf{0}$$

(2)

where ${\mathbf{x}}_i^{\mathrm{TSFE}}=$ position vector of the $i$th TSFE node.

Recent research has paved the way to EB-I models that allow for capturing the coupling between ${\mathrm{\Omega }}^s$ and ${\mathrm{\Omega }}^b$ in a more general manner compared to EB-L methods (Turello et al. 2016; Smulders et al. 2019; Granitzer and Tschuchnigg 2021, 2023c). The basic idea of the underlying embedding method is to dispatch the coupling contributions onto a larger set of $S_k$ [cf. Truty (2023)]. Instead of employing geometrically reduced 1D-to-3D configurations, as incorporated by EB-Ls, EB-Is are equipped with a surface-to-volume (2D-to-3D) coupling scheme. In this way, the numerical SSI behavior is defined at multiple peripheral ${\mathbf{x}}_i$’s spread along the physical soil–structure interface [Fig. 2(b)]. Depending on the mesh topology, the coupling points are therefore enclosed by $n_{\mathrm{TSFE}}$ nonidentical peripheral TSFEs. As a consequence, the kinematical constraints have to be modified and incorporated in nonlocal form:

$${\mathbf{u}}_b(B_k)=\sum _{j=1}^{n_{\mathrm{TSFE}}}{w}_j\cdot \sum _{i=1}^{n_{S_k}}{N}_i^{\mathrm{TSFE}}({\xi }_s^p){\mathbf{u}}_i^{\mathrm{TSFE}}$$

(3)

where $w_j=$ global weight that accounts for the area fraction of the peripheral interaction surface segment ${\mathrm{\Gamma }}_{c,j}$ [Fig. 2(b)]; and ${\xi }_s^p=$ local position vector of ${\mathbf{x}}_i$ described in the local coordinate system of the respective peripheral TSFE. From Eq. (3), it becomes obvious that the motion of $B_k$ is coupled to the nodal displacements of multiple peripheral TSFEs, as opposed to EB-Ls. Favorable effects of this enhanced 2D-to-3D coupling scheme are highlighted in a number of publications (Turello et al. 2016; Smulders et al. 2019; Granitzer and Tschuchnigg 2021; Granitzer and Felic 2022; Granitzer and Tschuchnigg 2023c). A particularly relevant advantage of EB-Is is the reduced mesh size sensitivity. This can be explained by the distributed mobilization of coupling forces that reduce the influence of singularities on the numerical predictions. As will be demonstrated in this paper, this aspect has an important implication on the calibration of ICMPs.

For the sake of brevity, this subsection has focused on a comparative discussion of EB-Is and EB-Ls to facilitate understanding of the interface constitutive models, as well as the comparative study discussed later in this work. For a more comprehensive description of implementational details, including a mathematical description of alternative methods for the treatment of contact constraints between ${\mathrm{\Omega }}^b$ and ${\mathrm{\Omega }}^s$, interested readers can refer to Turello et al. (2018), Steinbrecher et al. (2022b), and Granitzer et al. (2024).

Numerical Description of Soil–Structure Contact

The soil–structure interface behavior of vertically loaded piles is characterized by displacement gradients that develop within a concentrated shear zone in the soil (Hu and Pu 2004; DeJong et al. 2006; Martinez et al. 2015). As shown in Fig. 3(a), pile displacements $u_p$ may thus be regarded as the sum of soil displacements $u_s$ (i.e., developing in the soil surrounding the local zone of shearing) and local shear displacements $u_{ls}$ (DeJong and Westgate 2009) [compare Lee and Xiao (2001)]. The latter are confined by particle size-dependent shear zones, which have a shear zone thickness of approximately $t_{\mathrm{sz}}=5\text{-}18.5\cdot D_{50}$, where $D_{50}$ is the mean grain size (Muir Wood 2002; Staubach et al. 2022). In this regard, it should be pointed out that this localized phenomenon forms within a zone spanning one FE in width (Ni et al. 2018b). To capture the shear localization with the finite-element method in a discrete manner, it is therefore required to use very thin solid FEs to represent the soil immediately adjacent to the interface (Potts and Zdravkovic 2001). Particularly in large-scale BVPs involving numerous piles (Tafili et al. 2023), however, this modeling technique leads to simulation runtimes that must be regarded as unbearable for most practical purposes (Ninić et al. 2014).

Alternatively, zero-thickness interface FEs could be used to overcome this problem (Ninić et al. 2014), as deployed in the SFEA [Fig. 3(b)]. This modeling approach aggregates the interface displacements $u_{\mathrm{int}}$, corresponding to $u_{ls}$, such that they become manifest in a boundary discontinuity where $t_{\mathrm{sz}}\to 0$. Particular advantages of this modeling approach are the ability to vary the interface constitutive behavior (e.g., the maximum friction angle) and the reduction of solid FEs with high mesh density close to the soil–structure interface leading to an improved simulation runtime. For this purpose, the SFEA is commonly used as a numerical benchmark solution to assess the credibility of novel EBFs; for example, see Engin et al. (2007), Tschuchnigg and Schweiger (2015), Turello et al. (2016), Granitzer and Tschuchnigg (2023a), and Truty (2023).

Contrary to the SFEA, EBFs represent the pile domain by means of dimensionally reduced 1D beam FEs and impose the coupling constraints between ${\mathrm{\Omega }}^s$ and ${\mathrm{\Omega }}^b$ in a pointwise manner at the implicit interaction domain (Fig. 2). In this regard, the readers should notice that there is no explicit geometry object in the 3D solid mesh to match the pile domain (Granitzer et al. 2024). While this pile modeling strategy proves beneficial in terms of modeling and computational efforts, it is unable to capture the discontinuities in the solid mesh. This limitation, also referred to as artificial continuity (Goudarzi and Simone 2019), applies for both embedding methods described in Figs. 3(c and d) and has important implications with respect to the numerical description of the soil–structure interface behavior (Granitzer and Tschuchnigg 2023c):

•

The accurate representation of local field variables (e.g., stresses and strains) within the local zone of shearing lies beyond the capabilities of EBFs.

•

Relative displacements $u_{\mathrm{rel}}$ developing at the interaction domain between ${\mathrm{\Omega }}^b$ and ${\mathrm{\Omega }}^s$ have a limited physical meaning since they are not manifested in the solid domain, as opposed to the SFEA where $u_{\mathrm{int}}$ can be numerically measured [Figs. 3(b–d)]. As a consequence, the representation of $u_{\mathrm{rel}}$ is, to some extent, restricted to the kinematic response of ${\mathrm{\Omega }}^b$. This concern applies for both embedding methods, albeit being considerably improved with the 2D-to-3D embedding method associated with EB-Is.

These theoretical considerations imply that the ICMPs of EBFs should be calibrated considering the global pile response (e.g., using load–settlement curves obtained from static pile load tests), instead of comparisons between experimentally measured and simulated relationships in stress–strain spaces, as typically carried out for the calibration of soil constitutive model parameters.

Nonlinear Interface Constitutive Model

Both EBFs considered in this work (EB-L* and EB-I*) are equipped with advanced interface constitutive models, built on the principle of displacement coupling (Steinbrecher et al. 2022a). Contrary to simplified EBFs that consider full adherence between ${\mathrm{\Omega }}^s$ and ${\mathrm{\Omega }}^b$, as shown in Fig. 1, they allow for capturing nonlinear effects of practical relevance, e.g., gap formation under tensile normal stresses and shear failure at the interaction domain. To realistically describe the mobilization of skin tractions along the shaft, the energy-conjugated pairs of relative displacements $u_{\mathrm{rel}}$ and skin tractions $t$ in the axial direction are related in a pointwise manner using empirical relationships. These are specified on the basis of extensive studies and have been found suitable for numerous problems of practical relevance (Tschuchnigg 2012). With regard to the EB-L*, the constitutive equations read (Tschuchnigg and Schweiger 2015)

$$t=\{\begin{array}{cc}50\cdot G_{el}^{av}\cdot {\mathrm{\Gamma }}_s\cdot u_{\mathrm{rel}},& \mathrm{\forall }{\mathbf{x}}_i\in {\mathrm{\Gamma }}_{c,s}\\ 50\cdot G_{el}^{av}\cdot {\mathrm{\Gamma }}_b\cdot R\cdot u_{\mathrm{rel}},& \mathrm{\forall }{\mathbf{x}}_i\in {\mathrm{\Gamma }}_{c,b}\end{array}$$

(4)

A similar approach has been adopted for the EB-I*, with the difference that the relations defined in Eq. (4) are additionally divided by the circumference $U^b$ and base area $A^b$ of the virtual beam geometry (Granitzer et al. 2024):

$$t=\{\begin{array}{cc}\frac{50\cdot G_{el}^{av}}{U^b}\cdot {\mathrm{\Gamma }}_s\cdot u_{\mathrm{rel}},& \mathrm{\forall }{\mathbf{x}}_i\in {\mathrm{\Gamma }}_{c,s}\\ \frac{50\cdot G_{el}^{av}\cdot R}{A^b}\cdot {\mathrm{\Gamma }}_b\cdot u_{\mathrm{rel}},& \mathrm{\forall }{\mathbf{x}}_i\in {\mathrm{\Gamma }}_{c,b}\end{array}$$

(5)

In Eqs. (4) and (5), $R$ is the virtual beam radius and $G_{el}^{av}$ is the average shear stiffness of the TSFEs belonging to the set of ${\mathbf{x}}_i$ located at the shaft ( ${\mathrm{\Gamma }}_{c,s}$) and base ( ${\mathrm{\Gamma }}_{c,b}$) interaction domain. It is important to note that the scalar-valued embedded interface stiffness multipliers ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ have been incorporated for calibration purposes. To the best of the authors’ knowledge, however, it is not publicly reported how these parameters affect the mechanical response of the EB-L* and EB-I*, respectively, and how the quality of numerical predictions can be improved by employing calibration procedures to select optimal values.

To shed light on these unsolved questions, Fig. 4 for the first provides insight into the influence of the embedded interface stiffness multipliers on the global response of an isolated pile modeled by means of the EB-L* and EB-I*. While the ultimate shaft resistance $R_s$, obeying to a Coulomb-type failure criterion (Tschuchnigg and Schweiger 2015; Smulders et al. 2019), remains relatively unaffected, ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ considerably influence the predicted pile resistance mobilization rates along the full range of loading. From this point of view, it appears imperative to select appropriate values for these ICMPs, mathematically expressed by the parameter vector $\mathbf{x}=[{\mathrm{\Gamma }}_s,{\mathrm{\Gamma }}_b{]}^T$, to approximate the benchmark pile response with high accuracy.

Provided that the ultimate shaft capacity $R_s$ is reasonably captured by the pile model, the aforementioned discussion renders a bivariate characterization problem (Cividini et al. 1981), which can be stated as follows: evaluate the values of the unknown components in $\mathbf{x}$ (i.e., ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$) such that the numerical predictions are as close as possible to the benchmark solution. For instance, the latter can be determined on the basis of field data (Granitzer et al. 2022), practical experience documented in the pertinent literature (German Geotechnical Society 2014), numerical benchmark simulations (Turello et al. 2016), or closed-form solutions (Randolph and Wroth 1978). As of now, however, default values $\mathbf{x}=[{\mathrm{\Gamma }}_s=1,{\mathrm{\Gamma }}_b=1{]}^T$ are considered due to the lack of numerical evidence. Note that this default configuration implies general applicability of Eqs. (4) and (5) for describing the SSI behavior of slender structures in numerical simulations, regardless of the embedding method employed. Nevertheless, due to the dimensional change between the embedding methods of the EB-L* (1D-to-3D) and EB-I* (2D-to-3D), as visually described in Fig. 2(b), convergence to the same optimal solution is not generally expected, i.e., $\{{\mathrm{\Gamma }}_{s,{\text{EB-L}}^{\ast }},{\mathrm{\Gamma }}_{b,{\text{EB-L}}^{\ast }}{\}}_{\mathrm{opt}}\nleftrightarrow \{{\mathrm{\Gamma }}_{s,{\text{EB-I}}^{\ast }},{\mathrm{\Gamma }}_{b,{\text{EB-I}}^{\ast }}{\}}_{\mathrm{opt}}$.

The aforementioned discussion raises obvious questions about the credibility of numerical simulations involving default values used for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ and the calibration of ICMPs associated with EBFs in general. Therefore, to the best of the authors’ knowledge, these aspects, along with the deployment of an automatic calibration procedure, are discussed for the first time in this work.

Numerical Tool Set

Working Scheme of the Optimal Parameter Vector Identification Strategy

As a starting point, let us regard the numerically predicted pile response obtained with the EB-L* or EB-I*, denoted by $R_{\mathrm{EBF}}$, as a function of the unknown parameter vector $\mathbf{x}$ (Fig. 4). From a mathematical point of view, finding the optimal parameter vector ${\mathbf{x}}_{\mathrm{opt}}=[{\mathrm{\Gamma }}_{s,\mathrm{opt}},{\mathrm{\Gamma }}_{b,\mathrm{opt}}{]}^T$ that minimizes the discrepancy between $R_{\mathrm{EBF}}$ and the user-defined benchmark results $R_{\mathrm{br}}$ can be expressed by

$$\underset{\mathbf{x}\in {\mathbb{R}}^2}{min}E(R_{\mathrm{br}},R_{\mathrm{EBF}})\text{s.t.}{R}_{\mathrm{EBF}}\to R_{\mathrm{br}}$$

(6)

where $E$ denotes the objective function that quantifies the discrepancy between $R_{\mathrm{br}}$ and $R_{\mathrm{EBF}}$. The proposed ACF solves this bivariate minimization problem by adopting the so-called direct approach (Cividini et al. 1981), i.e., ${\mathbf{x}}_{\mathrm{opt}}$ is identified by linking the results of parallel forward calculations to an optimization routine that iteratively updates the trial values in $X$. Employing direct procedures to solve related back-analysis problems for the determination of suitable numerical parameters in geotechnical simulations has previously found wide application, for example, in large-scale mass movements (Meier et al. 2008, 2013), deep foundations (Pitteloud and Meier 2019), and deep excavations (Calvello and Finno 2004; Lin et al. 2015); for a comprehensive overview of relevant geotechnical application cases, interested readers can refer to Ebid (2021). In all cases, the authors agree on one issue: Provided that the forward solution of the employed simulation model captures the relevant properties of the benchmark problem with reasonable accuracy, the direct approach yields a valuable concept to increase the numerical fidelity of finite-element simulations. With reference to this work, this applies also for the determination of suitable values for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$.

As visually explained in Fig. 6, the basic idea of the direct procedure adopted in this work concerns the identification of ${\mathbf{x}}_{\mathrm{opt}}$ by solving a user-defined reference BVP involving an EBF in an iterative manner, followed by the successive evaluation of the objective function for $\mathrm{\forall }\mathbf{x}\in X$ on the basis of Eq. (6). In the current version, the suitability of the candidate solutions is assessed by comparing $R_{\mathrm{EBF}}$ with $R_{\mathrm{br}}$, either in terms of a user-defined load–settlement curve or pile normal force distribution. Depending on the computed discrepancy between the benchmark data and the EBF simulation results, the trial values in $X$ are updated in parallel using a constrained optimization algorithm, i.e., ${\mathbf{x}}_{\mathrm{opt}}$ is bound to a predefined search space $\mathrm{\Omega }$ representing the range of admissible values for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$. This sequence of cycles continues until a termination criterion is fulfilled. With regard to the ACF configuration deployed in this work, the latter is defined by means of a user-defined maximum number of iteration cycles. Details concerning the mathematical expression of the objective function, optimization algorithm, and search space constraints are given in the next subsections.

Objective Function

The objective function driving the calibration of $\mathbf{x}$ has to quantify the discrepancy between $R_{\mathrm{br}}$ and $R_{\mathrm{EBF}}$ based on Eq. (6) (Fig. 6). Moreover, it provides knowledge about the calibration state. The readers should notice, however, that convergence to the global optimum solution is not generally warranted. For example, this may be attributed to an inappropriate selection of the similarity measure (Machaček et al. 2022) or the incorporation of reference data types with limited physical meaning (Kadlíček et al. 2022).

To circumvent related issues in the calibration process, the evaluation of the objective function carried out in this work is restricted to the load–settlement ( $F_z-u_z$) curve (Fig. 6), a widely accepted representation of the physical pile behavior. In addition, the developed ACF offers the functionality to specify two conceptually different similarity measures to quantify the discrepancy between $R_{\mathrm{br}}$ and $R_{\mathrm{EBF}}$, namely, the relative squared error ( $E_{\mathrm{RSE}}$) (Jadon et al. 2022) and the robust and interpolation-free technique ( $E_{\mathrm{RIFT}}$) (Lin et al. 2015). Both similarity measures are graphically described in Fig. 7(a). The dimensionless $E_{\mathrm{RSE}}$ constitutes a regression-type similarity measure that takes the general form:

$$E_{\mathrm{RSE}}=\frac{\sum _{i=1}^N(y_i-{\hat{y}}_i{)}^2}{\sum _{i=1}^N(y_i-\overline{y}{)}^2}$$

(7)

where

$$\overline{y}=\frac{1}{N}\sum _{i=1}^N{y}_i$$

(8)

In Eqs. (7) and (8), $y_i$ denotes the $i$th reference data point, ${\hat{y}}_i$ is the numerically predicted data point, and $\overline{y}$ is the mean of $N$ reference data points. Since the reference and numerical data points are not generally expected to be aligned in one dimension (e.g., if the pile resistance values are sampled at different settlements) and may differ with respect to the number of available data points, the evaluation of $E_{\mathrm{RSE}}$ requires the integration of an implicit regression procedure. The working principle of the underlying interpolation strategy adopted by the ACF is illustrated in Fig. 7(a). In this respect, it should be noted that the polynomial degree $p$ specified in the regression model should be specified with care because it may induce the common problem of underfitting or overfitting, leading to a poor approximation of the discrete data (Harrell et al. 1996) [Fig. 7(b)]. Alternatively, in cases where the discrete data describing the load–settlement curve is limited, it may be reasonable to deploy regularization techniques, such as Lasso regression (Ni et al. 2018a), to facilitate the evaluation of $E_{\mathrm{RSE}}$. However, it is pointed out that the use of related approaches may increase the complexity in the model selection due to the incorporation of additional hyperparameters, such as the shrinkage parameter associated with Lasso regression. In general, it is good practice to assess the suitability of the regression model employed in the ACF based on preliminary studies with varying $p$. In the scope of this work, $p=8$ was found to be adequate to achieve a reasonable performance of the calibration procedure.

As an attractive alternative to prevent the use of a regression model, $E_{\mathrm{RIFT}}$ has been implemented as alternative similarity measure. Instead of computing purely distance-based measures between (artificially) aligned data points, $E_{\mathrm{RIFT}}$ accounts for the accumulated area of successive triangles formed by the two data sets under study [Fig. 7(a)]. Since $E_{\mathrm{RIFT}}$ can be classified as an interpolation-free similarity measure, its application is particularly advised in cases where interpolation may introduce significant errors in the approximation of the candidate curve, for example, in the case of data sets with various levels of data density. To date, however, the performance of $E_{\mathrm{RIFT}}$ has been investigated on the basis of a rather small number of demonstration cases documented in two publications (Lin et al. 2015; Hsieh et al. 2017), both of which provide limited discussion about potential limitations. For this purpose, part of the comparative studies presented in this work focus on how the selection of the similarity measure affects the performance of the ACF in the context of a single-pile problem.

Optimization Method

The computational efficiency of the ACF is mainly governed by the number of forward calculations required to approximate ${\mathbf{x}}_{\mathrm{opt}}$ with sufficient accuracy. Accordingly, increasing the efficiency of the calibration procedure is equivalent to minimizing the number of iteration cycles (Fig. 6), notably at the cost of a higher parameter uncertainty. Consequently, it is imperative to select an optimization method that offers a reasonable balance between fast convergence and exploratory search, particularly because the objective functions may incorporate a number of locally optimal solutions. However, there is no general rule for selecting a proper optimization method since the performance depends on the problem characteristics that have to be assessed on a case-by-case basis (Meier and Schanz 2013).

In the present work, $\mathrm{\forall }\mathbf{x}\in X$ are iteratively updated employing the particle swarm optimization (PSO) algorithm, originally introduced by Kennedy and Eberhart (1995) and Shi and Eberhart (1998). The latter can be classified as the stochastic optimization algorithm (Wang et al. 2018), as opposed to deterministic algorithms that generally require gradient information with respect to the objective function; interested readers can refer to Hajihassani et al. (2018) for a comprehensive overview of PSO application cases in the field of geotechnical engineering, including hybrid configurations of PSO with conceptually different algorithms. Traditionally, PSO-based calibration frameworks greatly benefit from a high level of robustness, in the sense that broad variations in the hyperparameter selection do not prevent convergence (Clerc 2006), and quick exploration of large search spaces, possibly at the cost of a longer time consumption (Wang et al. 2018; Machaček et al. 2022). These features make the PSO algorithm particularly suitable for our low-dimensional characterization problem with two unknown calibration parameters ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$, which lacks profound numerical evidence about adequate value ranges for comparison.

Adopting the terminology originally introduced by Kennedy and Eberhart (1995), the PSO operates on the so-called swarm $X=\{{\mathbf{x}}^{(1)},\dots ,{\mathbf{x}}^{(n)}\}$ involving $n$ particles, whereas the embedded interface stiffness multipliers belonging to the $i$th particle ${\mathbf{x}}^{(i)}=[{\mathrm{\Gamma }}_s^{(i)},{\mathrm{\Gamma }}_b^{(i)}{]}^T\in \mathrm{\Omega }$ may be interpreted as a candidate solution to the minimization problem stated in Eq. (6). In the present work, the swarm size $n=20$ is defined following empirical guidelines documented in Shi and Eberhart (1998). Starting from random initial positions, a velocity update method is used to simultaneously transfer the individual ${\mathbf{x}}^{(i)}$ to updated positions at time step $t+1$ (Fig. 8):

$${\mathbf{x}}^{(i)}(t+1)={\mathbf{x}}^{(\mathbf{i})}(t)+{\mathbf{v}}^{(i)}(t+1)$$

(9)

where ${\mathbf{x}}^{(i)}(t)$ denotes the positions of the $i$th particle at time step $t$, corresponding to the last iteration cycle, and ${\mathbf{v}}^{(i)}(t+1)$ is the velocity of the $i$th particle at time step $t+1$, expressed by

$${\mathbf{v}}^{(i)}(t+1)=\underset{\text{cognitive component}}{\underbrace{c_1\cdot r_1\cdot ({\mathbf{x}}_p(t)-{\mathbf{x}}^{(i)}(t))}}+\underset{\text{social component}}{\underbrace{c_2\cdot r_2\cdot ({\mathbf{x}}_g(t)-{\mathbf{x}}^{(i)}(t))}}+\underset{\text{inertia component}}{\underbrace{w\cdot {\mathbf{v}}^{(i)}(t)}}$$

(10)

Therein, $c_1$ and $c_2$ are positive acceleration coefficients used to scale the contribution of the cognitive and social components, respectively. $r_1$ and $r_2$ denote random scalars between 0 and 1 sampled from a uniform distribution which introduce a stochastic element to the algorithm. ${\mathbf{x}}_p$ is the previous individual best position of the $i$th particle and ${\mathbf{x}}_g$ represents the previous global best positions of the swarm. $w$ is the inertia weight that controls the influence of the previous particle velocity ${\mathbf{v}}^{(i)}(t)$ in the course of exploratory search.

From Eq. (10), it can be inferred that the particle movement is governed by (1) its personal experience (cognitive component), (2) the experience of all particles in the swarm (social component) on the basis of the star topology [i.e., each particle disseminates information to all particles of the swarm; compare Ni and Deng (2013)], and (3) the inertia component in a role similar to friction in a physical setting [cf. Bratton and Kennedy (2007)]. In a biological context, this configuration emulates the social behavior of birds within a flock, particularly their ability to fly synchronously and continuously regroup in an optimal formation. A detailed discussion of the underlying social learning processes and the PSO hyperparameter characteristics, however, falls outside the scope of the present work. Details can be found in Engelbrecht (2007), Kennedy et al. (2009), and Helwig (2010).

In the current version of the ACF, the PSO hyperparameters documented in Table 1 are selected according to the specifications recommended in Eberhart and Shi (2000), originally derived on the basis of the constriction approach mathematically described in Clerc (1999) and Clerc and Kennedy (2002). The suitability of this hyperparameter set is additionally confirmed by empirical rules documented in Clerc (2006). Nevertheless, it should be noticed that the selection of PSO hyperparameters remains empirical to a large extent (Trelea 2003). This means that we may achieve a better optimization performance with a particular set of hyperparameter values. However, optimization experiments concerning the tuning of PSO hyperparameters, along with comparative studies employing alternative optimization algorithms, are beyond the scope of the present work and may be addressed in the future.

$c_1$ (1)	$c_2$ (1)	$w$ (1)
1.494	1.494	0.729

Search Space Constraints

The particle positions of the swarm are initialized uniformly at random in the search space $\mathrm{\Omega }\subseteq {\mathbb{R}}^2$ (Fig. 6). To handle invalid particle positions in the course of the optimization, the so-called reflect position handling strategy (Helwig 2010) was adopted, as conceptually explained in Fig. 8(a). Essentially, both features require the definition of suitable box constraints to limit the search space region, being a delicate issue given the limited experience in the selection of $\mathbf{x}$. In general, box constraint specifications have a considerable influence on the performance of PSO search behavior: While an overly large $\mathrm{\Omega }$ may lead to a high variance between the different candidate solutions, the specification of considerably narrow $\mathrm{\Omega }$ tends to produce clusters at the boundaries (Machaček et al. 2022). Needless to say, the constrained PSO algorithm is unable to evaluate appropriate solutions to Eq. (6) if ${\mathbf{x}}_{\mathrm{opt}}$ falls outside the feasible domain.

To evaluate suitable box constraints with adequate confidence, preliminary studies comprising 30 forward FE calculations with varying embedded interface stiffness multipliers were carried out employing the well-documented Alzey Bridge pile load test as a reference scenario (Sommer and Hambach 1974). Fig. 8(b) displays the adopted mesh configurations. Unless stated otherwise, the results are obtained with the medium mesh configuration comprising 20,964 tetrahedral solid FEs with quadratic element functions. The bored pile is modeled wished-in-place, that is, changes in the soil surrounding the pile induced by the installation process are not considered. The groundwater level is controlled by defined hydraulic pressure heads that are related to the groundwater table located 3.5 m below ground level. The numerical modeling parameters of the reference BVP, including the hardening soil small (HSS) (Benz 2007) constitutive model parameters of the soil domain conforming to slightly overconsolidated Frankfurt clay, combined with the material and geometrical parameters of the pile model are given in Table 2.

Category	Parameter	Value
Pile parameters (constitutive model: linear elastic)	$\gamma$ (kN/m ${}^3$)	25.0
	$E$ (MN/m ${}^2$)	10,000
	$R_{\mathrm{inter}}$ (1)	0.9
	$L$ (m)	9.5
	$R$ (m)	0.65
	$\nu$ (1)	0.2
Soil parameters (constitutive model: hardening soil small)	${\gamma }_{\mathrm{sat}}$, ${\gamma }_{\mathrm{unsat}}$ (kN/m ${}^3$)	20.0
	$E_{50}^{\mathrm{ref}}$ (kN/m ${}^2$)	45,000
	$E_{\mathrm{oed}}^{\mathrm{ref}}$ (kN/m ${}^2$)	27,150
	$E_{ur}^{\mathrm{ref}}$ (kN/m ${}^2$)	90,000
	$m$ (1)	1.0
	${\nu }_{ur}^{\prime }$ (1)	0.2
	${\phi }^{\prime }$ ( ${}^{\circ }$)	20.0
	$c$ (kN/m ${}^2$)	20.0
	POP (kN/m ${}^2$)	50.0
	$p_{\mathrm{ref}}$ (kN/m ${}^2$)	100.0
	$G_0^{\mathrm{ref}}$ (kN/m ${}^2$)	116,000
	${\gamma }_{0.7}$ (1)	0.00015

Source: Data from Granitzer and Tschuchnigg (2021).

All simulations were carried out utilizing the identical calculation phase sequence. The initial stress field was generated using the so-called K0 procedure. In this phase, only soil clusters are present and gravity loads are applied. These loads are derived from the soil unit weight and balanced by the effective vertical stresses in the soil. In analogy to Engin et al. (2007), where this pile problem has been adopted to validate an enhanced EB-L, the horizontal stresses are computed considering a preoverburden pressure $\mathrm{POP}=50$ kPa to realistically capture overconsolidation effects in the slightly overconsolidated Frankfurt clay. Interested readers can refer to Tschuchnigg (2012) for details concerning the influence of the adopted initial stress generation procedure on the yield surfaces of the HSS soil constitutive model. In the subsequent phases, the piles are activated and subjected to displacement-driven head-down loading.

A comparison of the experimental benchmark data obtained from field measurements (Sommer and Hambach 1974) and the load–settlement curves obtained using the EB-L* and EB-I* is plotted in Fig. 9. The readers should notice that the objective function $E_{\mathrm{RSE}}$, indicating the discrepancy between the results obtained from the measurements and the pile models, was evaluated by hand without using the ACF. From the results, three observations become striking: First, ${\mathrm{\Gamma }}_{b,\mathrm{opt}}$ evaluated for the EB-I* is considerably smaller than 1 (i.e., the default value) and differs by around one order of magnitude compared to results obtained with the EB-L*. This clearly demonstrates that optimal embedded interface stiffness multipliers cannot be interchangeably used for both implementations and need to be assessed on a case-by-case basis. Second, the EB-L* produces mesh-dependent values with relatively high ${\mathrm{\Gamma }}_{b,\mathrm{opt}}$ bandwidth, while the corresponding results obtained with the EB-I* remain almost constant in all cases considered. A third observation concerns the reproducibility of the benchmark solution: While the latter can be well captured by the EB-I*, irrespective of the mesh configuration, the EB-L* is unable to resemble the benchmark solution with the fine mesh discretization. An additional sensitivity study (not shown) has indicated similar but less pronounced tendencies with regard to ${\mathrm{\Gamma }}_{s,\mathrm{opt}}$.

To some extent, the latter two observations can be attributed to the unwanted mesh size dependency inherent to the EB-L*, which is almost eliminated with the EB-I* due to the enhanced coupling scheme. It should be emphasized that this notable difference in the mesh sensitivity has an important implication with regard to the general applicability of calibration procedures for the determination of ${\mathbf{x}}_{\mathrm{opt}}$: While the results indicate a high level of reproducibility with regard to the EB-I*, the mesh size effect on the results obtained with the EB-L* leads to calibrated parameter sets that are only valid for one particular mesh topology deployed in the reference BVP; that is, calibrated ICMPs associated with the EB-L* tend to be mesh size dependent. From a practical point of view, this concern greatly disqualifies the application of calibration strategies, such as the ACF, to improve the numerical fidelity of large-scale BVPs involving a high number of EB-Ls* and strategic mesh refinement zones. For this purpose, the EB-L* is not considered in further analyses. In contrast, the results obtained with the EB-I* allow one to infer mesh size-independent box constraints. These are summarized in Table 3 and adopted for the ACF demonstration cases presented in the next sections.

Parameter bounds	${\mathrm{\Gamma }}_s$ (1)	${\mathrm{\Gamma }}_b$ (1)
Min.	0.001	0.001
Max.	0.1	0.1

Reference Solution

The validation of the proposed ACF necessitates a careful evaluation of the reference solution ${\mathbf{x}}_{\mathrm{ref}}$ valid for the reference scenario. For this purpose, a grid search was adopted with evenly spaced intervals of the embedded interface stiffness multipliers. This involves taking a number of equally spaced trial solutions spread over the search space, thereby creating a total of 625 combinations to check. In Fig. 10(a), each candidate solution is represented by a black dot. As can be inferred from Fig. 10(b), the grid search-based reference solution with respect to $E_{\mathrm{RSE}}$ is ${\mathbf{x}}_{\mathrm{ref}}=[{\mathrm{\Gamma }}_{s,\mathrm{ref}}=0.005,{\mathrm{\Gamma }}_{b,\mathrm{ref}}=0.022{]}^T$. In this regard, it is interesting to note that the reference parameter set significantly differs from the default values. This observation clearly underlines that the choice of suitable values for the ICMPs is imperative to achieve a high level of numerical fidelity in simulations involving EBFs. In the next section, the grid search-based reference solution serves as a validation basis to explore the capabilities of the ACF to identify optimal values for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$.

Alzey Bridge Pile Load Test

Influence of Similarity Measure

The first set of ACF-enabled analyses is executed employing the two different similarity measures (RSE and RIFT). Fig. 11 studies the convergence behavior with regard to the mean error (i.e., the mean objective function value of the swarm) and the minimum error (i.e., the best particle solution) found within each iteration cycle of the automatic calibration. Since the similarity measures have different mathematical definitions, a quantitative comparison in light of their influence on the performance of the ACF is not possible. Instead, the similarity measure has to produce curves with stagnant decrease and approximate ${\mathbf{x}}_{\mathrm{ref}}$ with sufficient accuracy to be regarded as acceptable to solve the optimization task. In both cases, the convergence of the automatic calibration procedure is proven by a reduction of around 0.5–1.0 orders in the magnitude of the mean error. Moreover, it is demonstrated that the best solutions found by the ACF (i.e., ${\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RIFT}}$ and ${\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RSE}}$) after 12 and 19 iteration cycles are in remarkable agreement with the grid search-based reference solution ${\mathbf{x}}_{\mathrm{ref}}$ reported in Fig. 10(c).

The rate of convergence observed in Fig. 11 is largely controlled by the sensitivity of the objective function to the calibration parameters (i.e., ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$), whereas a higher level of sensitivity leads to a faster convergence. As an intuitive measure to quantify the parameter sensitivity and describe the relevance of each calibration parameter to achieve a good fit of the simulations with the measured load–settlement curve, it is common practice (Mendez et al. 2021) to compare the coefficient of variation (COV) with respect to the tested calibration parameters. This metric is defined as the ratio of the standard deviation normalized by the mean. As can be seen from the results documented in Table 4, the COV shows similar values in the range of 66.5%–74.6% for each embedded interface stiffness multiplier, irrespective of the similarity measure. This implies that the adequate estimation of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$, respectively, is of similar importance to achieve a good fit of the numerically derived EB-I* load–settlement curve with the measured reference data.

	${\mathrm{\Gamma }}_s$			${\mathrm{\Gamma }}_b$
Similarity measure	$\mu$ (1)	$\sigma$ (1)	COV (%)	$\mu$ (1)	$\sigma$ (1)	COV (%)
RIFT	0.043	0.028	66.5	0.040	0.028	69.9
RSE	0.039	0.029	74.6	0.036	0.026	72.3

From a user perspective, the most relevant aspect in the performance assessment of the ACF is probably the comparison of the measured load–settlement curve with the simulation results. In this context, Fig. 12 displays the effect of using calibrated values for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$. A visual inspection of the results obtained with ${\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RIFT}}$ and ${\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RSE}}$ confirms the significantly improved agreement with the measured load–settlement curve compared to the EB-I* simulations with noncalibrated (default) values for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$. From this observation, it becomes obvious that the numerical predictions obtained with the EB-I* greatly benefit from the ACF.

Despite the good agreement of EB-I* predictions obtained with calibrated ICMPs observed in Fig. 12, the computed COV values may be considered relatively high. This can be attributed to the wide limits selected for the (unknown) search space constraints (Table 3) and the lack of physical interpretation of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$, both of which represent fitting parameters that control the kinematic behavior of ${\mathrm{\Omega }}^b$ (see Fig. 4). This implies a certain degree of parameter uncertainty that is, to some extent, reflected in the nonnegligible variance of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ in the scatter plots illustrated in Fig. 13, showing the mutual particle distributions after 1, 10, and 20 iteration cycles. A lower scatter in the candidate solutions and hence an improved performance of the ACF in terms of reproducibility could be achieved by tuning the PSO hyperparameters in Table 1 or by selecting narrower bounds with respect to the search space constraints. The latter option, however, would require the examination of search space constraints on a case-by-case basis, representing a time-consuming task that impedes general employment of the ACF in practice. Considering the corresponding Gaussian kernel density estimates of the calibration parameters presented in Fig. 13, however, this appears nonessential because the results clearly indicate the desired gradual movement of particles toward ${\mathbf{x}}_{\mathrm{ref}}$. This concern is substantiated by the concentration of high-density regions close to ${\mathrm{\Gamma }}_{s,\mathrm{ref}}$ and ${\mathrm{\Gamma }}_{b,\mathrm{ref}}$ and holds for both similarity measures. Although the results infer the desired gradual increase of density peaks around ${\mathrm{\Gamma }}_{s,\mathrm{ref}}$ and ${\mathrm{\Gamma }}_{b,\mathrm{ref}}$ with increasing iterations cycle number and hence an adequate convergence of the calibration procedure toward the optimal solution, it appears somewhat premature to claim general applicability of the termination criterion adopted in this study. This aspect is therefore revisited in the next section considering an isolated pile in multilayer soil as a reference scenario.

The objective function topologies shown in Fig. 14 shed light on an additional interesting detail concerning the capabilities of the different similarity measures: While $E_{\mathrm{RSE}}$ shows a relatively smooth pattern with a global minimum close to ${\mathbf{x}}_{\mathrm{ref}}$, $E_{\mathrm{RIFT}}$ produces a clearly nondifferentiable function with multiple local minima. It could therefore be expected that the use of $E_{\mathrm{RIFT}}$ leads to a relatively reduced stability of the ACF. This can be explained by the different mathematical descriptions of the discrepancy between the measured load–settlement curve and the simulation results, particularly because $E_{\mathrm{RSE}}$ is more sensitive to large deviations between corresponding pairs of reference and numerically predicted data points. From the previous discussion, it can be concluded that the ACF is more prone to the occurrence of unwanted effects during the calibration procedure, such as premature convergence and local optimization (Qian et al. 2021), if $E_{\mathrm{RIFT}}$ is selected as a similarity measure. Consequently, automatic calibration by means of $E_{\mathrm{RSE}}$ is recommended and selected for further analyses.

Influence of Loading Magnitude

The results presented in the previous subsection are obtained considering the full load–settlement measuring interval $u_z\in [0,3]$ cm of the field test. Nevertheless, it needs to pointed out that this specification is user-defined and may vary depending on the available reference data. This step is comparable to the selection of experimental value ranges in oedometric and triaxial tests required for the calibration of constitutive model parameters (Knabe et al. 2013; Kadlíček et al. 2022). Considering the simplified numerical description of the soil–structure contact inherent to EBFs visually desribed in Fig. 3, in combination with the bivariate nature of the calibration problem underlying the ACF, it is nearly impossible that the EB-I* matches the measured load–settlement curve in a perfect manner. Even if the perfect interface constitutive model for the EB-I* existed, scattering of experimental results, measurement errors and remaining mesh size effects (albeit significantly reduced compared to the EB-L*) would prevent a perfect fit. From a calibration point of view, this is likely to result in nonunique solutions for ${\mathbf{x}}_{\mathrm{opt}}$ that are specific to the user provided load–settlement range. It is noteworthy that the accurate anticipation of relevant loading ranges for large-scale simulations, however, is subject to human influences, especially when dealing with practical problems at early project phases. This appears to be a step backward because it is the human factor that is pursued to become dispensable by the development of the ACF. Nevertheless, to this date, this path is inevitable [compare Machaček et al. (2022)].

To assess the impact of the selected loading conditions on the variability of ${\mathbf{x}}_{\mathrm{opt}}$, the automatic calibration procedure was executed using two loading magnitudes: $u_z/D_{\mathrm{pile}}=0.34\mathrm{\%}$, to capture the end of the initial linear region in the load–settlement curve (Chen and Fang 2009), and $u_z/D_{\mathrm{pile}}=1.0\mathrm{\%}$. As could be expected from the aforementioned discussion, the loading conditions influence the calibration results, with the optimal parameter vectors computed after 20 iteration cycles being ${\mathbf{x}}_{\mathrm{opt}}^{0.34\mathrm{\%}}=[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.006,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.038{]}^T$ and ${\mathbf{x}}_{\mathrm{opt}}^{1.0\mathrm{\%}}=$ $[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.007,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.015{]}^T$ (Fig. 15). Correspondingly, the loading magnitude has a notable effect on the objective function topology. Moreover, it can be seen that ${\mathrm{\Gamma }}_{s,\mathrm{opt}}$ yields similar values, as opposed to ${\mathrm{\Gamma }}_{b,\mathrm{opt}}$ values that are sensitive to the loading magnitude. This is attributed to the different contributions of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$, respectively, to the mobilization of EB-I* pile resistance (Fig. 4). As long as the ultimate shaft resistance reported at around $R_s=\mathrm{1,900}$ kN (Sommer and Hambach 1974) is not reached, the pile resistance mobilization rate of an EB-I* is dominated by ${\mathrm{\Gamma }}_s$. Consequently, ${\mathrm{\Gamma }}_b$ values specified for $u_z/D_{\mathrm{pile}}=0.34\mathrm{\%}$ are of subordinate importance to achieve a good fit with the measurements because this loading magnitude falls inside the initial loading range. In contrast, a notable base resistance is mobilized with $u_z/D_{\mathrm{pile}}=1.0\mathrm{\%}$. This has the consequence that the parameter uncertainty of ${\mathrm{\Gamma }}_b$ simultaneously decreases with increasing loading magnitude.

Despite the variance in the calibrated parameters triggered by a varying loading magnitude, Fig. 16 underlines that the ACF is a powerful tool to increase the credibility of simulations where the EB-I* is used for the modeling of pile-type structures. It should be noted, however, that the loading magnitude provided as input to the ACF may influence the calibration results, particularly with respect to the base resistance mobilization, i.e., overly small values may yield inaccurate estimates of ${\mathrm{\Gamma }}_{b,\mathrm{opt}}$. It is therefore advised to provide reference data with sufficiently high loading magnitude, such that the maximum loading levels experienced by the calibrated EB-I* in the envisaged large-scale BVP are captured.

Sony Center Pile Load Test

Finite-Element Model

All geometrical and constitutive model parameters utilized in this study were adopted from Wehnert and Vermeer (2004). Fig. 17 illustrates the discretized domain of the reference boundary value problem comprising 21,631 tetrahedral solid FEs with quadratic element functions. The soil stratification is assumed horizontally layered and consists of medium dense sand (Sand 1) underpinned by a dense sand (Sand 2) formation with a narrow interlayer of Marl; the corresponding hardening soil model (Schanz et al. 1999) and Mohr–Coulomb constitutive model parameters are listed in Table 5, along with the pile properties applied for the EB-I*. The calculation phase sequence is identical to the one utilized in the Alzey Bridge pile load test analysis scenario.

Category	Parameter	Value
Pile parameters (constitutive model: linear elastic)	$\gamma$ (kN/m ${}^3$)	25.0
	$E$ (MN/m ${}^2$)	10,000
	$R_{\mathrm{inter}}$ (1)	0.3/1.0
	$L$ (m)	23.0
	$R$ (m)	0.44
	$\nu$ (1)	0.2
Soil parameters Sand 1/Sand 2 (constitutive model: hardening soil)	${\gamma }_{\mathrm{sat}}$, ${\gamma }_{\mathrm{unsat}}$ (kN/m ${}^3$)	20.0/21.0
	$E_{50}^{\mathrm{ref}}$ (kN/m ${}^2$)	45,000/60,000
	$E_{\mathrm{oed}}^{\mathrm{ref}}$ (kN/m ${}^2$)	45,000/60,000
	$E_{ur}^{\mathrm{ref}}$ (kN/m ${}^2$)	180,000/270,000
	${\phi }^{\prime }$ ( ${}^{\circ }$)	32.5/37.5
	$\psi$ ( ${}^{\circ }$)	2.5/7.5
	$c$ (kN/m ${}^2$)	0.1
	$m$ (1)	0.6
	${\nu }_{ur}^{\prime }$ (1)	0.2
	$p_{\mathrm{ref}}$ (kN/m ${}^2$)	100.0
Soil parameters Marl (constitutive model: Mohr–Coulomb)	${\gamma }_{\mathrm{sat}}$, ${\gamma }_{\mathrm{unsat}}$ (kN/m ${}^3$)	21.0
	$E$ (kN/m ${}^2$)	200,000
	${\phi }^{\prime }$ ( ${}^{\circ }$)	27.5
	$c$ (kN/m ${}^2$)	30.0
	$\nu$ (1)	0.33

Source: Data from Wehnert and Vermeer (2004).

Influence of Iteration Cycle Number

As of now, there is little experience regarding the choice of the maximum iteration cycle number required to ensure satisfactory convergence of the optimal parameter vector to the sought solution. It should be noticed that this setting has a crucial impact on both, the time consumption of the calibration procedure and the quality of the identified calibration parameters. To provide insight into the relevance of the ACF termination criterion, the maximum iteration cycle number is increased to 50, as opposed to the previous examples where 20 iteration cycles are considered. Fig. 18(a) collects the main results of the automatic calibration as function of the iteration cycle number. In accordance with the numerical evidence obtained from the Alzey Bridge pile load test reference scenario, the results indicate reasonable convergence rates of the calibration procedure, underlined by a reduction of more than one magnitude of order in the mean error. Fig. 18(b) allows for a more detailed analysis with respect to the quality of the calibration parameters identified by the best particle position in the swarm (i.e., ${\mathbf{x}}_p$). It can be inferred that the first 25 iteration cycles yield a notable scatter in the objective function, followed by a relatively constant distribution until the termination criterion is triggered. In the present case, this means that the contribution of iteration cycle numbers higher than 25 to the increase in quality of the identified calibration parameters can be regarded insignificant.

From the results in Fig. 18(b), it can be also deduced that the best fit of the EB-I* simulations with the measured load–settlement curve, corresponding to the particle position with the smallest $E_{\mathrm{RSE}}$ value found during the iterative search (i.e., ${\mathbf{x}}_{\mathrm{opt}}=[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.005,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.018{]}^T$), is qualitatively comparable to the results obtained with ${\mathbf{x}}_p$ in Iteration cycle 6. This is particularly reflected by the sharp drop of the minimum error in Fig. 18(a). The readers should notice, however, that this observation is largely triggered by the initialization of the particle positions and velocities and hence should not be generalized to reduce the runtime of the automatic calibration because the final solution quality may be affected in an adverse manner [cf. Helwig (2010)]. It is therefore essential to notice that the formulation of general guidelines regarding the maximum number of iteration cycles is out of reach for this work because it depends on numerous of factors, including but not restricted to the initial configuration of the swarm, box constraints, similarity measures, and reference data. In addition, it may be advised to introduce a supplemental termination criterion that defines a minimum lowering of $E_{\mathrm{RSE}}$ required after a number of iteration cycles to proceed with the iterative search (Meier et al. 2008). To date, we can only suggest to produce a study similar to the one depicted in Fig. 18, in combination with a detailed visual inspection of critical areas in the load–settlement curves to ensure calibration parameters with sufficient accuracy at reasonable computational expense. In the present case, the numerical predictions depicted in Fig. 19 show that the simulated EB-I* load–settlement curve is considerably improved between Iteration cycles 1 and 20, particularly in the initial loading range. Inversely, higher iteration cycles tend to yield an insignificant change in the predicted EB-I* response but affect the computational expense associated with the automatic calibration in an unfavorable manner.

Pile Group Analysis

Investigated Scenario

Fig. 20 illustrates the problem geometry and loading conditions adopted from Franza and Sheil (2021). In their work, a stick-slip EBF (Granitzer et al. 2024) has been used for the pile modeling, notably without providing details about the interaction domain geometry and relevant ICMPs that control the pile displacement behavior prior to the onset of plastic slip. Table 6 lists the linear–elastic pile and HSS soil parameters adopted for Vienna sand, originally published by Tschuchnigg and Schweiger (2010). In analogy to Franza and Sheil (2021), the ACF performance was examined based on comparisons with a numerical benchmark solution where the piles are modeled by means of the SFEA, in which zero-thickness interface elements (Day and Potts 1994) are introduced at the pile–soil interfaces with the same properties of the soil. From the displayed mesh topologies, it can be observed that the FE number and hence the number of unknowns that has to be solved are substantially reduced with the EB-I* model. The elevated rigid cap is modeled weightless and attached to the piles by means of fixed connections. It is pointed out that the cap is not in contact with the underlying soil. This geometrical configuration renders the BVP particularly suitable for examining the ACF performance because the structural response is fully controlled by interactions between the piles and the soil that are, in turn, bound to the values selected for ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$. The calculation phase sequence is identical to the ones utilized in the previous analysis scenarios. A salient difference is the selected point displacement magnitude imposed on the cap center $u_z/D_{\mathrm{pile}}=10\mathrm{\%}$ that is defined in accordance with design guidelines commonly used in pile group analyses (Witt 2018).

Category	Parameter	Value
Pile parameters (constitutive model: linear elastic)	$\gamma$ (kN/m ${}^3$)	25.0
	$E$ (MN/m ${}^2$)	10,000
	$R_{\mathrm{inter}}$ (1)	0.9
	$L$ (m)	12.5
	$R$ (m)	0.25
	$\nu$ (1)	0.2
Soil parameters (constitutive model: hardening soil small)	${\gamma }_{\mathrm{sat}}$, ${\gamma }_{\mathrm{unsat}}s$ (kN/m ${}^3$)	20.0
	$E_{50}^{\mathrm{ref}}$ (kN/m ${}^2$)	25,000
	$E_{\mathrm{oed}}^{\mathrm{ref}}$ (kN/m ${}^2$)	25,000
	$E_{ur}^{\mathrm{ref}}$ (kN/m ${}^2$)	62,500
	$m$ (1)	0.65
	${\nu }_{ur}^{\prime }$ (1)	0.2
	${\phi }^{\prime }$ ( ${}^{\circ }$)	32.5
	$\psi$ ( ${}^{\circ }$)	2.5
	$c$ (kN/m ${}^2$)	2.0
	POP (kN/m ${}^2$)	600.0
	$p_{\mathrm{ref}}$ (kN/m ${}^2$)	100.0
	$G_0^{\mathrm{ref}}$ (kN/m ${}^2$)	78,125
	${\gamma }_{0.7}$ (1)	0.0002

Source: Data from Tschuchnigg (2012).

Parameter Calibration

In accordance with the studies presented in the two preceding sections, ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ are calibrated based on a single-pile reference BVP (Fig. 21). The constitutive model parameters, pile geometry and loading conditions, are identical to the piled cap reference scenario documented in Fig. 20 and Table 6. In the absence of field measurements, however, the benchmark load–settlement curve was recovered from the numerical SFEA benchmark model. The PSO hyperparameters and search space constraints considered in the ACF were adopted from Tables 1 and 3, respectively. To ensure a reasonable convergence of the automatic calibration, $E_{\mathrm{RSE}}$ was used as a similarity measure.

As a key output of the calibration procedure, Fig. 21 shows the EB-I* results evaluated with the best parameter vector ${\mathbf{x}}_{\mathrm{opt}}=[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.010,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.012{]}^T$ observed during 20 iteration cycles and the load–settlement curve of the numerical benchmark model. A comparison with the EB-I* results corresponding to the worst parameter vector ${\mathbf{x}}_{\mathrm{worst}}=$ $[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.026,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.010{]}^T$ highlights, on the one hand, the remarkable ability of the EB-I* to emulate the SFEA benchmark response, provided that adequate IMCPs are selected. On the other hand, it confirms the high potential of the adopted ACF to increase the credibility of EB-I* predictions by solving this calibration task. In the next subsection, we investigate whether these observations also extend to the piled cap problem involving numerous piles.

Foundation Displacements

To date, the majority of simulations involving EB-Is* reported in the literature has been constrained to single-pile analyses. In this study, the EB-I* is for the first time used to analyze a capped pile group (Fig. 20). To assess the significance of using calibrated ICMPs to achieve EB-I* predictions with reasonable accuracy, the simulations were carried out with default and calibrated parameter sets, that is, $\mathbf{x}=[{\mathrm{\Gamma }}_s=1,{\mathrm{\Gamma }}_b=1{]}^T$ and ${\mathbf{x}}_{\mathrm{opt}}=[{\mathrm{\Gamma }}_{s,\mathrm{opt}}=0.010,{\mathrm{\Gamma }}_{b,\mathrm{opt}}=0.012{]}^T$.

Fig. 22 illustrates the pile load–settlement curves obtained with the EB-I* models and the numerical SFEA benchmark. It can be inferred from the results that the load-bearing behavior of the piles differs with the respective pile position, whereas the corner pile displays the highest pile resistance mobilization rates and the central pile shows the least. This observed tendency is typical for low-settlement bored pile groups and becomes more determinant with increasing pile spacing, pile number, pile slenderness ratio ( $L/D$), and soil stiffness; interested readers can refer to Poulos and Davis (1980) for details. Although not observed in this study, it should be noted for completeness that settlement magnitudes $u_z/D_{\mathrm{pile}}>10\mathrm{\%}$ could lead to a reversal of this tendency due to interlocking effects. The results presented in Fig. 22 additionally highlight the high relevance of an adequate selection of ICMPs. While the default parameter set results in an overly stiff response of the piles, the remarkable agreement of the calibrated EB-I* results with the numerical SFEA benchmark observed along the full range of loading is striking.

With reference to the selection of the reference BVP in Fig. 6, the ACF could also be theoretically deployed using a full representation of the piled cap, leading to a considerable increase in computational expense for the calibration procedure. However, the results presented in this study underline the generalizability of the ACF in identifying optimal parameter sets based on single-pile reference BVPs, in the sense that they are equally applicable to large-scale BVPs involving multiple EB-Is*. Readers should note that this is an important observation of this study because it allows one to restrain the computational effort during the calibration procedure to stay beyond practical limits.

Conclusions

Data Availability Statement

Acknowledgments

Notation

$A^b,U^b$

base area (m ${}^2$) and circumference (m) of virtual beam geometry;

$B^k,S^k$

beam nodes and virtual solid nodes;

$c_1,c_1,w$

acceleration coefficients and inertia weight (PSO hyperparameters);

$D_{\mathrm{pile}}$

pile diameter (m);

$d_{50}$

medium value of the particle size distribution (m);

$E_{\mathrm{RIFT}}$

error obtained with robust and interpolation-free technique;

$E_{\mathrm{RSE}}$

relative squared error;

EB-L

embedded beam FE with interaction line (1D-to-3D coupling scheme);

EB-L ${}^{\ast }$

EB-L after Tschuchnigg and Schweiger (2015);

EB-I

embedded beam FE with interaction surface (2D-to-3D coupling scheme);

EB-I ${}^{\ast }$

EB-I after Granitzer et al. (2024);

$F_z$

vertical pile resistance (kN);

$G_{el}^{av}$

average shear stiffness of TSFE (kN/m ${}^2$);

$N_i^{\mathrm{TSFE}}$

element function associated with the $i$th node of TSFE;

$n_{S_k}$

number of solid nodes belonging to TSFE;

$n_{\mathrm{TSFE}}$

number of TSFE at the EB-I interaction surface;

$R$

virtual beam radius (m);

$R_{\mathrm{br}},R_{\mathrm{EBF}}$

benchmark response and numerically predicted response;

$R_s$

ultimate shaft resistance (kN);

$r_1,r_2$

random scalar values between 0 and 1;

$t$

traction component in the axial direction (kN/m, kN/m ${}^2$);

$t_{\mathrm{sz}}$

shear zone thickness (m);

${\mathbf{u}}_b$, ${\mathbf{u}}_s$, ${\mathbf{u}}_i^{\mathrm{TSFE}}$

displacement vectors of beam and solid domain/ $i$th node of TSFE (m);

$u_{\mathrm{int}},u_{ls},u_p$

interface, local shear, and pile displacements (m), respectively;

$u_{\mathrm{rel}},u_s$

relative and soil displacements (m), respectively;

$u_z$

vertical pile head settlement (m);

${\mathbf{v}}^{(i)}(t),{\mathbf{v}}^{(i)}(t+1)$

velocities of the $i$th particle at time steps $t$ and $t+1$, respectively;

$w_j$

global weight of the peripheral interaction surface segment;

$\mathbf{x}$

parameter vector with values of embedded interface stiffness multipliers;

${\mathbf{x}}^{(i)}(t),{\mathbf{x}}^{(i)}(t+1)$

positions of the $i$th particle at time steps $t$ and $t+1$, respectively;

${\mathbf{x}}_g,{\mathbf{x}}_p$

global and individual best particle positions, respectively;

${\mathbf{x}}_i$

position vector of the $i$th coupling point (m);

${\mathbf{x}}_i^{\mathrm{TSFE}}$

position vector of the $i$th node belonging to TSFE (m);

${\mathbf{x}}_{\mathrm{opt}}$

optimal parameter vectors with calibrated values of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$ identified by ACF;

${\mathbf{x}}_{\mathrm{opt}}^{0.34\mathrm{\%}},{\mathbf{x}}_{\mathrm{opt}}^{1.0\mathrm{\%}}$

optimal parameter vector for $u_z/D_{\mathrm{pile}}=0.34\mathrm{\%}$ and $u_z/D_{\mathrm{pile}}=1.0\mathrm{\%}$, respectively;

${\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RIFT}},{\mathbf{x}}_{\mathrm{opt}}^{\mathrm{RSE}}$

optimal parameter vectors considering RIFT and RSE, respectively;

${\mathbf{x}}_{\mathrm{ref}}$

grid search-based parameter vector with reference values of ${\mathrm{\Gamma }}_s$ and ${\mathrm{\Gamma }}_b$;

${\mathrm{\Gamma }}_b,{\mathrm{\Gamma }}_s$

embedded interface stiffness multipliers considered at the base and shaft, respectively;

${\mathrm{\Gamma }}_{c,b},{\mathrm{\Gamma }}_{c,s}$

interaction zones at the base (b) and shaft (s), respectively;

${\mathrm{\Gamma }}_{c,j}$

peripheral interaction surface segment;

${\mathrm{\Gamma }}_{(\cdot ),{\text{EB-L}}^{\ast }},{\mathrm{\Gamma }}_{(\cdot ),{\text{EB-I}}^{\ast }}$

embedded interface stiffness multipliers associated with EB-L* and EB-I*, respectively;

$\mu$

mean value;

${\xi }_s^c,{\xi }_s^p$

position vectors in the local reference system of TSFE;

$\sigma$

standard deviation;

$\mathrm{\Omega }$

search space; and

${\mathrm{\Omega }}^s,{\mathrm{\Omega }}^b$

3D solid (s) and 1D beam (b) domains, respectively.

References