**Contents**


#### **Tomislav Kiˇsiˇcek, Tvrtko Reni´c, Damir Lazarevi´c and Ivan Hafner**


## *Article* **Shear Strength Estimation of Reinforced Concrete Deep Beams Using a Novel Hybrid Metaheuristic Optimized SVR Models**

**Mosbeh R. Kaloop 1,2,3 , Bishwajit Roy 4 , Kuldeep Chaurasia 4 , Sean-Mi Kim <sup>1</sup> , Hee-Myung Jang 1 , Jong-Wan Hu 1,2,\* and Basem S. Abdelwahed 5**


**Abstract:** This study looks to propose a hybrid soft computing approach that can be used to accurately estimate the shear strength of reinforced concrete (RC) deep beams. Support vector regression (SVR) is integrated with three novel metaheuristic optimization algorithms: African Vultures optimization algorithm (AVOA), particle swarm optimization (PSO), and Harris Hawks optimization (HHO). The proposed models, SVR-AVOA, -PSO, and -HHO, are designed and compared to reference existing models. Multi variables are used and evaluated to model and evaluate the deep beam's shear strength, and the sensitivity of the selected variables in modeling the shear strength is assessed. The results indicate that the SVR-AVOA outperforms other proposed and existing models for the shear strength prediction. The mean absolute error of SVR-AVOA, SVR-PSO, and SVR-HHO are 43.17 kN, 44.09 kN, and 106.95 kN, respectively. The SVR-AVOA can be used as a soft computing technique to estimate the shear strength of the RC deep beam with a maximum error of ±3.39%. Furthermore, the sensitivity analysis shows that the deep beam's key parameters (shear span to depth ratio, web reinforcement's yield strength, concrete compressive strength, stirrups spacing, and the main longitudinal bars reinforcement ratio) are efficiently impacted in the shear strength detection of RC deep beam.

**Keywords:** reinforced concrete; deep beam; shear strength; support vector regression; metaheuristic optimization

#### **1. Introduction**

In many high-rise reinforced concrete (RC) buildings, as the use of areas is changed from one story to another, some columns in the upper stories are not permitted to reach the foundation. To solve this conflict, transfer girders with a considerable thickness named deep beams are required [1,2]. Furthermore, deep beams are used in many other critical structures and play a significant role in delivering heavy loads to the bearing elements [2,3]. Deep beams have high flexural stiffness, and the shear diagonal failure is the predominant mechanism as the loads are mainly transferred from their action points to the supports locations through a direct diagonal strut [1–4]. Concrete compressive strength, the provided top/bottom reinforcements, and web reinforcements in terms of amount and spacing all form the shear resistance of these deep beams [1,2,4,5]. In literature, plenty of analytical and numerical studies have been focused on the ultimate shear strength assessment of such beams with large depths compared to their spans [2,3,6–9]. Unavoidable discrepancies were found with the implementation of both analytical/numerical methods compared to the experimental results [4]. This study aims to propose a novel soft computing approach

**Citation:** Kaloop, M.R.; Roy, B.; Chaurasia, K.; Kim, S.-M.; Jang, H.-M.; Hu, J.-W.; Abdelwahed, B.S. Shear Strength Estimation of Reinforced Concrete Deep Beams Using a Novel Hybrid Metaheuristic Optimized SVR Models. *Sustainability* **2022**, *14*, 5238. https://doi.org/10.3390/su14095238

Academic Editors: Chiara Bedon, Maged A. Youssef, Mislav Stepinac, Marco Fasan and Ajitanshu Vedrtnam

Received: 22 March 2022 Accepted: 20 April 2022 Published: 26 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

that can be used to predict an accurate shear strength of RC deep beams based on a wide range of test results collected from different experimental studies.

Many researchers have used different regression methods to estimate the shear strength of RC deep beams [8,10]. Recently, soft computing techniques have been proposed and improved the prediction techniques of shear strength of beams [3,7,11–15]. A conventional artificial neural network (ANN) was applied to estimate the shear strength of the RC deep beam, and the accuracy of the designed model was high [7]. The ultimate shear strength of the RC deep beam was also estimated using ANN and compared to different building codes, and the proposed model provided an accurate prediction of shear capacity [6]. ANN, adaptive network-based fuzzy inference system (ANFIS), and group method of data handling (GMDH) approaches were used in predicting the shear strength of RC beam-column joints, and the performance of these models was high at a different range of shear strength [15]. A probabilistic model was applied to estimate the shear strength of beams, and the determination of shear strength was shown to be accurate [14]. Integrated genetic programming and simulated annealing (GSA) outperformed American concrete institute (ACI) and Canadian standard association (CSA) codes in modeling the shear strength of RC deep beams [2]. The shear strength of beams reinforced by fiber was calculated using hybrid support vector regression (SVR) and firefly optimization algorithm (FFA), and the designed model was shown to be robust in shear strength prediction [13]. ANN was integrated with the adaptive harmony search optimization (AHS) technique for modeling the shear strength of RC walls, and the proposed model accuracy was high [12]. Multivariate adaptive regression splines (MARS) and artificial bee colony (ABC) were also integrated to design a model for predicting the shear strength of RC deep beams, and the performance of MARS-ABC was higher than different building codes in shear strength estimation [16]. Generally, parameters of machine learning (ML) models are tuned using metaheuristic algorithms to improve the prediction efficiency of ML models [17–22].

Meanwhile, novel optimization algorithms have recently been developed, such as the African Vultures optimization algorithm (AVOA), particle swarm optimization (PSO), and Harris Hawks optimization (HHO). Although these techniques are used in different engineering applications [17–19], the AVOA optimization technique is not applied yet in shear strength prediction based on our literature. PSO was integrated with an adaptive neuro-fuzzy inference system (ANFIS) to predict the shear strength of high strength concrete for a slender beam, and the ANFIS-PSO attained the best modeling accuracy over ANFIS -ant colony optimizer (ANFIS-ACO), -differential evolution (ANFIS-DE), and -genetic algorithm (ANFIS-GA) [20]. Teaching–learning-based optimization (TLBO), PSO, and HHO were integrated with SVR, and the results of SVR-PSO, SVR-HHO, and SVR-TLBO were robust and can be used to estimate an accurate shear strength prediction of RC shear walls [21].

Based on the above literature, the hybrid SVR models are more robust for predicting the shear strength of RC deep beams [3,16,21,23]. This study aims to evaluate a new hybrid technique (SVR-AVOA) in predicting the shear strength of RC deep beams. To benchmark the proposed SVR-AVOA model, the hybrid known models SVR-PSO and SVR-HHO are proposed and compared; in addition, the recent mathematical studies and building codes are compared to the proposed model to assess its accuracy of it in modeling shear strength of RC beam. SVR-AVOA, SVR-PSO, and SVR-HHO are developed using different scenarios of input variables. For this study, 202 datasets, including 19 variables of experimental studies, were collected from literature to design and evaluate the proposed models. The sensitivity analysis of optimum input variables is proposed and evaluated.

#### **2. Background of Variables Impacts the Shear Strength of RC Deep Beams**

Figure 1i presents a real case of deep beam function in load transfer of buildings and variables that impact the shear strength value. Figure 1ii demonstrates the main parameters of the deep beam. Figure 1iii illustrates the failure mode of deep beams. In the figures, *V* represents the deep beam shear capacity, *a* is the horizontal distance from the load to the support, and *d* denotes the deep beam depth. Here, *V* depends on (1) Concrete quality *f* ′ *c* (for the diagonal strut), (2) Main steel yield strength *f<sup>y</sup>* (for the main tie), and (3) Web reinforcement (horizontal and vertical). As some variables have a big role in forming the deep beam's shear strength, the main variables considered in this study are the shear span to depth ratio, the main reinforcement, ratio and yield strength, concrete compressive strength, and web reinforcement characteristics.

**Figure 1.** (**i**) Real case of using deep beams (**left**) and deep beam terminology (**right**); *V*: Shear strength, *a*: Shear span, *d*: Effective depth of beam. (**ii**) Basic reinforcement details of simple RC deep beam. (**iii**) Failure pattern in deep beam (**left**), Different mechanism components (**middle**), and Flow of forces in deep beam (**right**); where, *C* = compression force in concrete, *T* = tensile force in the main steel, *Hs* = tensile force in horizontal web reinforcement, *Vs* = tensile force in vertical web reinforcement; 1: Main diagonal strut, 2: Splitting tension force, 3: Main tensile force.

Many researchers investigated the deep beam's shear strength evaluation and predictions [1–4,24,25]. They concluded that the compressive strength of concrete (*f* ′ *c* ), the shear span to the beam's depth ratio (*a*/*d*), bottom longitudinal reinforcement ratio (*ρ*), and web reinforcement ratio (both vertical *ρv* and horizontal *ρh*) are the main key parameters. The previous analytical/experimental studies focused on the role of *f* ′ *c* in forming *V*. In the case of beams with smaller *a*/*d* ratios, the role of the truss mechanism in transferring loads to the support location diminishes, and the direct diagonal strut is the main transferring load mechanism to the support's location. The effectiveness of such struts has a significant impact on *V* values [3]. Smith and Vantsiotis [26] and Ahmed [27] observed an increase in *V* of the deep beam with increasing fc of concrete, but the relationship was not linearly proportional. In addition, they observed no improvement in *V* if *f<sup>c</sup>* was above a certain limit. The non-proportional increase in shear strength compared to the increase in concrete compressive strength can be attributed to two reasons. First, the limited contribution of the aggregate interlock mechanism in members with high strength concrete compared to the one with normal strength concrete, as the cracks cross the aggregate particles in high strength concrete and do not go around them as in normal strength concrete. Second, the formed tensile strains perpendicular to the main diagonal strut work on reducing the benefits of using high strengths. Oh and Shin [28] noticed a brittle failure of deep beams with concrete of 74 MPa without any warning, which is different from the failure of other beams with 23 MPa. They also observed a decrease in the rate of increase in the ultimate strength of beams with high-strength concrete.

The inclination angle of the main diagonal strut plays a significant role in determining the concrete efficiency in the diagonal strut. This angle is directly dependent on the shear span to depth ratio *a*/*d*. As this angle increases, the forces can go directly inside the diagonal strut to the support. The previous studies [3,29] noticed that the increase in shear strength could be detected by decreasing the *a*/*d* ratio. Kim and Park [30] found a trivial impact of this ratio on the shear strength of beams with ratios greater than three, and the contrary was noticed for beams with ratios less than three. Oh and Shin [28] concluded that the ratio of *a*/*d* is the governing key parameter in determining the shear strength of a deep beam with high-strength concrete. In addition to resisting the induced tensile force of the main horizontal tie, the main reinforcement bars play an important role in controlling the width enlargement of the diagonal main cracks by dowel action mechanism and enable the aggregate interlock to work more effectively. Many researchers investigated the impact of the main reinforcement ratio on the deep beam's shear strength [3,31,32]. They observed a significant increase in the shear strength with increasing the main reinforcement ratio but up to a certain limit. Above a ratio of 1.5%, Ashour et al. [33] noticed a local concrete crushing damage due to compressive stress concentration at the top strut far from the main diagonal strut without enabling the diagonal strut to reach its ultimate resistance.

Web reinforcement has an important role in confining the concrete and delivering the tensile stresses at the main shear diagonal crack to the intact zones around the crack, which consequently increases the deep beam's shear strength [2]. Both vertical and horizontal web reinforcement has a key role in resisting shear stresses and limiting the enlargement of the width of the main diagonal crack [3]. In deep beams with higher *a*/*d* ratios, the contribution of vertical reinforcement is more obvious than the horizontal reinforcements, and the contrary is noticed for beams with *a*/*d* less than 1.0. As the deep beam's shear strength is dependent on many parameters, the increase of horizontal and vertical reinforcement above a certain limit does not influence the ultimate shear strength as other parameters may govern the situation without reaching the maximum capacity of the provided web reinforcement.

Table 1 presents the existing models used in this study compared to the developed models. The performance of these models was used significantly in shear strength determination for the RC deep beams of structures.


where: *f* ′ *c* is the compressive strength of concrete, *b* is the beam width, *d* is the beam effective depth, *Av* is the vertical web reinforcement, *fyv* and *fyh* are the yield strength of vertical and horizontal web reinforcement respectively, *s* is the spacing between the vertical web reinforcement, *ρ<sup>h</sup>* and *ρ<sup>v</sup>* are the ratio of horizontal and vertical web reinforcement respectively, *n* is the modular ratio, *ρ* is the ratio of the main longitudinal bars.

#### **3. Material and Data Collection**

The Supplementary Material (Table S1) presents the data collected from the literature. For this study, 202 datasets were collected from the literature. In the current study, 19 input variables were used and divided into two categories, the main (8 variables) and other (11 variables) variables, as presented in Table S2. Here, the main variables were considered based on our literature in the previous section, Section 2; the other variables were considered while the impact on shear strength calculation was high. The direct relationship between each variable and the ultimate shear strength (*Vu*) of the RC deep beam is presented in Figure S1. Exponential, linear, logarithmic, and power functions were used to estimate the best direct relationship equation between *Vu* and input variables. Table 2 presents the summary of these functions.


**Table 2.** Direct relationship functions between Vu and input variables.

From Table 2, it can be observed that the power function has the best correlation with Vu in the case of using *a*/*d* and *fyv* of the main variables. The relationship of *ρ<sup>v</sup>* with *Vu* is exponential. The linear correlation can be detected between *Vu* and (*s* and *ρh* ). *ρ*, *fy*, and *f* ′ *<sup>c</sup>* are correlated with *Vu* based on logarithmic functions. The best R<sup>2</sup> between *Vu* and the main variables is 0.39 for the *f* ′ *<sup>c</sup>* variable. These results indicate that the relationship between the main variables and *Vu* cannot be estimated directly, and a complex relationship may be detected by using all main variables. Similarly, for the other variables, the relationship between *Vu* and variables varies. The best R<sup>2</sup> is estimated using a number of main bars, R <sup>2</sup> = 0.37. The variation in R<sup>2</sup> indicates the complex relationship between *Vu* and all variables. Table 2 and Figure S1 show the increase of the resistance with beam effective width and height, number of main bars, and concrete strength up to a certain limit.

The statistical evaluation, range (RA), average (M), standard deviation (SD), kurtosis (KU), and skewness (SK) of the used datasets is presented in Table 3. From the table, the range of datasets varies and will affect the models' performances, so the normalized datasets are used to overcome the range change of variables. The data is normalized between 0 and 1 in this study. In addition, in the proposed models, it is recommended to use the given ranges of input variables. The average and standard deviation values show that the distortion of datasets is high. The kurtosis and skewness values indicate that the distribution of datasets is not normal. Figure 2 presents the histogram and distribution of main variables and *Vu*. The figure shows positive skewness for whole variables is observed; the distribution is also supported by the presented values in Table 3.

**Figure 2.** *Cont*.

**Figure 2.** Distribution histogram of input and output variables.

**Table 3.** Statistical analysis of input and output variables.


#### **4. Methods and Development Models**

*4.1. Support Vector Regression*

Pal and Deswal [23] and Mozumder et al. [34] proposed the SVR formulas and theory in shear strength to predict RC beams. It was found to be a powerful computation technique for predicting the shear strength of deep beams [3,23]. Here, a summary of SVR is presented. SVR is the regression category of support vector machine (SVM), aiming to find a function that represents the relationship between inputs features to forecast the corresponding value when a new input is used. In SVR, "a fixed mapping procedure to map its input to n-dimensional feature space; then nonlinear functions are used to fit the high-dimensional features [35]". Vapnik [35] proposed a loss function to allow the concept of SVM margin

to be used for regression solutions. The following equation represents the mathematical formula for the SVRs' approximation function [23,34]:

$$f(\mathbf{x}) = w\varphi(\mathbf{x}) + b \tag{1}$$

$$\mathcal{C} = \frac{1}{2}w + \mathcal{C}\frac{1}{n}\sum\_{i=1}^{n}L(\mathbf{x},d) \tag{2}$$

where *w* and *ϕ* represent the weight vector and transformation functions, respectively; *x* and *d* are the input and output vectors, and *b* denotes a scalar. In Equation (1), the *w* and *b* are used to determine the normal and scalar vector, respectively, for the high-dimensional space, which is determined through *ϕ*(*x*). Terms <sup>1</sup> <sup>2</sup>*w* and *C* 1 *n n* ∑ *i*=1 *L*(*x*, *d*) in Equation (2) are the standard error and the penalty terms, respectively. The *ε*-insensitive loss function introduced by Vapnik [34] is commonly used to estimate the Equation (1) parameters through the following minimization function [23,34]:

$$\text{minimize } 0.5 \|\ w \parallel ^2 + \mathbb{C} \sum\_{i=1}^n (\mathfrak{f}\_i + \mathfrak{f}\_i^\*) \quad \text{subject to } \begin{cases} \displaystyle d\_i - (w\chi\_i) - b \le \varepsilon + \mathfrak{f}\_i \\ (w\chi\_i) + b - d\_i \le \varepsilon + \mathfrak{f}\_i^\* \\ \mathfrak{f}\_i \mathfrak{f}\_i^\* \end{cases} \tag{3}$$

where *C* > 0, which controls the trade-off between the model complexity and the amount up to which deviations larger than ε are tolerated. Equation (3) can be transformed to a dual space using Lagrange multipliers solution. This solution can be expressed as follows [23,34]:

$$\begin{aligned} \text{maximize } L &= -\varepsilon \sum\_{i=1}^{n} \left( \mathbf{a}\_{i}^{\*} + \mathbf{a}\_{i} \right) + \sum\_{i=1}^{n} d\_{i} \left( \mathbf{a}\_{i}^{\*} - \mathbf{a}\_{i} \right) \\ &- 0.5 \sum\_{i=1}^{n} \sum\_{j=1}^{n} \left( \mathbf{a}\_{i}^{\*} - \mathbf{a}\_{i} \right) \left( \mathbf{a}\_{j}^{\*} + \mathbf{a}\_{j} \right) \left( \mathbf{x}\_{i} - \mathbf{x}\_{j} \right) \quad \text{subject to } \begin{cases} \sum\_{i=1}^{n} \left( \mathbf{a}\_{i} - \mathbf{a}\_{i}^{\*} \right) = \mathbf{0} \\ \mathbf{0} \le \mathbf{a}\_{i}^{\*} \le \mathbf{C} \\ \mathbf{0} \le \mathbf{a}\_{i} \le \mathbf{C} \end{cases} \end{aligned} \tag{4}$$

where *L* denotes the Lagrangian and *α<sup>i</sup>* and *α* ∗ *i* represent the Lagrange multiplier. Once Equation (4) is used to estimate the parameters of Equation (2), Equation (1) can be rewritten as [23,34]:

$$f(\mathbf{x}) = \sum\_{\mathbf{n} \mathbf{s} \mathbf{v}} (\mathbf{a}\_i^\* - \mathbf{a}\_i)(\mathbf{x}\_k.\mathbf{x}) + b \tag{5}$$

where *nsv* represents the number of support vectors (*x<sup>r</sup>* , *xs*). Here, the solution of this equation depends on the training pattern of Lagrange multipliers, which are only applied to estimate the *w* and *b*. Therefore, the kernel function is commonly used to solve the nonlinear regression problems in SVR. The Kernel functions can transform the nonlinear problems into linear problems, as presented in Yaseen et al. [36], which allows the SVR to solve more complex problems. The nonlinear SVR can be expressed as follows:

$$f(\mathbf{x}) = \sum\_{n \ge v} (\mathbf{a}\_i^\* - \mathbf{a}\_i) \mathbf{K}(\mathbf{x}\_i \mathbf{x}) + b \tag{6}$$

where *K* is the kernel function; *K*(*xix*) = (*ϕ*(*x<sup>i</sup>* )*ϕ*(*x*)). In the current study, the radial basis kernel (RBF) is applied.

It should be mentioned that the SVR is built with statistical theory and based on the minimization of structural risk principle [37]. It is a popular method for a small count of data, high dimensional, and non-linear problems. Therefore, SVR is used in many applications [13,38,39] for prediction tasks. Generally, SVR is a type of convex optimization technique to search a local solution within a problem domain [37]. The tuning of learning parameters of SVR greatly impacts the evaluation quality. Therefore, finding optimal values of SVR parameters from global searched cost is a difficult task [40]. Nature-inspired algorithms are proved to be successful in finding the local best solution from the global one.

This work applied and showed PSO, HHO, and AVOA for optimizing SVR parameters with faster convergence capability to provide better prediction accuracy of the SVR model in the deep beam's shear strength prediction.

#### *4.2. Optimization Methods*

To optimize the best parameters (*C*, *ε*, and KernelScale (*γ*)) of SVR, the PSO, HHO, and AVOA algorithms are used separately in the current work.

#### 4.2.1. PSO

Kennedy and Eberhart (1995) proposed a metaheuristic population-based optimization algorithm named particle swarm optimization (PSO). The mechanism of the PSO is inspired by the fish schooling and foraging of the flock of birds while exploring an unknown region. Further, PSO is defined as a swarm of particles moving nearby the problem space by the influence of its global (*Gbest*) and local best (*Pbest*) position [41]. Therefore, a population search algorithm is used in PSO in a nature-inspired manner to analyze the input data features. The best position of the whole swarm is required to optimize the parameters. These can be performed through nature-inspired behaviors and learning experiences of population particles. PSO was found to be a robust integrated technique with SVR and ANN to model the shear strength of concrete [36,41,42]. PSO algorithm may be summarized as follows:


$$w\_{ik} = wv\_{ik} + \alpha e f\_1 rand\_1 (P\_{\text{best},ik} - y\_{ik}) + \alpha e f\_2 rand\_2 (G\_{\text{best},ik} - y\_{ik}) \tag{7}$$

$$y\_{ik} = y\_{ik} + v\_{ik} \tag{8}$$

where *y<sup>i</sup>* denotes the *i*-th particle, *k* = the *k*-th dimension of the particle, *coe f* <sup>1</sup> and *coe f* <sup>2</sup> represent the acceleration coefficients, *w* refers to the inertia weight, *rand* <sup>1</sup> and *rand* <sup>2</sup> represent the random coefficients, which are randomly limited between zero and one. More details for PSO theory can be found in [43,44].

#### 4.2.2. HHO

Heaidari and Mirjalili (2019) proposed a gradient-free, population-based optimization technique named Harris hawks optimization (HHO) [45]. The main inspiration behind HHO is surprised pounce, i.e., the chasing style and cooperative behavior of Harris' hawks in nature. According to the HHO optimization technique, numerous hawks cooperate to surprise a prey by pouncing it from multiple directions. They display a variety of pursuit patterns based on different scenarios and escaping patterns of the prey. The mechanism of HHO is that it mimics the Harris' hawk behavior in that trace, encircle, flush out, and attack the prey. It has been integrated with SVR and ANN to model the concrete characteristics and other engineering applications [21,46–49]. The main phases in the attacking of hawks are exploration (phase 1), transferring (phase 2), and exploitation (phase 3). In phase 1, the hawk depends on its position from the prey based on his waiting, seeking, and discovering; this can be expressed as follows:

$$\mathbf{Y}(iter+1) = \begin{cases} \mathbf{Y}\_{num}(iter) - r\_1|\mathbf{Y}\_{num}(iter) - 2r\_2\mathbf{Y}\_{num}(iter)|\text{ if } n \ge 0.5\\ \left(\mathbf{Y}\_{prev(iter) - \mathbf{Y}\_n(iter)}\right) - r\_3(LB + r\_4(IL - LL)) \text{ if } n < 0.5 \end{cases} \tag{9}$$

where Y*ranm* and Y*prey* indicate the random position for the selected hawk and prey's position, respectively. *UL* and *LL* indicate the upper and lower range; *r<sup>i</sup>* indicates a random number, having a value between 0 and 1; and Y*<sup>m</sup>* = 1/*N* ∑ *N* <sup>1</sup> Y*i*(*iter*).

In phase 2, the prey energy is modeled as *E* = 2*E*<sup>0</sup> <sup>1</sup> <sup>−</sup> *iter T* , where *T* and *E*<sup>0</sup> ∈ (−1, 1), and they indicate that energy falls for the prey with their escapes. Thus, the hawk can decide the solution based on the E computation and starting in phase 3 when |*E*| ≥ 1, and exploiting the neighborhood when <sup>|</sup>*E*<sup>|</sup> <sup>&</sup>lt; 1. Once starting phase 3, hawks decide to apply a soft or hard besiege. |*E*| ≥ 0.5 indicates the prey still has enough energy to escape, but maybe some misleading jumps occur in it to fail, so a soft besiege works. On the other hand, in the case of <sup>|</sup>*E*<sup>|</sup> <sup>&</sup>lt; 0.5, the prey is too fatigued to escape, so hard besiege works. Here, the HHO is used to optimize the SVR parameters.

#### 4.2.3. AVOA

Abdollahzadeha and Gharehchopogh (2021) recently introduced a metaheuristic algorithm named African vultures' optimization algorithm (AVOA) [17]. The inspiration behind the development of the AVOA algorithm is the competing and searching behavior of vultures to acquire a large amount of food. To acquire a large amount of food, these vultures, '*N*' (*N* denotes the population of vulture), were divided into categories based on their fitness to find food and eat. The solution with the highest fitness value is treated as the first-best vulture and the second-best solution as the second-best vulture. The rest of the vultures were trying to approach the best vulture. This is formulated as follows.

Step 1: To determine the best vulture in the group. Fitness of all solutions is determined, and the best solution is selected as the best vulture of the group and other solutions will move towards the best solution using:

$$R(i) = \begin{cases} \
multicolumn{3}{c}{f} \
ormalfont \text{ $p\_{best1}$ } \
if \ p\_i = \mathcal{K}\_1 \\ 
multicolumn{3}{c}{f} \
with \ p\_i = \mathcal{K}\_2 \end{cases} \tag{10}$$

where the value of *K*<sup>1</sup> and *K*<sup>2</sup> lies between 0 and 1 with their sum equal to 1.

Step 2: Starvation rate of vultures. The starvation rate is the rate at which the vultures are satiated or hungry. The satiated rate has a declining trend, and to model, behavior Equation (11) is used,

$$F = (2ram\_1 + 1)P\left(1 - \frac{iter\_i}{iter\_{\max}}\right) + t \tag{11}$$

$$t = h \left( \sin^w \left( 0.5 \pi \frac{iter\_i}{iter\_{\max}} \right) + \cos \left( 0.5 \pi \frac{iter\_i}{iter\_{\max}} \right) - 1 \right) \tag{12}$$

where *F* indicates the vultures are satiated. If the value of |*F*| is greater than 1, vultures search for food in different areas, and the algorithm enters the exploration phase, whereas if the value is less than 1, AVOA enters the exploitation phase, and vultures search for food in the neighborhood. *iter<sup>i</sup>* Indicates the iteration number, *itermax* indicates the total number of iterations, *ranm*<sup>1</sup> has a random value between 0 and 1. Here, *z* and *h* indicate random numbers with values lying between −1 to 1 and −2 to 2, respectively. If the value of *z* is less than 0, it indicates the vulture is starved, and if it increases to 0, it indicates the vulture is satiated. Here, w indicates the optimization operation disrupts the exploration and operation phases. By increasing the value of *w*, the probability of entering the exploration phase in the final optimization stages increases, and vice versa for decreasing the value of *w*.

Step 3: Exploration phase. In this phase, different random areas can be examined using two different strategies. To select the strategies in the *ranm*<sup>1</sup> exploration phase, a random number between 0 and 1 is generated. This procedure is shown in Equation (13).

$$P(i+1) = \begin{cases} \begin{array}{c} \mathcal{R}(i) - |XR(i) - P(i)| \, \text{if } P\_1 \ge \text{ramm}\_{P1} \\ \mathcal{R}(i) - \mathcal{F} + \text{ramm}\_2((\text{LB} - \text{LB}) \, \text{ramm}\_3 + \text{LB}) \, \text{if } P\_1 < \text{ramm}\_{P1} \end{array} \end{cases} \tag{13}$$

where *P*(*i* + 1) indicates the vulture position vector in the next iteration, *F* indicates the rate of vulture being satiated, and *R*(*i*) indicates one of the best vultures, which is selected at step 1 in the current iteration. *X* indicates a coefficient vector that increases the random motion by changing in each iteration and is obtained using the formula *X* = 2 × *ranm*, where *ranm* is a random number between 0 and 1.

Step 4: Exploitation phase. In this phase, the efficiency stage of the algorithm is investigated. If the value of |*F*| is less than 1, the algorithm enters the exploitation phase. Here, the exploitation phase is categorized based on the |*F*| value. If |*F*| is between 1 and 0.5, the rotating flight strategy of the vulture is processed based on parameter *P*2. This can be processed as follows:

$$P(i+1) = \begin{cases} |XR(i) - P(i)|(F + num\_4) - (R(i) - P(i)) \text{ if } P\_2 \ge num\_{P2} \\\ R(i) - \left[ \left( R(i) (\cos(P(i))) \left( \frac{num\_4 P(i)}{2\pi} \right) \right) + \left( R(i) (\sin(P(i))) \left( \frac{num\_6 P(i)}{2\pi} \right) \right) \right] \text{ if } P\_1 < num\_{P2} \end{cases} \tag{14}$$

If |*F*| is less than 0.5, the two vultures' movements accumulate several types of vultures over the food sources, and the siege and aggressive strife to find food are implemented using parameter *P*3. This can be defined as follows:

$$P(i+1) = \begin{cases} \frac{A\_1 + A\_2}{2} \text{ if } P\_3 \ge rann{\,}m \,\updownarrow\,\, \text{m} \\\ \,\, R(i) - (|R(i) - P(i)|) (F) (LF(d)) \,\, \text{if } P\_3 < rann{\,}\,\text{m} \end{cases} \tag{15}$$

where,

$$A\_1 = 
overline{\text{best}\_{\text{best1}}(i)} - \text{F}\frac{\text{vulltrue}\_{\text{best1}}(i)P(i)}{\text{vulltrue}\_{\text{best1}}(i) - P(i)^2}; \text{and } A\_2 = 
overline{\text{best2}}(i) - \text{F}\frac{\text{vulltrue}\_{\text{best2}}(i)P(i)}{\text{vulltrue}\_{\text{best2}}(i) - P(i)^2} \tag{16}$$

$$LF(x) = 0.01 \frac{\mu \sigma}{|v|^{1/\beta}}, \sigma = \left(\frac{\Gamma(1+\beta)\sin(\pi\beta/2)}{\Gamma(1+\beta^2)\beta\mathbb{2}((\beta-1)/2)}\right)^{1/\beta} \tag{17}$$

In which, *vulturebest*1(*i*) and *vulturebest*2(*i*) are the best vulture of the first and second groups, respectively, in the iteration *i*; *d* is the problem dimensions, *µ* and *v* denotes a random number between 0 and 1, and *β* = 1.5.

In this work, AVOA is used to tune the SVR parameter set to find the efficient performance of the SVR model.

#### *4.3. Models' Development and Accuracy Assessment*

This study proposed a new hybrid AVOA metaheuristic algorithm-based SVR model to find the optimal parameters (*C*, *ε*, and *γ*) of SVR and compare its performance with other two metaheuristics nature-inspired algorithms (called PSO and HHO)-based SVR model. Generally, the robustness of the SVR model depends on an appropriate selection of the parameters named as the penalty parameter/"BoxConstraint" (*C*), insensitive loss function/"epsilon" (*ε*), and the kernel parameter/"KernelScale" (*γ*/gamma). The range of these parameters is large, and it is difficult to search for the optimal set of values for these three parameters. Therefore, this optimization task may be solved using optimization algorithms. The authors of this article used three metaheuristic algorithms (AVOA, PSO, and HHO) to find the optimal parameter set of the SVR model. Figure 3 shows the process flow of the proposed technique. Initially, the missing value from the dataset is replaced using the k-nearest neighbor (KNN) imputer method. This study used a Euclidean distance measure to fill the missing value. The whole dataset is bifurcated into a train (80%) and test (20%) set. The three-metaheuristic algorithm is used to train SVR parameters for all/selected featured datasets separately. Root means squared error (RMSE) is selected as the fitness function of all algorithms. Metaheuristic algorithms are sensitive to their different parameter set. The Hit and trail approach is used to select the initial parameter set of metaheuristic algorithms. Table 4 shows the initial value of all parameters of three algorithms for SVR training. Since the number of epochs and population size affect the

convergence rate of metaheuristic algorithms, this research work aims to find a faster convergence rate of metaheuristic algorithm with 15 epochs and five population sizes (Table 4).

**Figure 3.** Process flowchart of the proposed method.

**Table 4.** Initial parameters of metaheuristic algorithms to train the SVR model.


To evaluate the accuracy of the proposed models, eight statistical indices are used: coefficient of determination (*R* 2 ), variance account factor (*VAF*), variance inflation factor (*VIF*), mean absolute error (*MAE*), root mean square error (*RMSE*), performance index (*PI*), mean bias error (*MBE*), and percentage error (*PE*). The *R* <sup>2</sup> and *VAF* are used to measure the correlation between the measured and predicted values. *VIF* is used to evaluate the collinearity between the measured and predicted values; *VIF* > 10 indicates high collinearity. The models' errors are evaluated using *MAE*, *RMSE*, and *MBE*, and the *PE* is used to estimate the accuracy of the proposed model error in predicting the shear strength of *RC* deep beams. The mathematical expression of these indices can be expressed as follows:

$$\mathcal{R}^2 = \frac{\sum\_{i=1}^{N} (V\_i - V\_{mean})^2 - \sum\_{i=1}^{N} \left(V\_i - V\_{pi}\right)^2}{\sum\_{i=1}^{N} \left(V\_i - V\_{mean}\right)^2} \tag{18}$$

$$VAF = 100 \left( 1 - \frac{var\left(V\_i - V\_{pi}\right)}{var\left(V\_i\right)} \right) \tag{19}$$

$$VIF = \frac{1}{1 - R^2} \tag{20}$$

$$RMSE = \sqrt{\frac{\sum\_{i=1}^{N} \left(V\_i - V\_{pi}\right)^2}{N}} \tag{21}$$

$$MAE = \frac{\sum\_{i=1}^{N} |\left(V\_i - V\_{pi}\right)|}{N} \tag{22}$$

$$MBE = \frac{1}{N} \sum\_{i=1}^{N} \left(V\_i - V\_{pi}\right) \tag{23}$$

$$PI = adj \ R^2 + (0.01 VAF) - RMSE \tag{24}$$

$$PE = 100 \frac{RMSE}{V\_{\max} - V\_{\min}} \tag{25}$$

where *V<sup>i</sup>* and *Vpi* represent the measured and predicted shear strength, *Vmean*, *Vmax*, and *Vmin* are the average, maximum, and minimum, respectively, of measured values, *adj R*<sup>2</sup> is the adjustment *R* 2 , and *N* is the number of the data sample.

#### *4.4. Sensitivity Analysis*

Cosine Amplitude Method (CAM) is used to analyze the strength of the relation between input the parameter and power factor [50]. It can also be used to determine the express similarity relation between correlated parameters. To apply CAM, all the data pairs were stated in common *X*-space. The data pairs used to construct a data array defined *K* as:

$$\mathbf{K} = \{ \mathbf{K}\_1, \mathbf{K}\_2, \mathbf{K}\_3, \dots, \mathbf{K}\_n \} \tag{26}$$

Every elements i.e., *K<sup>i</sup>* , in the data array *K* is a vector of lengths *j*, i.e.,:

$$K\_{i} = \left\{ k\_{i1}, k\_{i2}, k\_{i3}, \dots, \mathbb{K}\_{ij} \right\}\_{1} \tag{27}$$

Therefore, each of the data pairs is represented as a point in *m*-dimensional space, where each point requires *j*-coordinates for its complete description.

#### **5. Results and Discussion**

Two scenarios are presented in this section. The first is the evaluation of the proposed models based on all variables and the study of the effect of all variables on *Vu* estimation. Second, the main variables are considered in modeling and evaluating the sensitivity of

input variables on the best selection model. The best solution is compared to existing models in the shear strength determination of RC deep beams.

#### *5.1. All Variables Impact on Vu Estimation*

Table 5 presents the statistical indices of the proposed models. The numerical investigation of the statistical indices shows that the performance of the AVOA-SVR model is better to estimate *Vu* with R<sup>2</sup> = 0.98 and RMSE = 32.20 kN in the training stage. The comparison between all models using performance indices of all statistical indices in the training and testing stages shows that the AVOA-SVR model has a high index. However, the performance of HHO-SVR is better in terms of PI, RMSE, and PE in the testing stage. Moreover, Figure 4 illustrates the linear correlation between experimental and predicted *Vu* values for the proposed models. From Figure 4, it is shown that the performance of the AVOA-SVR model is more accurate than other proposed models in the training and testing stages. The distortion of data points around best fitting is small in modeling *Vu* with the AVOA-SVR model, and the VIF is higher than in other models, as presented in Table 5. This means that when we used all variables, the AVOA-SVR model can be used to estimate *Vu* with a model error approach of 6.95%.


**Table 5.** Statistical evaluation of the proposed models.

The comparison between the AVOA-SVR model and previous studies is presented in Table 6. The statistical evaluation of relative predicted shear strength (*Vu* measured/predicted (*Vm*/*p*)) is presented in Table 6; COV is the coefficient of variation. From this table, it can be observed that the AVOA-SVR model performance is better and slightly better than ACI and Russo algorithms, respectively. Although the distortion around the mean for the proposed model is lower than for the Russo technique, the range of datasets for Russo is better than the proposed models. The Liu technique is better than previous and proposed models when considering the whole variables in modeling shear strength through the AVOA-SVR model. The selected variables are used and evaluated in the next section to estimate more accurate *Vu* values.

**Table 6.** *Vm*/*p* statistical evaluation for AVOA-SVR and previous studies.


**Figure 4.** Scatter plot of model's performances in the training (**upper row**) and testing (**lower row**) stages.

#### *5.2. Selected Variables Impact on Vu Estimation*

Table 7 and Figure 5 present the performance evaluation of the proposed models. The performance of AVOA-SVR is high in the training stage. A high correlation, R<sup>2</sup> = 0.97, and low model error, PE = 2.25%, are observed. In the testing stage, a low distortion around best fitting is observed with the AVOA-SVR. In addition, the statistical correlation factors are high, R<sup>2</sup> = 0.97 and VAF = 94.46, as presented in Table 7 and Figure 5. The VIF values of AVOA-SVR in the training and testing stages are higher than other models. This means the accuracy of AVOA-SVR is acceptable with low distortion around the observed values. Although the PI and RMSE of the HHO-SVR models are lower than for the AVOA-SVR model, the distortion of HHO-SVR datasets is high, as presented in Table 7 and Figure 5. Meanwhile, the performance of the proposed models is shown to be low to estimate the high shear strength (as presented in Figures 4 and 5). However, the performance of AVOA-SVR is seen as more robust. This indicates that AVOA-SVR can overcome the variation change in the data used. Figure 6 also shows a faster convergence rate of the AVOA-SVR model compared to the other two hybrid models. Therefore, the AVOA-SVR can be used to estimate the shear strength of RC deep beams with 3.4% model accuracy. The statistical comparison indices in Tables 5 and 7 show that the performance of proposed models with the selected variables is better than that for using all variables in modeling the proposed techniques. This means that the selected variables are more influential in the shear strength of RC deep beams.

**Figure 5.** Scatter plot of model's performances in the training (**upper row**) and testing (**lower row**) stages.

**Table 7.** Statistical evaluation of the proposed models.


**Figure 6.** Convergence rate of three models.

#### *5.3. Comparison with Previous Studies and Codes*

Table 8 and Figure 7 show the performance of the proposed model compared to the previous studies' formulas. As seen in Table 8, the used selected variable improved the AVOA-SVR performance by 60% in terms of COV. This indicates that the selected variables are significantly affected by the *Vu* values of the RC deep beams. The comparison between previous formulas and AVOA-SVR shows that the proposed model accuracy is high in estimating the shear strength of RC deep beams. The small range is estimated with AVOA-SVR, and the range is 0.57 kN. The small SD and COV of the statistical indices are observed with AVOA-SVR. This means the accuracy of AVOA-SVR is high compared to other models.


**Table 8.** *Vm*/*p* statistical evaluation for AVOA-SVR and previous studies.

The boxplot in Figure 7a shows that the median, red horizontal line of Russo, is close to the true value "1", followed by the AVOA-SVR and Liu models, respectively. The low interquartile range (IQR) value, the height of the box, is observed to be small with the AVOA-SVR model, followed by Liu and Russo models, respectively. The maximum and minimum quartiles, the black horizontal solid lines, are small with the AVOA-SVR model, followed by Liu and Russo models, respectively. The outliers are observed near the median of the AVOA-SVR and far to the median of the ACI model. From the visualization of boxplot results, it can be concluded that the performance of the AVOA-SVR model is more accurate than the previous studies for modeling the shear strength of RC deep beams. In addition, the following model is the Liu model, as this model considers more shear resistance mechanisms and shows a higher normal distribution and lower error than Russo's and ACI's models. The quantile-quantile (Q-Q) plot is presented in Figure 7b for the Liu and AVOA-SVR models for further investigation. The relative shear strength is presented versus the standard normal distribution. From this figure, both models have approximately followed the normal distribution; this indicates that both models can be used to estimate the shear strength of the RC deep beam. The AVOA-SVR model has more correlation with the standard normal distribution, and the VAOA-SVR model is more accurate than the Liu technique in modeling the shear strength of RC deep beams. The scatter plot presented in Figure 7c,d shows that the worst model for estimating the shear strength is the ACI's model. The variation in the best solution "1" is shown as small for Russo, Liu, and AVOA-SVR models, respectively. The most of relative shear strength of the AVOA-SVR model falls within ±20%. The comparison of the AVOA-SVR model and previous models shows that the developed model can be used accurately to model the shear strength of the RC deep beams. Therefore, the AVOA-SVR model is a potential soft computing technique that can be used in predicting the shear strength of RC deep beams.

**Figure 7.** Comparison of model's performances, (**a**) boxplot, (**b**) Q-Q plot, (**c**) scatter plot of relative shear strength with measured shear strength, and (**d**) zoom in for upper plot with +20% limits for the best models.

#### *5.4. Sensitivity Analysis of Input Variables*

Figure 8 presents the most influential input variables in modeling the shear strength of the RC deep beam. The impact of the input variables is presented for the three models. From Figure 8, it can be noticed that the significant impact of the ratio of vertical and horizontal web reinforcements is low compared to other variables. The sensitivity of the shear span to depth ratio is high, followed by the yield strength of the main steel, the ratio of the main tensile bars, yield strength of vertical web reinforcement, stirrups spacing, and concrete compressive strength, respectively. The impact of the input variables on output for the other developed models is similar. These results imply that the shear strength of the RC deep beam is highly influenced by the beam geometry, concrete strength, and yield

strength of the steel bars. The stirrups spacing also has a large effect on the shear strength of RC deep beams.

**Figure 8.** Sensitivity of input variable on model's prediction.

#### **6. Conclusions**

This study investigated the use of new metaheuristic optimization algorithms integrated with SVR to model the shear strength of RC deep beams and evaluate the sensitivity of input parameters. SVR-AVOA, -PSO, and -HHO were designed and compared to existing models in the current study. In this study, 202 datasets, including 19 variables of experimental studies, were collected from literature to design and evaluate the proposed models. The common eight parameters (shear span to depth ratio, the ratio of the main tensile bars, yield strength of main bars, concrete compressive strength, the ratio of vertical web reinforcement (stirrups), stirrups spacing, yield strength of vertical web reinforcement, and the ratio of horizontal web reinforcement) are also used to evaluate the performance of the proposed models' in predicting the shear strength of RC deep beams. The performance of SVR-AVOA is high in the cases of the used 19 and 8 parameters for modeling the shear strength. The accuracy of SVR-AVOA is improved by 60%, in COV terms, using the common input variables. Thus, other parameters were found less significant in modeling the shear strength of RC deep beams. The comparison of the SVR-AVOA and the previous studies shows that the accuracy of the proposed model is higher than Liu [4], Russo [8], and ACI [10] by 46%, 63%, and 90%, respectively, in terms of COV. This indicates that SVR-AVOA is the more robust model and can be accurately used in modeling the shear strength of RC deep beams. The sensitivity of the input variables in modeling the shear strength of RC beams with the SVR-AVOA was assessed. This investigation shows the impact of the shear span on the beam's depth ratio, yield strengths of vertical and horizontal web reinforcement, concrete compressive strength, stirrups spacing, and the ratio of the main longitudinal bars on the deep beams' shear strength.

Furthermore, the sensitivity of AVOA algorithm parameters can be tested to balance between exploitation and exploration side for enhancing the SVR performance. In the future, to check the efficiency of the proposed model should be tested on other datasets and other civil engineering application areas. The AVOA algorithm can be combined with other machine learning models like an extreme learning machine, random forest, etc., for prediction tasks.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/su14095238/s1, Figure S1: Direct relationship between inputs and output; Table S1: Data used in Modeling; Table S2: Modeling variables.

**Author Contributions:** Conceptualization, M.R.K. and B.S.A.; methodology, M.R.K. and B.R.; software, B.R.; validation, M.R.K., B.R. and B.S.A.; formal analysis, M.R.K. and B.S.A.; investigation, M.R.K., B.R. and B.S.A.; resources, B.S.A.; data curation, M.R.K. and B.S.A.; writing—original draft preparation, M.R.K., B.R. and B.S.A.; writing—review and editing, M.R.K., B.R., K.C. and B.S.A.; visualization, M.R.K., B.R., K.C., S.-M.K., H.-M.J. and B.S.A.; supervision, M.R.K. and J.-W.H.; project administration, M.R.K., J.-W.H. and B.S.A.; funding acquisition, J.-W.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by Incheon National University Research Concentration Professors Grant in 2021.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used are available in the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Deformation Capacity of RC Beam-Column Joints Strengthened with Ferrocement**

**M. Zardan Araby 1,2 , Samsul Rizal 3 , Abdullah 2 , Mochammad Afifuddin <sup>2</sup> and Muttaqin Hasan 2, \***


**Abstract:** Beam-column joints constructed in the pre-seismic building code do not provide transverse reinforcement and good reinforcement detailing within the region. These cause the occurrence of brittle shear failure, which is one of the factors affecting the number of reinforced concrete (RC) moment resistance building structures collapsing during an earthquake. Therefore, in this study a brittle beam-column joint with a non-seismic building code was designed and strengthened by a ferrocement. Four layers of wire mesh with a diameter of 1 mm and a mesh size of 25.4 mm were installed on both sides of the beam-column joint and cement mortar was cast on it. As a comparison, a ductile beam-column joint was also designed following the current building code, which considers seismic effects. The test results by applying reversed cyclic loading at the beam tip showed that strengthening using ferrocement prevents crack propagation, increasing the deformation capacity, ductility, stiffness, and energy dissipation of beam column joint which are higher than those of the beam-column joint which is designed following the current building code. However, the strengthening does not improve the load carrying capacity significantly.

**Keywords:** beam-column joint; ferrocement; crack; ductility; displacement

#### **1. Introduction**

Reinforced concrete (RC) buildings constructed in the 1970s and 1980s or earlier lacked transverse reinforcement installations in the beam-column joint region. Furthermore, they did not adhere to the requirements of a detailed reinforcement to withstand seismic loads because the building code at that time did not accommodate these effects on the beamcolumn joint structures. In Indonesia, the building code used to design certain structures during that period was NI-2 [1]. The absence of special provisions regarding reinforcement detailing and transversal reinforcement in beam-column joint region due to the ductile design philosophy has not been adopted yet, led to the failure of the beam-column joints. Globally, it is one of the causes of collapsed buildings during earthquakes [2–14]. During that era, many buildings were also designed without the strong column–weak beam design philosophy which results in the appearance of column hinges in their collapse pattern. However, the column hinges were also found although the structures are designed based on the strong column–weak beam criterion [15].

Generally, the beam-column joint experiences a brittle shear failure. Therefore, efforts are needed to strengthen the structures that were built in the 70s and 80s for their sustainability [16–19]. Several structural strengthening methods have been proposed to withstand seismic loads. These include the use of Fiber Reinforced Polymer (FRP) [20–40], steel jacketing [20,41–44], pre-fabricated composite blocks [45], diagonal steel bars [46], steel haunch strengthening and confining joint reinforcement [44,47], injection of cracks with epoxy and

**Citation:** Araby, M.Z.; Rizal, S.; Abdullah; Afifuddin, M.; Hasan, M. Deformation Capacity of RC Beam-Column Joints Strengthened with Ferrocement. *Sustainability* **2022**, *14*, 4398. https://doi.org/10.3390/ su14084398

Academic Editors: Maged A. Youssef, Chiara Bedon, Mislav Stepinac, Marco Fasan and Ajitanshu Vedrtnam

Received: 10 March 2022 Accepted: 4 April 2022 Published: 7 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

concrete jacket [20,48,49], textile-reinforced engineered cementitious composites [50,51], installing reinforced concrete wing walls [52], and modified reinforcement technique [53]. However, those strengthening systems are expensive and need to be professionally installed.

Ferrocement is a type of thin reinforced concrete usually made of hydraulic cement mortar reinforced with a metallic mesh or similar materials [54] and has been used as a material for the strengthening of reinforced concrete structures. It possesses greater tensile strength and resistance to cracking than conventional reinforced concrete. Despite having similar durability, ferrocement is extremely elastic. In this study, ferrocement was used to strengthen the beam-column joint, because this method is cheap compared to the above-described methods, and its materials are always available, simple, and easy to install. Several studies have been carried out on its use as a beam strengthening material, and it was discovered to increase flexural capacity, stiffness, ductility, and energy dissipation [55–61]. Similarly, as a column strengthening material, it also significantly increases ductility and energy dissipation as well as load-bearing capacity and stiffness [62–65]. However, very limited studies have been carried out on beam-column joint strengthening using ferrocement [66–68].

This study was carried out with the aim to understand the deformation capacity of a beam-column joint designed with a non-seismic building code strengthened with ferrocement under reversed cyclic loading. The reversed cyclic loading was chosen to simulate the earthquake action in building structures since the strengthening method proposed in this study will be applied in building structures in a seismic prone area. The strengthened beam-column joint performance is then compared to the unreinforced one and used to ascertain the increase in its deformation capacity. Furthermore, it is also compared with the performance of beam-column joint designed with the current Indonesian building code [69], which considers the seismic loads effect to discern the strengthening efficiency to withstand such impacts. Therefore, it is expected that this study recommends an economical and practical structural strengthening method that is applicable when strengthening existing buildings. It is important to note that the behaviors of strengthened beam-column joint presented in this paper are only based on the experimental results. The numerical analysis with reliable models of such structure is considered for further study.

#### **2. Experimental Program**

#### *2.1. Detail of Specimens*

Three beam-column joint specimens were prepared as shown in Table 1. The beams and columns are designed to have similar cross-sectional width and longitudinal reinforcements. The beam cross-section is 300 mm × 400 mm with 8D14 mm longitudinal reinforcements. Furthermore, the transverse reinforcements is in the form of D10 mm stirrups, which are installed every 100 mm distance from the column face to as far as 200 mm in its front. In the other part of the beam, D10 mm stirrups are installed within a distance of 50 mm. The column cross-section is 300 mm × 300 mm with longitudinal reinforcements of 8D14 mm. Conversely, the transverse reinforcements in the form of D10 mm stirrups are installed every 100 mm distance from the beam face to as far as 200 mm in its front. In the other part of the column, D10 mm stirrups are installed within a distance of 50 mm. The installation of a tighter stirrup in those areas is to ensure that cracks do not occur in this area when a load is applied. This is due to the fact the beam-column joint is the only part observed in this study. The weak column–strong beam was selected to represent buildings constructed in the 70s and 80s. The mix proportion and the concrete compressive strength as well as reinforcing bar yield strength and Young's modulus used for the beam-column joints are shown in Tables 2 and 3. The concrete compressive strength presented in Table 2 is the average compressive strength tested on 20 cylinder specimens with diameter of 150 mm and height of 300 mm at the age of 28 days, while the yield strength and elastic modulus of steel bars presented in Table 3 were obtained from the test results on one specimen for each bar diameter.



**Table 2.** Mix proportion of concrete.


**Table 3.** Yield strength and Young's modulus of reinforcing bars.


BCJ1 is a beam-column joint specimen designed according to NI-2 [1] without using transverse reinforcement in the joint region. The beam longitudinal reinforcing bars are not continuously bent towards the upper and lower columns. Figure 1 shows details of specimen BCJ1. Specimen BCJ2 has a similar size and longitudinal reinforcement as BCJI. The difference is that specimen BCJ2 added stirrups in the joint region using bars with a diameter of 10 mm placed at a spacing of 100 mm according to the current Indonesian building code [69]. The beam longitudinal reinforcing bars are also bent towards the upper and lower column with a length of 500 mm therefore their functions as anchors. Comprehensive details of specimen BCJ2 are shown in Figure 2. The specimen BCJ3 was designed in a similar manner as BCJ1. However, it was strengthened using ferrocement on both sides of the beam-column joint, as shown in Figure 3. Subsequently, the strengthening was provided in the following way. The first step involves disassembling the concrete cover of the specimen in the joint area, up to 400 mm in front of the beam and column, with a thickness until the reinforcing bars of the specimen is visible with the thickness of 40 mm. Furthermore, a T-shaped wire mesh with 1 mm diameter and 25.4 mm mesh size, similar to the beam-column joint area was cut. Four layers of wire mesh were installed and mounted on both sides of the specimen with the orientation angle of 0◦ to beam and column longitudinal axis as shown in Figure 4. Afterward, the wire mesh was anchored to the specimen using 4 dynabolts. Finally, a cement mortar was re-cast (Figure 5). Cement and sand with a volume ratio of 1:4 and water to cement ratio of 0.5 was used for the mortar mixture. A polycarboxylate-based superplasticizer having a specific gravity of 1.06 with a content of 1% of cement weight was also added to the mixture thereby enabling it flows easily when filling the cavities of the wire mesh during casting. Furthermore, the strengthening of the beam-column joint with ferrocement was carried out on both sides of the specimen by strengthening it on one side first, and after the mortar has hardened, the same procedure was performed on the other side.

**Figure 1.** Detail specimen BCJ1.

**Figure 2.** Detail specimen BCJ2.

**Figure 3.** Detail specimen BCJ3.

**Figure 4.** Wire mesh already installed on one side of the beam-column joint.

**Figure 5.** Re-casting the beam-column joint after installing the wire mesh.

## *2.2. Loading Procedure*

The specimen was first set on the loading frame, as shown in Figure 6. The column was placed horizontally above the loading frame by anchoring it at a distance of 200 mm from both ends using bolts. Furthermore, at both ends of the column, L shape steel is installed, welded to its longitudinal reinforcement, and anchored to the loading frame using bolts. The load is applied through an actuator driven by a hydraulic jack placed on the beam tip. A steel plate is installed on the beam surface to fasten the actuator and specimen, thereby providing a reversed cyclic loading under deformation control. Unloading and reloading were performed at the displacement of 0.75 mm, 1.5 mm, 3 mm, 6 mm, 12 mm, and 24 mm while two loading cycles were applied for each displacement as shown in Figure 7. After 12 cycles, the loading was continued with monotonic until the specimen failed.

**Figure 6.** Loading set up.

**Figure 7.** Displacement controlled loading history.

tance of 600 mm from the column face, as shown in Figure 6. In the joint area, two π gages Beam displacement was measured by installing 2 transducers on both sides at a distance of 600 mm from the column face, as shown in Figure 6. In the joint area, two π gages were mounted to measure the width of the crack that occurs. Furthermore, strain gages were installed at the longitudinal reinforcement as well as on the stirrups of the beam and column to measure the strains. Moreover, the applied load was measured with a load cell. All information acquired during loading was recorded with a data logger and entered into a computer. In addition, the crack pattern during the loading process was also observed and drawn until the failure of the specimen.

#### **3. Results and Discussion**

#### *3.1. Crack Propagation and Failure Mode*

Figure 8 clearly shows an image of the cracks that occurs at the BCJ1 specimen when the displacement was relatively 3 mm, 6 mm, 12 mm, and 24 mm. The crack did not occur until the displacement was approximately 1.5 mm. The first crack occurred at the beamcolumn joint corner in the longitudinal direction of the column towards the center of the beam in the direction of the push load. Furthermore, when the displacement was relatively 2.3 mm, the length of the first crack was increased to 300 mm. Another crack appeared at the same corner, however it was directed towards the center of the beam-column joint (inclined crack). When the displacement was approximately 3 mm, the crack that appeared first has turned towards the center of the column with a crack length of 100 mm and the inclined crack has reached the center of the beam-column joint as shown in Figure 8a. Figure 8b illustrates that when the displacement has reached 6 mm, 3 more cracks appeared, and with one starting at the corner of the beam-column joint in the longitudinal direction of the beam, it reached the other side of the column. Furthermore, another one started from the opposite side of the column in a slightly inclined direction, and the final crack appeared in the beam on the opposite side with the load applied and propagated to its center. Meanwhile, the first and second appeared cracks only increased in width. Consequently, when the displacement is relatively 12 mm, the second crack, which was an inclined crack, was propagated to the other side of the column until relatively 150 mm in front of the beam, as shown in Figure 8c. This inclined crack, as well as the horizontal one at the column face, increased in width with increasing in displacement. In addition, when the displacement was approximately 24 mm, another inclined crack appeared at another beam-column joint corner and propagated to the reverse side of the column until it was relatively 150 mm measured from the face of the beam. This crack formed an X shape together with the previously occurred inclined crack as shown in Figure 8d. Furthermore, lack of stirrups in the joint region led to a rapid increase in the width of the X crack and horizontal cracks that first appeared, thereby

causing a decrease in load when the displacement reached 25.27 mm and structural failure in the subsequent cycle at 32.58 mm. The failure mode of the specimen BCJ1 is shown in Figure 9. This is a typical shear failure in the beam-column joints.

**Figure 8.** Crack propagation of specimen BCJ1 at displacement of: (**a**) 3 mm; (**b**) 6 mm; (**c**) 12 mm; and (**d**) 24 mm.

‐

**Figure 9.** The failure mode of specimen BCJ1.

‐

Crack propagation of specimen BCJ2 at a displacement of 3 mm, 6 mm, 12 mm, and 24 mm is shown in Figure 10. A crack occurred when a displacement exceeded 1.5 mm at the beam-column joint corner on the opposite side of the applied loading side. The crack propagated parallel to the longitudinal direction of the column. However, when the displacement was approximately 3 mm, the crack length was relatively 190 mm, as shown in Figure 10a. Another crack with a length of 30 mm simultaneously appeared on the side of the beam opposite the load. The application of a compressive load in the same cycle led to the formation of new cracks at an opposite corner which was also horizontal towards the previous one, thereby causing the two cracks to meet. At the time of applying the tensile load for the next cycle, a new crack with a length of 220 mm occurred, starting from the first one, which was located 50 mm from the face of the beam towards column's longitudinal axis. Meanwhile, when the compressive load was applied in the same cycle, a similar incident where another crack occurred in the longitudinal axis direction of the beam leading to the inner joint, which started at the previous one with a length of 190 mm, as shown in Figure 10b. Besides, when the displacement was greater than 6 mm, inclined cracks started to occur at the joint. Furthermore, there were also cracks on the other side of the column and jointed with previous cracks in the column axis at the joint region when the displacement is relatively 12 mm, as shown in Figure 10c. In addition, when the crack in the beam propagated, its length was relatively 180 mm apart from the appearance of the other two cracks with similar lengths. Subsequently, when the displacement was approximately 24 mm, the cracks at the joint formed the X shape, as shown in Figure 10d. However, those X shape cracks did not increase in width when the displacement was greater than 24 mm, due to the stirrups provided in the joint region. This was different from specimen BCJ1, where the crack got wider and caused structural failure. In the BCJ2 case, the structure was still capable of deforming up to 52.39 mm and had only failed recently. Structural failure was also not caused by X-shaped cracks that occurred in the joints. However, it was caused by the widening of the cracks at the beam-column joint corner. The failure mode of specimen BCJ2 is shown in Figure 11.

Crack propagation of the specimen BCJ3 at displacement of relatively 3 mm, 6 mm, 12 mm, and 24 mm are shown in Figure 12. The first crack started immediately after the displacement exceeded 1.5 mm. Cracks also started from the beam-column joint corner on the loading side. In contrast to the BCJ1 and BCJ2 specimens, which the first crack was initially directed horizontally, in the BCJ3 specimen, the first crack had an inclined shape directed towards the center of the beam-column joint. Changes in the load direction led to the formation of a second crack on the other beam-column joint corner, which also had an inclined shape towards the beam-column joint center. In the next cycle a displacement of 3 mm was reached, and the 2 inclined cracks met at the middle of the column, in a vertical direction, away from the center of the beam-column joints. There was a cracked branch at a depth of approximately 1/3 of the column height in a horizontal direction towards its center as shown in Figure 12a. This crack failed to propagate and increase its width until the displacement was relatively 6 mm due to the presence of wire mesh. This was installed as a strengthening to prevent the propagation and widening of the crack as shown in Figure 12b. However, when the displacement was approximately 6 mm, another fine crack with a length of 5 mm appeared on the side of the column. As the displacement increased, the newly emerged inclined crack propagated towards the center of the beam-column joints. The initial crack also propagated towards the other side of the column. Furthermore, 2 other inclined cracks were formed, starting with the horizontal one that appeared first. Additionally, vertical and horizontal cracks appeared on the side of the column and beam respectively. There were also 2 vertical cracks in the beam ferrocement section. The shape of the crack when the displacement was relatively 12 mm is shown in Figure 12c. The crack pattern remained the same as shown in Figure 12d till the displacement was relatively 24 mm. There was an increase in the number of cracks on the beam and an additional vertical one on the side of the column. In the beam-column joint area, there was no increase in the length and width of the cracks. It was due to the presence of wire mesh installed

on the beam-column joint as a strengthening way to prevent crack growth. This condition occurred in the next load cycle until the connection between the ferrocement and the old concrete surface was damaged at a displacement of 51.37 mm. This is almost similar to the maximum displacement of the BCJ2 specimen. The failure mode of the beam-column joints strengthened with ferrocement was delamination of the ferrocement from the old concrete as shown in Figure 13. Supposing the bond between the old concrete and ferrocement is made stronger, it is believed that the displacement achieved by this beam-column joint may increase, because at the time of failure the crack width in the joint area was less than 1 mm, thereby improves its deformation capacity and ductility.

the first one, which was located 50 mm from the face of the beam towards column's lon-

**Figure 10.** Crack propagation of specimen BCJ2 at displacement of: (**a**) 3 mm; (**b**) 6 mm; (**c**) 12 mm; (**d**) 24 mm.

**Figure 11.** The failure mode of specimen BCJ2.

**Figure 12.** Crack propagation of specimen BCJ3 at displacement of: (**a**) 3 mm; (**b**) 6 mm; (**c**) 12 mm; (**d**) 24 mm.

**Figure 13.** The failure mode of specimen BCJ3.

#### *3.2. Load and Displacement Relationship*

The relationship between the load and displacement measured with transducers mounted on the beam for specimen BCJ1 is shown in Figure 14. The maximum load of this specimen is 73.95 kN at a displacement of 25.74 mm. When displacement was 24 mm, there was an extremely wide crack on the column face with a length has reached along the beam height. This also includes an X-shaped crack in the beam-column joint area, as shown in Figure 8d. Therefore, in the next cycle, the load only reached approximately 51.6 kN, and the specimen failed at the displacement of 32.58 mm. The absence of the anchorage of beam longitudinal reinforcement to column as well as no transverse reinforcement in

beam-column joint region of this specimen led to the easy of crack propagation in this specimen, thereby leading to a small deformation capacity and ductility index. The first yield of beam's longitudinal reinforcement of specimen BCJ1 occurred at a displacement of 12.07 mm. yield of beam's longitudinal reinforcement of specimen BCJ1 occurred at a displacement yield of beam's longitudinal reinforcement of specimen BCJ1 occurred at a displacement

**Figure 14.** The hysteretic load-displacement curve for specimen BCJ1.

placement at the first yield of beam's longitudinal reinforcement of specimen BCJ2 was 13.55 irrups in the joint region and the anchorage of beam's longitudinal reinforcement to relacolumn joint region and the anchorage of beam's longitudinal The hysteretic load-deflection curve for specimen BCJ2 is shown in Figure 15. The displacement at the first yield of beam's longitudinal reinforcement of specimen BCJ2 was 13.55 mm, which is similar to that of BCJ1. This is because both specimens have similar cross-sectional size and reinforcement. Due to the fact that specimen BCJ2 has a stirrup in the beam-column joint region, the maximum load is greater than that of the specimen BCJ1, which is 83.48 kN at a displacement of 23.15 mm. It is important to note that the presence of stirrups in the joint region and the anchorage of beam's longitudinal reinforcement to relatively 560 mm into the column, causing delays in the widening and propagation of cracks, therefore the specimen was able to withstand any suddenly load drop as in the case of BCJ1. As a result, this specimen was able to sustain the maximum load in the next cycle. Furthermore, the specimen failure occurred at a displacement of 52.39 mm with a maximum load at failure of 77.89 kN, which was decreased by 6.7%. Therefore, the presence of transverse reinforcement in the beam-column joint region and the anchorage of beam's longitudinal reinforcement in the column cause the deformation capacity to increase. placement at the first yield of beam's longitudinal reinforcement of specimen BCJ2 was 13.55 irrups in the joint region and the anchorage of beam's longitudinal reinforcement to relacolumn joint region and the anchorage of beam's longitudinal

**Displacement (mm)**

**Figure 15.** The hysteretic load-displacement curve for specimen BCJ2.

Figure 16 shows the hysteretic load-deflection curve for specimen BCJ3, which is a strengthening of the BCJ1 model using ferrocement. Similar to the cases of specimens BCJ1 and BCJ2, BCJ 3 experienced its first yield of longitudinal reinforcing bars at a deformation of 12.00 mm. Installation of wire mesh in the beam-column joint area had an insignificant effect on the maximum load, although it increased the deformation capacity of specimen BCJ3. The maximum load was similar to BCJ1, namely 75.64 kN at a displacement of 16 mm. In addition, this load remained constant until the displacement was 24 mm. There was no increase or decrease in load after a displacement of 16 mm because the specimen had significant number of cracks without any crack propagation and widening. The propagation and widening of cracks was prevented by wire mesh which was installed as a structural strengthening. As a result of the ferrocement strengthening, in the following cycle, the load did not decrease immediately as was the case of BCJ1, rather BCJ3 was able to undergo deformation even without increasing the load till a failure occurred at displacement of 52.04 mm with a load reduction of only 5%, which was equivalent to 71.86 kN. However, in the opposite direction, the maximum load of BCJ3 was 81.68 kN, which exceeded that of BCJ2 relatively 74.56 kN in the same loading direction. Therefore, the strengthening of beam-column joint designed with non-seismic building code using ferrocement improves the deformation capacity and ductility. The deformation capacity and ductility of strengthened beam-column joint even was higher than those of beamcolumn joint designed with current code. However, in this case, as previously reported, a failure occurred due to the delamination of the ferrocement from the old reinforced concrete beam-column joint. Furthermore, it is necessary to establish a good bond between the old concrete and the ferrocement, either by increasing the dynabolt anchor numbers or by providing a bonding adhesive between these old concrete and new mortar. Cases with a higher number of anchors and bonding adhesive between the old concrete and new mortar were not reviewed in this study and need to be further investigated.

**Displacement (mm)**

**Figure 16.** The hysteretic load-displacement curve for specimen BCJ3.

The maximum loads, loads at first crack, loads at first yield of beam's longitudinal The maximum loads, loads at first crack, loads at first yield of beam's longitudinal reinforcing bars, and loads at failure together with their corresponding displacements of all specimens tested in this study are summarized in Table 4. Table 4 shows that the strengthening of beam-column joint with ferrocement did not improve the load carrying capacity, but enhanced the ability to deform after peak load. Therefore, the deformation capacity of strengthened beam-column joint was increased by 60%, which was almost similar to that was designed with current building code. To improve the load carrying capacity as well as the deformation capacity, it is recommended to use high strength heattreated steel with fine grains for wire mesh in the future study since such steel has better

performance in carrying static and dynamic loads [70,71]. This type of steel may be welded to reinforcing bars using welding heat input [72], so that the delamination of ferrocement that occurred in this study may be prevented.

**Table 4.** Maximum loads (Pmax), loads at first crack (Pcr), loads at first yield of reinforcing bars (Py) and loads at failure (Pfail ) and their corresponding displacements (∆). Δ


The enveloped load-displacement curves of the 3 specimens tested in this study were compared as shown in Figure 17. It also shows that strengthening the beam-column joints designed with a non-seismic building code using ferrocement increases the deformation capacity. Therefore it is similar to a beam-column joint designed with a new building code that takes into account the effects of seismic.

**Figure 17.** Envelope load-displacement curves.

#### *3.3. Structural Ductility*

μ Δ Δ Structural ductility is defined as the ability of a structure to undergo large inelastic deformation without experiencing significant loss of strength. In this study displacement ductility index was used for assessing the structural ductility. The displacement ductility index (µ) is defined as the ratio of ultimate displacement (∆u) to yield displacement (∆y) as follows [73]:

$$
\mu = \frac{\Delta\_{\rm u}}{\Delta\_{\rm y}} \tag{1}
$$

μ = ∆y – For specimen BCJ1, ultimate displacement is defined as displacement corresponding to 15% drop of maximum load [37,74–77]. Since the maximum load for specimens BCJ2 and BCJ3 at failure dropped only by 6.75 and 5%, respectively, then the ultimate displacement for those specimens is given as failure displacement.

The yield displacement was assessed by three different methods. For the first method, the yield displacement is defined as the displacement at the first yield of reinforcing bars as presented in Table 4. The second method is based on balance of energy. The detail description on how to calculate the yield displacement based on the balance of energy can be found in the references [68,74,75]. The third method is based on the general yielding. The detail description on how to calculate the yield displacement based on the general yielding can be found in the references [68,74]. Since the yield displacement obtained by those three different methods was almost similar, then the average value was used in calculating the displacement ductility index.

–

Table 5 presents the yield displacement, ultimate displacement, and displacement ductility index of all specimens tested in this study. The table shows that the structural ductility of the ferrocement-strengthened beam-column joint was improved by 91% which was higher than the ductility of beam-column joint designed with the new seismic building code. This significantly ductility improvement is due to the inhibition of the crack propagation by the resistance of the installed wire mesh. As a result, the failure of the beam-column joint, which was originally brittle, becomes more ductile. These results indicate that the structures constructed before the implementation of the seismic building code may be strengthened by using ferrocement. Therefore it is presumed to have a similar deformation capacity and ductility as building structures designed with the new seismic building code. This aids in preventing sudden collapse due to failure of the beam-column joint during an earthquake.

**Table 5.** Displacement ductility index.


#### *3.4. Stiffness Degradation*

Stiffness is one of the parameter that can show the seismic performance of reinforced concrete members which is the ability to resist deformation [75,76]. In this study, the stiffness was calculated as secant modulus of envelope load-displacement curves in positive direction shown in Figure 17. Figure 18 shows the comparison of the stiffness of all beamcolumn joint specimens tested in this study. It further shows that the presence of wire-mesh in the beam-column joint strengthened with ferrocement tends to increase its initial stiffness due to the higher elastic modulus of wire-mesh compared to the elastic modulus of concrete. Since there was no crack and the linear load-displacement relationship at the initial stage of loading, the stiffness was constant until the displacement of around 1.5 mm. The presence of cracks that occurred due to loading leads to stiffness degradation along with an increase in beam displacement, as shown in the figure. The stiffness degradation of all tested specimens is similar and closely related to the existent crack propagation.

**Figure 18.** Degradation of stiffness.

#### *3.5. Energy Dissipation*

The total energy given to a structure under loading is called the input energy. Some of the input energy given to a structure is absorbed (dissipated) by the structure. Energy dissipation is described as the ability of a structure to absorb energy through the yielding process in the plastic hinge region. It is used in designing an earthquake-resistant building structure that is ductile in the plastic hinge area. This, therefore, leads to plastic deformation that occurs before failure. The greater the energy dissipation of a structure, the greater it is able to withstand earthquake loads.

Energy dissipation at each cycle can be calculated from the enclosed area within load-displacement loop at this cycle. The cumulative energy dissipation is calculated by summating energy dissipated in previous cycles [68]. Figure 19 shows the cumulative energy dissipation of all beam-column joint specimens tested in this study as a function of displacement. Furthermore, when there was slight deformation, the energy dissipation of all the specimens was similar. As the displacement increased, the beam-column joint specimen strengthened with ferrocement provided greater cumulative energy dissipation than the unreinforced. The energy dissipation of ferrocement-reinforced specimens was almost similar to the specimens designed with the new seismic building code. Even at the displacement greater than 20 mm, the strengthened beam-column joint had a greater cumulative energy dissipation, which shows the effectiveness of beam-column joint strengthening using ferrocement in resisting earthquake loads.

**Figure 19.** Cumulative energy dissipation.

#### *3.6. Comparison with Previous Studies*

Table 6 presents the comparison of improvement of maximum displacement, ductility, initial stiffness, and energy dissipation obtained from this study and the previous studies [66,68] with the difference in strengthening scheme. This table shows that the different strengthening scheme affects the improvement of deformation capacity. The number layer of wire mesh and its orientation angle affect the improvement of maximum displacement and ductility significantly. Meanwhile, the addition of diagonal reinforcement improves energy dissipation significantly.


**Table 6.** Comparison of improvement of maximum displacement, ductility, initial stiffness, and energy dissipation obtained from this study and the previous studies [66,68].

NA \* = not available.

#### **4. Conclusions**

In this research, the deformation capacity of the reinforced concrete beam-column joint designed with a non-seismic building code was investigated. In addition, structural strengthening was analyzed and compared with the deformation capacity of a similar structure designed with the new seismic building code. Strengthening was carried out using ferrocement provided on both sides of the beam-column joint using 4 layers of wire mesh with a diameter of 1 mm and mesh of 25.4 mm. The reversed cyclic loading test results show that the beam-column joint strengthened with ferrocement improved the deformation capacity and ductility. The beam-column joint, which initially experienced brittle shear failure after being strengthened, increased its ductility index from 2.23 to 4.26. This was greater than the ductility index of the beam-column joint designed with the new seismic building code, which was 3.84. The beam-column joint, which initially failed at a displacement of 32.58 mm after being strengthened with ferrocement, failed at a displacement of 52.04 mm and was almost the same as the displacement, which was designed with the new seismic building code. Moreover, the strengthening also significantly improved its stiffness and energy dissipation. Meanwhile the load carrying capacity of the strengthened beam-column joint was 75.64 kN, a slight higher than that of non-strengthened one which was 73.95 kN, but still lower than that was designed with new seismic building code which was 83.48 kN. The failure of the beam-column joint strengthened with ferrocement occurred due to the delamination of ferrocement from the old reinforced concrete beam-column joint.

**Author Contributions:** Conceptualization, M.Z.A. and A.; methodology, M.Z.A., A., M.A. and M.H.; investigation, M.Z.A.; validation, A., M.A. and M.H.; resources, M.Z.A.; writing—original draft preparation, M.H.; writing—review and editing, M.H.; supervision, S.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data of this research is available upon request to the first author.

**Acknowledgments:** The authors are grateful to Delfian Masrura, Maulana Lirad, Muhammad Farizki, Mahlil, and M. Nasir for the assistance provided in obtaining the experimental activities needed for data collection.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Behavior of RC Beam–Column Joints Strengthened with Modified Reinforcement Techniques**

**Aditya Kumar Tiwary 1 , Sandeep Singh 1 , Jasgurpreet Singh Chohan 2 , Raman Kumar 2 , Shubham Sharma 3, \* , Somnath Chattopadhyaya 4 , Farid Abed <sup>5</sup> and Mislav Stepinac 6, \***


**Abstract:** Using a significant number of transverse hoops in the joint's core is one recognized way for achieving the requirements of strength, stiffness, and ductility under dynamic loading in a column joint. The shear capacity of a joint is influenced by the concrete's compressive strength, the anchoring of longitudinal beam reinforcement, the number of stirrups in the joint, and the junction's aspect ratio. Seismic motion on the beam may produce shear capacity and bond breaking in the joint, causing the joint to fracture. Furthermore, due to inadequate joint design and details, the entire structure is jeopardized. In this study, the specimens were divided into two groups for corner and interior beam–column joints based on the joint reinforcement detailing. The controlled specimen has joint detailing as per IS 456:2000, and the strengthened specimen has additional diagonal cross bars (modified reinforcement technique) at the joints detailed as per IS 456:200. The displacement time history curve, load-displacement response curves, load-displacement hysteretic curve, and load cycle vs. shear stress were used to compare the results of the controlled and strengthened specimens. The findings show that adding diagonal cross bars (modified reinforcing techniques) to beam–column joints exposed to cyclic loads enhances their performance. The inclusion of a diagonal cross bar increased the stiffness of the joint by giving an additional mechanism for shear transfer and ductility, as well as greater strength with minimum cracks.

**Keywords:** beam–column joints; shear capacity; cyclic loading; joint's numerical modeling; interior joint; corner joint; modified reinforcement technique (MRT)

#### **1. Introduction**

The beam–column joint is the crucial zone in a reinforced concrete (RC) frame subjected to large forces during several ground shaking events, and its behavior has a significant influence on the response of the structure. Beam–column joints are the link between horizontal and vertical structural elements, and therefore, the joints are directly involved in the transfer of seismic forces [1]. The strength of the joint's component materials is restricted, and the joint's force-carrying capacity is also restricted. As a result, joints can be severely damaged or even destroyed during an earthquake [2,3]. The primary cause of joint failure is insufficient joint shear strength, which occurs due to insufficient and inadequate reinforcing details at the junction region [4]. Since fixing a fractured joint is challenging, the damage

**Citation:** Tiwary, A.K.; Singh, S.; Chohan, J.S.; Kumar, R.; Sharma, S.; Chattopadhyaya, S.; Abed, F.; Stepinac, M. Behavior of RC Beam–Column Joints Strengthened with Modified Reinforcement Techniques. *Sustainability* **2022**, *14*, 1918. https://doi.org/10.3390/ su14031918

Academic Editor: Francesca Tittarelli

Received: 29 December 2021 Accepted: 4 February 2022 Published: 8 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

level needs to be minimized at the construction stage using a variety of techniques. An earthquake-resistant frame's column should be stronger than its beam. Joint panel stirrups contribute to the confining pressure and shear strength required to prevent early brittle collapses; the movement between the columns and beams can be transmitted correctly with enough transverse reinforcement [5]. However, structural elements constructed against gravity loads or in line with seismic standards in the Mediterranean region usually lack transverse reinforcements in the junction [6–9].

Consequently, beam–column junctions have been recognized among the principal causes of damage in pre-existing reinforced concrete (RC) structures in various studies conducted in the context of prior significant earthquakes. In most countries vulnerable to seismic events, pre-seismic codes do not meet present standards for reinforced concrete structures [10]. Earthquakes that have occurred recently (e.g., M7.4 Oaxaca (Mexico), M7.0 Aegean Sea (Turkey–Greece), M6.4 Croatia) have mostly caused failures of masonry buildings, but RC buildings were also, in some cases, heavily damaged [11–13]. One reason is that the beam–column joints in moment-resisting RC frame constructions have enough shear strength and ductility [14].

A design guideline for achieving the appropriate strength for beam–column joints is included in the existing standards. These specifications include enough anchoring for both beam and column bars traveling through or terminating in the joint area, and appropriate flexural strength for the beam and column to ensure beam-failure mechanisms [15,16]. Tsonos et al. [17] gave an overview of modern design codes (Eurocode family of codes) for the seismic performance of RC beam–column joints and compared them with the older standards. A simplified model for strengthening RC beam–column internal joints is given by Bossio et al. [1], which could be used by designers of new joints to quantify the performance of new structures, as well as by designers of external strengthening of existing joints to compute the benefits of the retrofit and shift the initial failure mode to a more desirable one.

In order to improve the structural behavior of precast beam–column connections, Hanif and Kanakubo [18] studied the use of fiber-reinforced cementitious composite (FRCC) to grout the joint region. Under reversed cyclic stress, two full-scale internal beam–column junctions were evaluated. In the joint region of the first specimen, aramid fibers cementitious composite was used. At the same time, polypropylene fibers cementitious composite was utilized to grout the junction in the second specimen. The findings were compared to previous research (no fibers, PVA fibers, and steel fibers). In contrast, steel fibers had considerably increased shear capacity and the largest hysteretic region. According to the test findings, PVA fibers outperformed others in terms of fracture width. In contrast, steel fibers had considerably increased shear capacity and the largest hysteretic region.

Li et al. [19] studied the cyclic behavior of joints built using prefabricated beams and columns composed of engineered cementitious composite (ECC). Two large-scale joint specimens made of conventional concrete and three large-scale joint specimens constructed of ECC were fabricated and tested under cyclic stresses till failure. One specimen was manufactured monolithically and served as a control. The effects of bar splicer sleeves and the connection location on the load-carrying capacity, failure mode, and ductility of produced joints were evaluated using three alternative assembly strategies. ECC enhanced the load-carrying capacity and ductility of constructed joints, according to the findings. The inclusion of longitudinal bars and splicer sleeves increased the load-carrying capacity but reduced ductility because the failure mechanism shifted from flexural to shear. The cyclic behavior was indifferent to connection location when ECC was employed.

The volume percentage of steel fibers in concrete and the detailing of reinforced steel in external beam–column joints were studied by Oinam et al. [20]. Regardless of transverse reinforcement in the joints, specimens of longitudinal beam bars positioned diagonally in the beam–column joints revealed interfacial shear fractures, according to the test findings. Those with a straight longitudinal bar in beams, on the other hand, showed a flexural plastic hinge away from the joint location. Steel fiber-reinforced concrete (SFRC) in the joint's

region demonstrated outstanding ductility, energy absorption, and consistent hysteresis response, despite increasing the spacing of the hoops in the beam–column joints region.

Ravichandran [21] tested fourteen specimens under cyclic loads, one of which was built according to seismic code IS 13920 [22]. In contrast, the others were built without seismic details according to American Concrete Institute (ACI 318) [23], with HyFRC substituting regular concrete in the joint's location. For each volume fraction, two specimens were cast using the same concrete grade with hybrid fibers (80% steel + 20% polyolefin and 60% steel + 40% polyolefin) (0.5%, 1%, 1.5%, and 2%). The findings revealed that high-strength concrete containing 80% steel and 20% polyolefin improved ductility, energy absorption, and overall strength across all volume fractions. However, when compared to the seismic detail specimen, the hybrid fibers specimen with a volume fraction of 2% (80% steel/20% polyolefin) outperformed the seismic detail specimen in terms of energy absorption capacity and ductility.

Six beam–column knee joint specimens were constructed utilizing five created hybrid synthetic fibers, and one control specimen was evaluated under lateral cyclic stress by Zainal et al. [24]. Ferro-Ultra (F6U3), Ferro-Super (F6S3), Ferro-Econo (F6E3), and Ferro-Nylo (F6N3) were used to cast the hybrid fiber-reinforced concrete (HyFRC) joint area (FFC). According to the findings, the HyFRC joints showed substantial improvements in energy dissipation capacity, stiffness degradation rate, and displacement ductility. Compared to the reference specimens, the F6U3 generated the most augmentation, while the FFC produced the least. All hybrid specimens were numerically simulated using the finiteelement method. The average margin error for peak load capacity, peak load displacement, and maximum displacements was 25.89%, 3.45%, and 0.18%, respectively.

Dehghani et al. [25] built a 3D finite-element model to analyze the impact of employing ECC in various patterns of beam–column connections. The model's validity was determined based on Yuan et al.'s findings [26]. The ECC material enhanced the load-carrying capacity and ductility of the beam–column connections but did not affect their initial stiffness. Furthermore, employing ECC outside of the plastic hinge areas was useless, as most tensile and shear cracks are found within the joints. The findings also showed that ECC was ineffective for preventing the diagonal shear crack in the joint area. Alwash et al. [27] and Bossio et al. [28] studied the corrosion of RC joints when exposed to bending moment and axial forces. Their study was oriented on the loss of integrity, a decrease of load-bearing capacity, stiffness, and serviceability due to corrosion, and the on joints' rehabilitation with the patch repair technique (PRT).

The steel plate energy absorption device (SPEAD) system, presented by Giuseppe Santarsiero et al. [29], is a revolutionary strengthening approach that aims to boost the flexural strength of beam and column components in RC frame constructions. The updated SPEAD model produced a 50-percent increase in strength, as well as a significant reduction in bond-slip effects in the joint panel region. This, in turn, resulted in an increase in ductility, which was good.

The concrete compressive strength, joint aspect ratio, and a number of lateral connections inside the joint are the most critical parameters determining the shear capacity of RC beam–column joints. A modified reinforcing technique to increase the shear capacity of cyclically loaded RC beam–column junctions is a viable option. The primary concerns discovered in the literature examined are the anchoring length requirements for beam bars, the provision of transverse reinforcements, and the involvement of stirrups in shear transmission at the joint. Research evaluating the use of extra cross-inclined bars at the joint core found that the inclined bars contribute a novel method of shear transfer, reducing the risk of a diagonal cleavage fracture at the joint. According to major international standards, diagonal cross bars have no effect on the shear strength of a joint. The goal of this study was to enhance core concrete confinement while avoiding reinforcing congestion in joints. The inclusion of diagonal bars adds another mechanism for shear transmission.

#### *1.1. Interior Beam–Column Joint*

In most instances, the breakdown of inner beam–column junctions triggered the failure mechanism in buildings. The statement mentioned above is proven by a thorough analysis of numerous collapsed or seriously damaged pre-seismic code-designed RCframed structures after moderate or major earthquakes [2,3]. As a result, the weakest connection in existent RC movement-resistant frames was identified as the beam–column joints. These joints' inadequate shear strength has been established as the primary cause of joint failure. This lack in strength is commonly caused by insufficient and poorly specified joint reinforcements [30–32].

Furthermore, joint brittleness develops due to deficient reinforcement, especially the joint's transverse reinforcement causing a reduction in the overall ductility of the construction. The modified reinforcing technique aids in transferring shear or provides an extra shear-transfer mechanism at the joint. The force exerted on bars at the column faces causes the bond force to be dispensed by one of the top beam bars.

$$\text{Bond Force} = \frac{\text{IId2}}{4} (\text{fy} + \text{f's}) \tag{1}$$

where f′ s is the compression steel stress at the far face of the joint, fy is the yield stress of steel, and d is the diameter of steel.

#### *1.2. Corner RC Beam–Column Joint*

The joints are usually situated near the roof level in movement-resistant RC structures. Suppose these joints are simply intended for gravity loads and are built according to pre-seismic regulations. In that case, they may sustain significant harm during seismic events because of insufficient shear strength in the corner beam–column junction [10,33,34]. Internal forces created at corner joints may induce joint failure before the beam or column whichever is weaker—reaches its maximum strength. In earthquake-prone nations such as Japan, Mexico, and China, several approaches for repairing and strengthening corner beam–column joints that earthquakes have damaged have been documented [35].

Only the following requirements may be expected to provide appropriate strength for the corner joint [35]:


$$
\rho \le 6\sqrt{\mathbf{f}\_{\mathbf{c}}^{\prime}/\mathbf{f}\_{\mathbf{y}}} \tag{2}
$$

The stresses are measured in pounds per square inch (psi).

#### **2. Detailing Recommendations for Joints**

The below suggestions are provided in regard to the need for anchorage, confinement, and shear inside the core of joints in earthquake-resistant structures [22]:

i. Anchorage: Due to loss of bond at the inner face of an exterior joint, the development length of the beam reinforcement should be computed from the beginning of the 90◦ bend, rather than the face of the column. In wide columns, any portion of the beam bars within the outer third of the column could be considered for the computed development length. For shallow columns, the use of stub beams will be imperative. A large-diameter bearing bar fitted along the 90◦ bend of the beam bars should be beneficial in distributing bearing stresses. In deep columns, and whenever straight beam bars are preferred, mechanical anchorage could be advantageous. Joint ties should be so arranged that the critical outer-column bars and the bent-down portions of the bars are held against the core of the joint;

ii. Shear Strength: When the computed axial compression on the column is small, the contribution of the concrete shear resistance should be ignored, and shear reinforcement for the entire joint-shearing force should be provided. In exterior joints, only the ties that are situated in the outer two-thirds of the length of the potential diagonal failure crack, which runs from corner to corner of the joint, should be considered to be effective. The joint shear to be carried by the ties is calculated as:

$$\mathbf{A\_{V}} = \frac{1.5 \,\text{VsS}}{\text{dFy}}\text{,}\tag{3}$$

where Vs = joint shear carried by the ties, A<sup>v</sup> = the total area of tie legs in a pair that makes up one layer of shear reinforcement, and d = the beam's effective depth. For preventing the excessive diagonal compression of core concrete, an upper bound for joint shear, usually stated as a nominal shearing force, must be imposed. The value between 10√ f ′ <sup>c</sup> and 11.5<sup>√</sup> f ′ <sup>c</sup> (psi) is recommended for beams;

iii. Confinement: Horizontal tie legs are ineffectual for providing constraint against the concrete core volumetric expansion, while shear reinforcement restricts only the joint's corner regions. As a result, extra confining bars at right angles to the shear reinforcement are required. The distance between these bars should not exceed 150 mm.

The IS 13920:1993 (Ductile Detailing) [22] gives the following provisions:


#### **3. Experimental Program**

Interior and corner beam–column joint specimens were divided into two groups: one with changed reinforcing techniques and another without. For both reinforcement frameworks, the specimens were cast with reinforcement details according to IS 456:2000 [36], including extra diagonal cross-bracing reinforcement at the two faces of the joints for joint confinement, in accordance with the reinforcement category.

#### *Details of Specimens*

The beam and column dimensions at the corner and interior beam–column joints were similar. The columns were 300-mm deep and 400-mm broad, while the beams were 400-mm deep and 300-mm wide. The beam's span was 3000 mm, while the column's height was 3500 mm. Figure 1 shows the inner beam–column joint reinforcement detailing without changed reinforcement techniques, whereas Figure 2a,b shows the interior beam–column joint reinforcement detailing with modified reinforcement techniques at the joint location. The concrete mix consisted of ordinary Portland cement (43 grade), sand passing through a 4.75-mm IS sieve, and coarse aggregate ranging in size from 10 to 18 mm. The concrete cube's compressive strength after 28 days was 25 N/mm<sup>2</sup> . The main reinforcement was made of steel bars with a yield stress of 415 N/mm<sup>2</sup> . To account for the pull-out force, the longitudinal beam bars and cross-bracing bars were given enough development lengths, as required by the code. Inside a steel mold, the specimens were cast horizontally.

**Figure 1.** Reinforcement details of controlled specimens.

All of the specimens were put through their paces with a constant axial load and cyclic loading at the beam's end, as shown in Figure 2c. To duplicate the gravity force on the column, a 0–500-kN hydraulic jack was attached vertically to the loading frame and applied a constant column axial load, as shown in Figure 2d. The external hinge support was attached to one end of the column and anchored to the strong reaction floor, while the other end was restrained laterally by roller support. To apply reverse cyclic loading, two 200-kN hydraulic jacks were employed, one connected to the loading frame at the top and the other to the strong reaction floor. At a distance of 50 mm from the free end of the beam section of the assembly, the cyclic load was applied. In a load-controlled test, the specimen was exposed to an increasing cyclic load until failure. The load increment was set at 1.962 kN. The specimens were equipped with a linear variable differential transformer (LVDT) with a least count of 0.1 mm to measure the deflection at the loading point [10]. Figures 1 and 2a,b illustrate the schematic diagrams of the controlled and reinforced specimens, respectively.

**Figure 2.** (**a**) Schematic diagram of strengthened specimens; (**b**) reinforcement detailing of strengthened specimens; (**c**) schematic diagram of test set-up; (**d**) specimens loading.

#### **4. Numerical Modeling and Analysis of Beam–Column Joints**

In order to appropriately replicate the tested joint sub-assemblies, symmetry boundary conditions are used. Solid 65, Solid 45, and Link8 elements were used to model the beam– column junction. The concrete was modeled with the Solid 65 element, while the hinge support at the base was modeled with the Solid 45 element. There are eight nodes in these elements, each with three degrees of freedom-translations in the nodal x, y, and z directions. The reinforcement was modeled using the Link8 element. There are two nodes in this three-dimensional spar element, each with three degrees of freedom-translations in the nodal x, y, and z directions. For this, the finite-element software ANSYS Workbench V12 was employed [37]. Following that, each material's element specifics are introduced in Tables 1 and 2. The major goal is to stiffen the column to emulate the beam–column junction behavior on the beam under cyclic loads. The finite-element method [30] converts partial differential equations into a series of algebraic linear equations:

$$\mathbf{[K]}\mathbf{(q)} = \mathbf{(F)},\tag{4}$$

where K = stiffness matrix, q = the nodal displacement vector, and F = the nodal force vector.

#### **Table 1.** Concrete characteristics.


**Table 2.** Steel characteristics.


Concrete: An 8-noded solid element, or Solid 65 element, is used to simulate the concrete [38]. Every node in the corner and inner beam–column junction solid elements have translations in the nodal planes x, y, and z with a degree of freedom of three. Therefore, plastic deformation, three-dimensional cracking, and crushing are all possible with this element. The concrete's characteristics are shown in Table 1, below.

Steel: Standard Grade Fe 415 Mpa steel is used for the steel reinforcement in both the corner and interior beam–column connections. A Link8 component characterizes the steel reinforcement. This element necessitates the use of two nodes. Every node has a degree of freedom of three that correlates to the translations of the node's x, y, and z coordinates. Table 2 shows the characteristics of the steel.

In the engineering data utilized for the FE analysis of corner and interior joints, the geometric properties were as stated in Table 3.


**Table 3.** Geometric properties of the corner and interior joints.

ANSYS software's geometry tools model the interior and corner beam–column junction specimens as a 3D model. Figures 3a,b and 4a,b demonstrate the developed geometry as well as usual reinforcing (with steel) details of regulated and strengthened specimens, respectively, where Figure 3a,b represents an interior joint and Figure 4a,b represents a corner joint. In an ideal circumstance, the true binding force between reinforcing steel and concrete needs to be envisioned. However, a perfect bond between the two materials is postulated in this investigation. The Link8 component symbolizes the reinforcing steel bars connected to the nodes of each adjacent solid component of concrete to provide a consistent bonding; hence, both the materials contribute to the same node. A square mesh is used to obtain good results from the Solid65 element. As a result, the meshing is configured to produce square or rectangular mesh segments. This ensures that the dimensions of the components in the concrete support are compatible with the components and nodes in the model's concrete sections. The specimen is modeled with a square concrete element and a mesh size of 50 mm [35].

**Figure 3.** Geometric model and detail of reinforcement for interior joint; (**a**) controlled specimen, (**b**) strengthened specimen.

**Figure 4.** Geometric model and detail of reinforcement for corner joint; (**a**) controlled specimen, (**b**) strengthened specimen.

#### **5. Results Obtained by Numerical Modeling**

FE analysis findings for controlled and strengthened specimens were obtained utilizing ANSYS Workbench [37]. The specimens were subjected to cyclic loads ranging between 0–500 kN. For both specimens, maximum strain, shear stress, and total deformation are compared. According to the results, the above three interior and corner beam–column junction parameters are managed with improved reinforcing methods at the joint region (Tables 4 and 5). At the joint, the deformation that occurred in the controlled specimen is reduced in a strengthened specimen (Figure 5b), and the same is also transmitted to the CB1 beam of the corner beam–column junction, according to the controlled specimen's overall deformation model (Figure 5a).

**Table 4.** Analyzed findings of interior beam–column joint.


**Table 5.** Post analysis findings of corner beam–column joint.


**Figure 5.** Total deformation for (**a**) control specimen with four beams (ICS) and (**b**) strengthened specimen with four beams (ISS).

From the total deformation model of the control specimen of an interior joint (ICS) (Figure 5a), it is observed that the minimum deformation is in the column and the maximum deformation is in the beam with no modified reinforcement techniques. Furthermore, the total deformation after strengthening with the modified reinforcement technique (ISS) was controlled, as shown in Figure 5b.

From Figure 6a, the maximum shear stress for the control specimen (i.e., without MRT) is observed to be equal in all sections of the interior joint. On the other hand, for the strengthened specimen, as shown in Figure 6b, it is found that there is minimum stress in the column and maximum stress in the beam; the stress in the beam–column joint was also controlled using the MRT technique.

**Figure 6.** Maximum Principal stress in the interior joint; (**a**) control specimen with four beams (ICS) and (**b**) strengthened specimen with four beams (ISS).

From Figure 7a, the maximum principal strain for the control specimen on the interior beam–column junction (without MRT) is observed at the center of the column; while the minimum strain at the bottom support of the column is the strengthened specimen from Figure 7b, it is observed that the maximum principal strain for strengthening the specimen (with MRT) is at the support of the beam and the minimum principal strain is in the column. The strain in the joint is controlled by introducing MRT at the joint.

**Figure 7.** Maximum principal elastic strain in the interior joint for (**a**) the control specimen with four beams (ICS) and (**b**) the strengthened specimen with four beams (ISS).

Figure 8a shows that the highest deformation is at the corner joint and the minimum deformation is at the bottom support of the column for the controlled specimen. Figure 8b shows that the total deformation in the corner joint is controlled, compared to the controlled specimen after being strengthened with the modified reinforcement technique.

**Figure 8.** Total deformation in the corner joint for (**a**) the control specimen (CCS) and (**b**) the strengthened specimen (CSS).

From Figure 9a, it is observed that the maximum principal stress for the control specimen in the corner joint (without MRT) is at the support of the beam, and for the strengthened specimen from Figure 9b, the maximum and minimum stress is found at the support of the beam. It can also be observed that the stress in the beam–column junction is controlled.

In Figure 10a, the maximum shear strain for the control specimen of the corner beam– column junction (without MRT) is observed at the support of the beam, while the minimum strain is at the corner joint. For the strengthened specimen, it is found that the strain in the corner joint is controlled by introducing MRT at the joint.

**Figure 10.** Maximum principal strain in the corner joint for (**a**) the control specimen (CCS) and (**b**) the strengthened specimen (CSS).

#### *5.1. Validation of Results*

The results obtained through finite-element analysis for interior and corner beam– column joints were validated with experimental results in load-deformation behavior, as shown in Figures 11a,b and 12a,b, for the control and strengthened specimens. The load-deformation behavior found in the simulation was very similar to the findings of the experimental studies, with the variation of load ranging from 3% to 5%.

**Figure 11.** (**a**) Load-displacement response of controlled specimen (ICS). (**b**) Load-displacement response of strengthened specimens (ISS).

#### *5.2. Load-Displacement Behaviour for Beam–Column Joints*

In Figures 11a and 12a the load-displacement behavior curves are used to compare the results obtained through experimental and finite-element analysis of controlled specimens. The comparison of strengthened specimens of interior and corner beam–column joints are shown in Figures 11b and 12b, respectively. The results obtained through finite-element analysis were in great concurrence with the experimental results. The load-deformation behavior shown in the simulation was extremely close to that observed in experimental studies, with load variations ranging from 3% to 5%. Compared to the controlled specimens,

the strengthened specimens displayed elastic behavior at the beginning stage in both the cases (interior and corner joint). Analysis determines the load-displacement characteristics indicated the better performance of strengthened specimens featuring cross-inclined reinforcement at the junction, which resulted in overall managed deformation and raised the ultimate loading capability compared to the controlled specimen in both types of joints, whether corner or inner.

**Figure 12.** (**a**) Load-displacement response of controlled specimen (CCS). (**b**) Load-displacement response of strengthened specimens (CSS).

Thus, considering the ultimate load-carrying capacities from numerical studies, the specimens with diagonal confining bars (modified reinforcement technique) performed well for both cases of column axial loads. Furthermore, it can be observed that the displacement is more controlled for the ISS and CSS specimens by using cross-inclined bars at the joint than that of the ICS and CCS specimens for both the column axial load cases.

#### **6. Results and Discussions**

Though the reinforcement detailing of structures conform to the general construction code of practice, it may not adhere to modern seismic provisions. Current seismic code specifications for reinforced concrete-framed constructions are frequently deemed unrealistic by structural experts. They lose structural efficiency when a beam–column

junction is subjected to significant lateral stresses, such as strong winds, earthquakes, or explosions. To satisfy the requirements of strength, stiffness, and ductility under cyclic loading, significant percentages of transverse hoops in the cores of joints are required in these locations. Provisions with a high percentage of hoops generate steel congestion, which causes construction difficulties.

On three fronts, researchers are looking into both kinds of beam–column joints (corner and interior).

The factors influencing the behavior of cyclically loaded corner and interior beam– column junctions are examined in the first approach. IS 456:2000 [36] was used to detail the joints. This method measures the maximum principal elastic strain and shear stress and overall deformation of controlled specimens (without MRT) under cyclic loading.

The controlled specimens are strengthened at the joint area in the second approach by using a modified reinforcement method (MRT) on both sides of the column, having a 12-mm diameter crossbar of length 450 mm (as per IS 456:2000) [36] installed. Testing of the joints under the same cyclic loads as the controlled specimens was performed. According to the findings, the improved reinforcing approach boosted the joint's shear resistance capability while simultaneously limiting overall deformation.

Comparison of all the FEM findings of both the control and strengthened specimens of corner and interior beam–column junctions is completed in the third approach. The cyclic response of the corner, as well as the inner beam–column connection, is found to be improved by utilizing updated reinforcing techniques concerning the maximum principal elastic strain and stress as well as overall deformation. The findings are compared using the lateral-loading vs. lateral-displacement curve, loading vs. deflection hysteretic curve, and deflection time history curve along with shear stress vs. load-cycle curve. The findings of the examination of the corner and interior joints are displayed in Table 6.


**Table 6.** Comparison of specimens with and without MRT.

#### *6.1. Hysteretic Behavior of Corner and Interior Beam–Column Junctions*

In consideration of shear capacity and deformation capability, the stress-strain behavior of both beam–column junctions, i.e., corner and interior, is investigated. Hysteretic curves in Figure 13a–d depicts load-displacement equations for fixed and strengthened specimens. The overall deformation in controlled specimens of both the interior and corner junctions (ICS and CCS) is higher than strengthened specimens (ISS and CSS). The loading capacities of ISS and CSS are significantly higher than ICS and CCS, as per the findings of the hysteretic analysis. The ductility is increased without compromising the stiffness. In general, specimens with diagonal crossbars perform better than their conventionally detailed counterparts.

#### *6.2. Shear Stress vs. Loading Cycle Behavior of Joints*

The ductility of specimens reinforced using the cross-inclined bar, as per IS 456:2000, at the corner and interior beam–column junctions outperforms the regulated specimen with no cross-inclined reinforcement. The addition of cross-inclined reinforcement boosted both the ultimate load-bearing capacity and ductility of the interior as well as of the corner junction in both load circumstances (downward and upward), according to the numerical investigation. The inclusion of slanted bars creates a new shear-transmission mechanism. The corner and interior beam–column junctions using the modified reinforcement method (MRT) have better strength, as shown in Figure 14a,b.

**Figure 13.** (**a**) Load vs. total deformation hysteretic graph for the controlled specimen (ICS) (without MRT) and (**b**) the strengthened specimen (ISS) (with MRT). (**c**) Load vs. total deformation graph for the controlled specimen (CCS) (without MRT) and (**d**) the strengthened specimen (CSS) (with MRT).

**Figure 14.** (**a**) Shear stress vs. load cycle for the controlled specimen (ICS) and the strengthened specimen (ISS). (**b**) Shear stress vs. load cycle for the controlled specimen (CCS) and the strengthened specimen (CSS).

#### *6.3. Displacement Time History Curve for Beam–Column Joints*

**Shear Stress (MPa)**

Figures 15a and 16a illustrate the lateral load-displacement time histories curve obtained through numerical analysis for the controlled specimens. Figures 15b and 16b represent the numerical findings that strengthened the specimen's lateral load-displacement time histories. All of the cycles progressed to the push motion after being started with the pull motion. Adopting cross-inclined bars at the joint location to reinforce beam–column joints offers more strength than the controlled specimen in both cases.

**1 2 3 4 5 6 7 8 9 10 11 12 13**

**Shear Stress (CCS) 46.29 46.58 46.82 47.05 47.29 47.77 48.24 48.74 49.19 49.79 50.56 51.33 52.11 Shear Stress (CSS) 0.42 0.87 1.36 1.88 2.40 3.45 4.50 5.55 6.60 7.65 8.70 9.76 10.81**

**Load Cycle**

**Figure 15.** (**a**) Displacement time history of the controlled specimen (ICS) and of (**b**) the strengthened specimen (ISS).

**Figure 16.** (**a**) Displacement time history of the controlled specimen (CCS) and of (**b**) the strengthened specimen (CSS).

#### **7. Conclusions**

The behavior of beam–column joints in RC structures is of great importance for the seismic behavior of a whole structure. Hence, the investigation and research in the field is beneficial, and new methods, techniques, and procedures can help achieve a better understanding of the complex behavior of the joints themselves. Strengthening and upgrading such elements can be completed, e.g., by high-strength steel bars [39], steel jacketing [40], cementitious composites [41], FRP ropes [42], self-centering friction haunches [43], etc. The RC beam–column joints can also be predicted by using modern techniques such as machine learning [44]. In this study, the performance of interior and corner beam–column joints were analyzed through an experimental program, and the results obtained through tests were validated using finite-element software ANSYS. Similar studies were performed by Santarsiero [45]. The following findings may be derived:

i. Based on the present research, the most critical parameters influencing joint shear capacity are the stirrups quantity, the aspect ratio of the joint, the beam longitudinal reinforcement anchorage, and the compressive strength of concrete;


**Author Contributions:** Conceptualization: A.K.T., S.S. (Sandeep Singh) and S.S. (Shubham Sharma); methodology: A.K.T., S.S. (Sandeep Singh), J.S.C., R.K., S.S. (Shubham Sharma) and M.S.; formal analysis: A.K.T., S.S. (Sandeep Singh), J.S.C., R.K., S.S. (Shubham Sharma) and M.S.; investigation: A.K.T., S.S. (Sandeep Singh), S.S. (Shubham Sharma) and M.S.; resources: S.S. (Shubham Sharma), S.C., F.A. and M.S.; writing—original draft preparation: A.K.T., S.S. (Sandeep Singh), J.S.C., R.K. and S.S. (Shubham Sharma); writing—review and editing: J.S.C., R.K., S.S. (Shubham Sharma), S.C., F.A. and M.S.; supervision: S.S. (Sandeep Singh), J.S.C., R.K., S.S. (Shubham Sharma) and M.S.; funding acquisition: M.S. and S.S. (Shubham Sharma). All authors have read and agreed to the published version of the manuscript.

**Funding:** This work receives no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**


#### **References**


## *Article* **Influence of Friction on the Behavior and Performance of Prefabricated Timber–Bearing Glass Composite Systems**

**Vlatka Rajˇci´c , Nikola Perkovi´c \* , Domagoj Damjanovi´c and Jure Barbali´c**

Civil Engineering, Structural Department, University of Zagreb, 10000 Zagreb, Croatia; vlatka.rajcic@grad.unizg.hr (V.R.); domagoj.damjanovic@grad.unizg.hr (D.D.); jure.barbalic@grad.unizg.hr (J.B.) **\*** Correspondence: nikola.perkovic@grad.unizg.hr

**Abstract:** The basic concept of seismic building design is to ensure the ductility and sufficient energy dissipation of the entire system. The combination of wood and bearing glass represents a design in which each material transmits the load, and with the mutual and simultaneous interaction of the constituent elements, it is also earthquake resistant. Such a system has been developed so that the glass directly relies on the wooden frame, which allows the load to be transferred by contact and the friction force between the two of materials. Within the seismic load, friction between glass and wood is an important factor that affects both the behavior and performance of a wood–glass composite system. The set-up system consists of a single specimen of laminated or insulating glass embedded between two CLT elements. The friction force was determined at the CLT–glass contact surface for a certain lateral pressure, i.e., normal force. Friction depends on the way the elements (especially glass) are processed, as well as on the lateral load introduced into the system. Conducted experimental research was accompanied by numerical analyses. Experimental research was confirmed by numerical simulations.

**Keywords:** composites; timber; CLT; load-bearing glass; earthquake; friction; FEM analysis

#### **1. Introduction**

In the current situation of increasingly acknowledging climate change as a threat to our environment and human society, binding agreements have been made during the COP26, taking place in Glasgow in 2021. The building sector has a huge impact and must provide answers on how to tackle climate change, develop a circular economy, and provide a sustainable environment. The building sector should base future technologies on environmentally friendly materials and construction processes. Timber is the leading biobased material and, through newly designed engineered wood-based materials, the material of the future. One innovative engineering wood product, known as cross-laminated timber (CLT), was introduced in the early 1990s in Austria and Germany. Due to its good mechanical properties, good fire resistance as well as advanced durability, and rheological properties, it has been seen as a potential material to replace reinforced concrete in lowand high-rise buildings. On the other hand, there has been significant development and increase in the use of glass as a load-bearing material. Rajˇci´c et al. concluded that loadbearing glass combined with a timber frame represents a load-bearing composite element, which has very good potential for excellent behavior under normal and seismic loads; it is cost-effective, energy-efficient, and aesthetically acceptable [1–3]. The lack of experience and scientific research as well as non-existing standards and codes covering the design of structural components made as composites from laminated glass and laminated timber limit the implementation of such structural elements in practice. Admittedly, there is a national guideline for the design of glass elements [4]. The purpose of these instructions is to seek to provide an overview that is as complete as possible for the various aspects that must be considered in the design, construction, and control of glass elements with regard

**Citation:** Rajˇci´c, V.; Perkovi´c, N.; Damjanovi´c, D.; Barbali´c, J. Influence of Friction on the Behavior and Performance of Prefabricated Timber–Bearing Glass Composite Systems. *Sustainability* **2022**, *14*, 1102. https://doi.org/10.3390/su14031102

Academic Editor: Syed Minhaj Saleem Kazmi

Received: 15 December 2021 Accepted: 13 January 2022 Published: 19 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to verifying their mechanical strength and stability. There are some existing harmonized standards for glass products, for instance [5,6], necessary for their CE marking, but there are no European harmonized standards that can serve as codes needed for the design of glass structures. Therefore, the European Committee for Standardization (CEN) started the preparation of a new code in the series of Eurocodes in order to clarify the design of safe glass-based structures [7]. Currently, they are in the form of technical specifications (CEN/TS 19100-1–CEN/TS 19100-3). Regarding composite elements of cross-laminated timber and laminated glass, extensive research work should be carried out, which will be a basis for implementation in code. Although there is a limited amount of research dealing with timber-glass composites, the need for large transparent surfaces led architects to use such elements. Adding the aesthetic value of timber and glass and the environmental friendliness of materials that can be fully recycled at the end of their life cycles, a new type of structural element was introduced as timber–glass composite structural systems. These types of structures are also built-in earthquake-prone areas (south Europe, Japan, China, USA). There is a significant concern that should be overcome. Generally, the opinion is that the glass has brittle behavior and cannot be used as a structural element in seismic zones. Antolinc, with his research [8], has contributed to the understanding of the behavior of hybrid structural components based on laminated glass and cross-laminated timber frames. Such a structural component, made of a cross-laminated timber frame infilled with load-bearing laminated glass, has a high potential for various applications. It may be used as for façade element timber houses, as a bracing element for newly built or existing frame structures, as a strengthening structural component in existing timber buildings, or as a supporting structural component in historic buildings during and after their retrofitting and restoration.

During the project financed by the Croatian scientific fund "VETROLIGNUM", led by Prof. Vlatka Rajˇci´c, the system was analyzed in terms of load-bearing capacity, stability, seismic performance, energy efficiency, water tightness, and airtightness. Building with wood is very fast and completely suitable for prefabrication in factories. The LCA (cradle to cradle) shows extremely good results in terms of cost-effectiveness and sustainable construction and a reduced CO<sup>2</sup> footprint [9]. Considering the complexity of the wood-bearing glass composite system and the intentions of presenting the most realistic performance and characteristics of such systems, the research is divided into two sections: laboratory testing and research on numerical models.

The contact between glass and wood, as well as a type of connection in the angles of the wooden frame, are the details that need to be given the utmost attention since it is precisely the manner of joining these elements that greatly determines the behavior of the entire composite system [10]. The most usual way of joining load-bearing glass and wooden structures is by adhesives and different types of mechanical fasteners [3]. Using steel mechanical fasteners results in complicated design solutions and details that damage the edge of the glass, which is the most sensitive part of the glass element. Using adhesives can provide a good connection, but that poses the question of the durability of such systems. In addition to said problems, the seismic load causes damage to the structure at the joint positions, and consequently, the failure of the entire load-bearing system [2,3]. In order to solve that problem and maximize energy dissipation during earthquake loading, a system was developed where the load-bearing glass was inserted into the wooden frame without additional mechanical fasteners and adhesives. The system is designed to allow the glass to move freely in a wooden frame while securing glass stability with additional wooden slats. Therefore, the load is transmitted by direct contact, i.e., friction between two elements [3,8,10,11].

During the experimental campaign, 45 cyclings (racking tests) of the composite panels were performed with four different types of connectors in the laminated timber frame. Tests have shown that failure of a composite panel occurs in a corner of a wooden frame [3]. Due to the partly free movement, the load-bearing glass panels remain intact, which is very

important in such a composite system since the wooden elements are easily replaceable, while the mechanical characteristics and properties of the load-bearing glass are retained.

The problem of friction was investigated and discussed in the already mentioned project "Vetrolignum" led by Prof. Vlatka Rajˇci´c Ph.D. The results of the research are presented here.

#### **2. Timber–glass Hybrid Elements: A Brief Literature Overview**

When designing hybrid timber–glass structural systems, special attention should be paid to the appropriate type of connectors to use. The main goal when choosing a connector and type of connection is to avoid the concentration of stress, such as the local crushing of glass on the edges as well as the occurrence of tensile stresses in the glass, which can lead to sudden failure. The usual way to deal with this problem is through the use of soft coating materials on the metal connectors. An overview of possibly used systems for joining laminated glass structural components is contained in Stepinac et al.'s research [12]. Hamm [13] discussed possible solutions for connecting timber with glass in a composite structural component as well as variants of possible practical application in buildings. Niedermaier, in his research [14], presented the connection of timber frame and glass sheets in a hybrid structural panel using two types of adhesives. In the first case, the glass panels were glued to the timber frame with polyurethane and silicone adhesive. In the second case, the bonding was achieved with an epoxy adhesive. The results of the panel racking test are also presented. The test results show that the distribution of tensile stresses and strains depends on the type of adhesive as well as the geometry of the specimens that are tested. Wellersho et al. [15] presented methods to stabilize the building envelope using glazing. A hinged steel frame was used in which a stabilizing glass sheet was inserted in the first case. The second model used the same steel frame but the glass sheet, in this case, was glued with acrylate and polyurethane adhesives. Weller et al. [16], along with Mocibob [17], further continued to examine the behavior of structural components composed of glass sheets inserted in steel frames using different types of adhesives. It was concluded that the thickness of the glass panel is very important because it determines the lateral inplane stiffness of the hybrid structural component. Different authors (Hochhauser et.al, Neubauer, and Winter et al.) [18–20] examined a hybrid panel in which the main timber frame is connected with a timber subframe by screws and the glass is glued to the subframe. An analysis of the in-plane loaded hybrid systems by using mechanical modeling was carried out and described by Cruz et al. in [21].

It was discussed that glass sheets significantly participate in the transfer of horizontal and vertical loads. Additionally, they participate in the prevention of excessive deformations and may substitute the usual type of bracings of steel and timber frames. Blyberg et al. [22], along with Nicklish et al. [23], presented the test results of timberlaminated glass panels where the connection was made by gluing. The authors present the characteristics of adhesives that may be used for structural bonds. A special focus on the analysis of failure mechanisms of timber–glass glued composite wall panels was presented by Ber et al. [24]. Amadio et al. [25] discussed the problem of glass panel buckling. It was analyzed using extended finite-element (FE) investigations and analytical methods for the effect of circumferential sealant joints and metal supporting frames. In [26] Štrukelj et al. presented results of the racking experimental tests of hybrid walls consisting of a timber frame and glass infill connected using polyurethane sealing. Ber et al. [27] used a parametric numerical analysis in their research. The racking stiffness of timber–glass walls is affected by different parameters, and their influence was reported in this study. In [28] Santarsiero et al. analyzed the potential use of glass in earthquake-prone areas as well as the lack of design codes and standards for the design of earthquake-resistant structures designed with glass. This paper concludes that during the design of the earthquake-resistant structures from glass it is necessary to ensure high ductility and dissipation capacity to glass components in buildings. In [29] Bedon et al. reported and demonstrated how the optimal design of glass components can be efficient and beneficial for multiple design configurations. Special mechanical joints were introduced to enhance the dynamic performance of the glass façade. It was reported that well-designed fasteners can introduce additional flexibility and damping capacities when using hybrid panels in a strengthening traditional building. The first published paper with results from the design and analysis of the innovative laminated timber frames infilled by the laminated glass, which is the main subject of this paper, was presented at the WCTE 2012, or the World Conference on Timber Engineering, 2012 [1]. The innovation in this hybrid element is the contact connection between the timber frame and the glass panels without an additional adhesive layer. The research started in 2007, and it was a collaborative research project between the University of Zagreb and the University of Ljubljana. Žarni´c et al., in their research, followed the conclusions of the EU JRC ELSA Italy feasibility study [30]. The cooperation was established within the former CEN TC250/WG3 and the current TC250/SC11 and is still ongoing. Stepinac et al. recently introduced glued-in steel rods as a standard connector because of their wide use all around the world [31]. Since the innovative element showed very good performance, further cooperation on new parts of structural glass codes and the new parts of the timber structure design will be necessary to upgrade Eurocode 8 to introduce such a new type of hybrid structure for retrofitting and strengthening the existing structures made from various materials (masonry, concrete, etc.). Generally, Neugebauer et al. emphasized that emerging laminated materials and hybrid structures are not sufficiently covered in the Eurocodes [32].

#### **3. Prototype of a Multifunctional Wood–Bearing Glass Composite System**

The main objective of the research was to develop, design, and construct a new composite system that will be used as an independent prefabricated structural component for construction in seismically active areas. The solution, in this case, is simplicity, where the desired behavior of the system is achieved by friction, and therefore, without the use of mechanical connectors or adhesives. The purpose of the research was to design composites and construct details of joints that do not adversely affect the load-bearing glass and to develop systems with a high degree of energy dissipation exclusively by friction between wood and glass. Preliminary research shows that certain composite systems can be used in seismically active areas (such as Croatia) [2,3,10,33–35], but system optimization and parameter analysis have not yet been carried out.

In recent years, thanks to the collaboration of the University of Zagreb and the University of Ljubljana, preliminary testing of the wood–bearing glass composite system at monotonous static and cyclic loading has been carried out. The research and development of energy-efficient composite systems are planned in a three-year project entitled VETROLIGNUM (prototype of a multifunctional wood–bearing glass composite system) funded by the Croatian Science Foundation. This project will build on structure dimensioning knowledge and explore new ways to connect load-bearing elements and prepare a study on optimizing certain parts of the panel to maximize energy efficiency. Žarni´c et al. [36] concluded that it is required to build a prototype of the wood–bearing glass composite system, which could be installed in a real building. Additionally, this type of hybrid element can be used as an independent element in the construction of wooden structures, as a temporary or permanent reinforcement, or to stabilize the elements of existing facilities and cultural heritage sites, and as an element for the construction of multi-purpose and adaptive façade systems (Figure 1).

One of the most important parameters when using the modal analysis is the horizontal stiffness of a building. Stiffness and mass determine the structure's vibration periods and hence the influence of an earthquake's frequency content on a structure's response. If a low-rise (only ground floor) building's vibration periods are overestimated (too long) the resulting seismic forces are too conservative. If the periods are underestimated (too short), the results are on the non-conservative side. The situation is just the opposite for higher buildings where the overestimated periods are on the non-conservative side and vice versa. Hence, great care must be taken in assigning the wall's correct stiffness. In the case of a

timber–glass panel, the stiffness is predominantly dependent on the shear and hold-down behavior of joints in frame corners. However, glass infill greatly increases panel stiffness, whereby the influence of friction between timber and glass also needs to be considered. For the static calculation of frame building reinforced with a timber glass panel, the whole system could be replaced with only one element.

**Figure 1.** 3D—the prototype of a composite system: (**a**) timber frame; (**b**) timber—glass system.

‐ ‐ ‐ ‐ ‐ The basic principle of panel installation is to connect the beams of timber frame with horizontal structural elements of the building. A connection derived by angular brackets has a higher bearing capacity than connections in frame angles (single glued-in rod). By such a solution horizontal force is transmitted directly to the panel without compromising the link between panel and frame structure. According to the analysis of the different types of frames with infill, such as concrete or steel frames with different types of infill (masonry infill as well as concrete, steel, and timber panels as infill), there are not many similarities to describe this system. However, certain similarities between the behavior of CLT panels and timber frame composite systems were found, where one of the important parameters is shear stiffness. The shear stiffness *kc, shear* of timber frame connections can be expressed with Equation (1) from [37]:

$$k\_{c,shear} = 4 \cdot K\_{c,shear} + \frac{0.6 \cdot q\_{vert} \cdot L\_{g,loss}^2 \cdot d\_{g,loss} \cdot c}{u\_{slip,Rd}} \tag{1}$$

‐ where *Kc, shear* is shear stiffness of a glued-in rod, *qvert,* is the vertical line load at the top of the panel, *uslip, Rd* is the slip of the weakest connector at the design strength, *c* is the dynamic friction coefficient of timber–glass contact, *Lglass* is the length of the glass panel, and *dglass* is total glass panel thickness. In order to confirm this hypothesis, it is necessary to know each of the above parameters. Because there is a lack of data in the literature on the value of these stiffnesses, as well as friction coefficients, the main goal of further studies is to determine them experimentally [37].

 Based on reverse-cyclic lateral loading tests on structural timber–glass panels [3], the data show a great way of spending seismic energy which contributes to the development of the forces of friction between wood and glass (Figure 2). The failure occurred in the timber frame corner, followed by the friction force between timber and glass, taking over a considerable amount of horizontal load, i.e., the seismic energy was dissipated through the sliding of glass on timber and activation of the joint in the corner of the timber frame.

‐

‐

‐ ‐ **Figure 2.** The relationship between the friction force and the bearing capacity of frame angles in horizontal loading (for a double-glazed panel). Reprinted with permission from Ref. [3]. Copyright 2015 Ph.D. thesis "Spojevi kompozitnih sustava drvo-nosivo staklo u potresnom okruženju" by Asst. Prof. Mislav Stepinac, Ph.D.

#### **4. Materials and Methods**

‐ ‐ ‐ ‐ ‐ Examining the friction between wood and glass is crucial to understanding the operation of the entire timber–bearing glass composite system, in which the glass panel can slide in a wooden frame. It is the sliding, that is, the friction between glass and wood, that is one of the factors that transfers part of the horizontal load [10]. The set-up system consists of asingle specimen of laminated or insulating glass embedded between two CLT elements. The positioning of the glass was achieved by making additional wooden slats that prevent the lateral displacement of the glass but do not press it laterally, and therefore, do not affect the force of friction. Based on the test, the friction force was determined at the wood–glass contact surface for a certain lateral pressure, i.e., the normal force). As a result, a coefficient of friction was obtained, which could be used to numerically model the contact between wood and glass in a calculation model. Numerical analysis was carried out with "Ansys" software support.

#### *Description and Preparation of Specimens*

‐ ‐ The friction testing system consists of CLT elements, glass specimens, and steel elements (lateral force introduction). In order to optimize the system, glass specimens of various types and thicknesses were prepared.

‐ ‐ ‐ ‐ Laminated safety glass is a "sandwich" of two or more glass surfaces that are glued together. "Glue" is a special transparent layer (PVB—polyvinyl butyral, EVA—ethylene vinyl acetate) with a thickness of 1–2 or sometimes more millimeters. In the event of glass breakage, shards and pieces of glass do not scatter but remain retained in the frame thanks to the plastic interlayer. This glass, too, absorbs wide-range sound vibrations and provides better sound insulation than float glass with the same thickness. It is most often single-laminated glass, which does not exclude the possibility of multiple laminating, i.e., joining several glass surfaces between which there is a transparent PVB (or some other) foil. Lamination is performed with PVB (polyvinyl butyral), EVA (ethylene vinyl acetate—transparent or opal), and TPU (thermoplastic polyurethane) foils. In our case, it was lamination with PVB foil.

‐ Insulating glass consists of two or more glass panels that are interconnected at the edge (spacing 6 mm, 9 mm, 12 mm . . . ). The connection allows for a flawless and long-lasting seal, and the interspace is filled with dry air or gas. The distance between the glass plates is provided by aluminum holders that are filled with drying agents. Insulating glass can be produced in combination with tempered or laminated glass. The properties of insulating

glass are that it reduces heat exchange, reduces energy consumption, and does not allow drafts or condensation, so larger glass surfaces can be used for a given room temperature without increasing energy costs.

For this research, 21 glass specimens measuring 200 mm × 400 mm were prepared (Table 1). The specimens were as follows: 3× laminated glass 2 mm × 6 mm, 3× laminated glass 2 mm × 10 mm (Figure 3), 3× insulating (IZO) glass with double-laminated glazing of 6 mm and a cavity width of 12 mm (Figure 4), 3× insulating (IZO) glass with doublelaminated glazing of 10 mm, and a cavity width of 12 mm, 3×—Laminated glass and wooden slat (2 mm × 6 mm) × 2, 3×—Laminated glass and wooden slat (2 mm × 10 mm) × 2, and 3× laminated glass 2 mm × 10 mm smooth ground edges.


**Figure 3.** Laminated glass panel.

**Table 1.** Specimens.

**Figure 4.** Insulating glass panel.

‐ ‐ ‐ All specimens were ESG-toughened glass according to the EN 12150-1 standard [38]. The manufacturing tolerance was within the permissible limits according to the EN 14179-8 standard [39]. The edges of the specimens were roughly sanded [40]. Laminated glass panes were bonded with a 0.76 mm thick PVB membrane. A total of 90% of the cavity in insulating glass was filled with argon. The spacer was 12 mm wide and made of aluminum with respective DC 3363 butyl and silicone layers. Mechanical characteristics of the glass panels can be seen in Table 2.

**Table 2.** Mechanical characteristics of glass.


‐ The timber CLT elements (Figure 5) were processed in the structural testing laboratory of the Faculty of Civil Engineering, University of Zagreb. The CLT consisted of 3 layers, and each layer was 30 mm thick. The timber was class CL24h according to [41]. The additional

wooden beams that support the glass were 30 mm × 30 mm. The material and mechanical properties (acc. to [42,43]) of CL24h timber are shown in Tables 3 and 4.

(**a**) (**b**)

**Figure 5.** CLT specimen: (**a**) wooden slats; (**b**) assembled CLT sample.

ν

**Table 3.** Material properties of CL24h timber.


ν **Table 4.** Mechanical properties of CL24h timber.


#### **5. Experimental Work**

The experiment was carried out by inserting a laminated or insulating glass specimen between two wooden elements. Before starting the experiment and introducing the vertical force, i.e., the force acting in line with the glass pane, it was necessary to secure certain lateral pressure between the glass and the wood. This way we could directly determine the contact point between wood and glass. The lateral force introduction system consists of six steel plates, four threaded rods with nuts, and four springs (Figure 6). ‐

(**a**) (**b**)

‐ **Figure 6.** Specimen set-up: (**a**) laboratory; (**b**) 3D model.

The system, dimensions, and positions of the elements are shown in Figure 7.

The introduction of lateral pressure of the desired amount (1 kN, 2 kN, or 3 kN) was achieved over a certain amount of spring displacement. Springs were positioned between metal plates. The spring displacement itself was achieved through the displacement of the metal plates that push the springs, that is, by controlled tightening and releasing of the nuts on the threaded rod. Such a system allows constant lateral pressure. In order to determine and control the lateral force that was introduced, a preliminary test was conducted to determine the spring stiffness, i.e., to obtain a force-displacement diagram. The diagram represents the force-displacement ratio for all four springs. The spring stiffness can be seen in Figure 8. The stiffness of one spring was determined in such a way that 25% of the amount of force in the diagram was read. The advantage of this system is its simplicity and accuracy. The distance between the two metal plates, that is, the length of the spring, determines the lateral force. The distance was controlled using a sliding caliper with an expanded measurement uncertainty of 20 µm. The simplicity was manifested in the ability to make spacers in desired dimensions that we could place between the two metal plates and then tighten the bolts.

‐ **Figure 7.** Dimension and positions of the test set-up.

‐ ‐

‐

‐

‐

‐

‐

**Figure 8.** (**a**) Spacer position; (**b**) spring stiffness.

‐

μ

μ

‐

‐

‐

After achieving the desired lateral force (Fn) and centering the specimen, the force was introduced to the glass panel. In order to prevent direct contact between the press (steel) and the glass, a thin rubber washer was placed on the edge of the glass, i.e., at the point of load introduction. The load was applied using a universal electromechanical Zwick/Roell testing machine equipped with force sensor class 0.5 in the range from 1 kN to 50 kN according to EN ISO 7500-1:2018 [44] and displacement sensor class 1 according to EN ISO 9513:2012 [45]. The load was applied by displacement control at a speed of 1 mm/min.

The specimen differed in the thickness and type of glass elements (Table 1). Eighteen samples had rough edges, while three samples had smooth ground edges. The sample with smooth ground edges was tested subsequently to see the impact of the glass treatment itself. Each of the samples was tested with a lateral compressive load of 1 kN, 2 kN, and 3 kN. Figure 7 schematically shows the dimensions of the sample and the place of load input. During the experiment, the relative displacement of the glass panels was measured, regarding the fixed CLT elements. Displacement was measured using two LVDTs (Figure 6a) with an expanded measurement uncertainty of 5 µm. The load on the glass panel tangentially to the contact between CLT and glass was measured for a certain normal force Fn. In all experiments, unloading (and then re-loading) was performed in order to eliminate local defects, irregularities, and gaps in the timber material, until the samples fit on the machine surface perfectly, in order to avoid the noise in the data results. The result can be graphically represented as a ratio between the friction force F<sup>t</sup> and the longitudinal displacement at a certain normal force (Fn), as shown in Figures 9–15. The friction force (Ft) is half of the force F required to move the glass panels. The coefficient of friction µ was obtained as the ratio of normal (lateral) force (Fn) and frictional force (Ft).

**Figure 9.** Laminated glass 2 mm × 6 mm.

‐

 ‐

**Figure 10.** Laminated glass 2 mm × 10 mm.

**Figure 11.** IZO glass 2 mm × 10 mm.

**Figure 12.** IZO glass 4 mm × 10 mm. **Figure 12.** IZO glass 4 mm × 10 mm.

**Figure 13.** (2 mm × 6 mm) ×2—Laminated glass and wooden slat.

**Figure 14.** (2 mm × 10 mm) ×2—Laminated glass and wooden slat.

**Figure 15.** Laminated glass 2 mm × 10 mm—smooth ground edges.

#### **6. FEM Research**

‐ The experimental studies carried out were accompanied by numerical analyses. The numerical analyses aimed to extend the knowledge of the behavior of the experimental research. Furthermore, the numerical simulations were performed to confirm and complement the experimental results. The analysis was conducted by Ansys software [46,47]. The entire geometry of the model was drawn by the software "Autodesk Inventor" and imported into "Ansys", where a finite element mesh was formed and in which further simulations were carried out (Figure 16). Element geometry, boundary conditions, loads as well as material characteristics were defined following the experiment. ‐

**Figure 16.** Schematic of numerical analysis procedures.

The model itself is composed of three different materials, namely CLT, glass, and PVB. Boundary conditions and lateral pressure were defined as can be seen in Figure 17.

‐ **Figure 17.** Numerical model: (**a**) Normal force (lateral pressure-1 kN, 2 kN or 3 kN); (**b**) displacement (0 mm-fixed). ‐

To discretize the model, the following Ansys mesh tools [48] were used: "edge sizing" and "sphere of influence". The methods MultiZone (allows the creation of models with a denser grid on the contacts) and Hex Dominant (allows the creation of models where the finite element mesh consists mostly of hexahedrons). Both methods were used to model the glass element (Figure 18).

**Figure 18.** Finite element mesh.

‐ ‐ For modeling contact surfaces, absolute stiffness for normal stresses was used, as well as the possibility of the tangential sliding of two surfaces (CLT and glass) with the corresponding coefficient of friction (Figure 19).

The load was introduced by displacement of the glass panel, according to the steps and data obtained from the laboratory. The result of the experiment, i.e., numerical analysis, is the friction stress that occurs on the contact surface.

#### *FEM Results*

The numerical analysis aimed to obtain a model and certain behavior legality, which would help the prediction of the behavior of such a system during a seismic event. The

main parameter for the control and comparison of numerical simulations and experimental work is the frictional stress that occurs on the contact surfaces. The result obtained from the laboratory was the frictional force required to shift the glass element. In order to compare and evaluate the results of the FEM analysis, the frictional stresses occurring at the contact surfaces were calculated manually, based on the frictional force obtained from the conducted laboratory test. The friction force F<sup>t</sup> is expressed as half the force F required for moving the glass element, as the frictional force occurs on the two surfaces where the glass and the timber connect. In addition to the results in the form of frictional stresses, the behavior of the sample (sliding) was obtained by numerical simulations as shown in Figure 20.

**Figure 19.** Defining contact surfaces.

**Figure 20.** Sample behavior.

For the 2 mm × 6 mm laminated glass sample with a lateral force of 2 kN, the mean frictional stress calculated from the experiment data was 0.25 MPa, while the mean frictional

‐

‐

‐ ‐ stress obtained by numerical simulations was also 0.25 MPa (Figure 21a). The results of the numerical analysis in form of frictional stress can be seen in Figure 21a. For the same type of specimen, but with a lateral force of 3 kN, the maximum frictional stress was 0.38 MPa, and the same value was obtained by numerical simulation (Figure 21b).

‐ ‐

‐

**Figure 21.** Stresses on contact surfaces of laminated glass 2 mm × 6 mm: **(a**) lateral load 2 kN; (**b**) lateral load 3 kN.

In order to confirm the FEM analysis, other types of specimens were subjected to numerical modeling as follows; for the laminated glass specimen 2 mm × 10 mm and a lateral force of 2 kN, the mean frictional stress calculated from the experiment data was 0.1 MPa, while the mean frictional stress obtained by numerical analysis was also 0.1 MPa. The results of the numerical analysis can be seen in Figure 22. For the same specimen, but with a lateral force of 3 kN, the mean frictional stress obtained by the experiment was 0.155 MPa, and the same value of frictional stress (0.155 MPa) was obtained by numerical simulation in Ansys.

**Figure 22.** Stresses on contact surfaces of laminated glass 2 mm × 10 mm: (**a**) lateral load 2 kN; (**b**) lateral load 3 kN.

‐

For all specimen types (lateral force of 2 kN), a comparative analysis is presented. Frictional stresses calculated in Ansys were compared with those obtained from the experimental tests. The maximum deviation in the results was 3.3% (Table 5). Furthermore, at different values of lateral pressure on specimens, an analogy can be established, which confirms the validity of the FEM analysis.


**Table 5.** Comparison of normal stresses between experimental tests and ANSYS.

#### **7. Discussion**

**Figure 23.** Friction load—comparison.

‐ ‐

‐

Based on the presented charts, it was concluded that the friction force increases linearly with increasing lateral force, as expected. Once the legality of the behavior has been determined, a coefficient of friction can be determined for each of the samples. However, to achieve the ultimate goal of the research, it is necessary to highlight and discuss the following findings related to the global behavior of the final product.

‐

‐ ‐

‐

**Figure 24.** Friction coefficient—comparison.

‐ The previous research [36] showed that the influence of glass infills on the lateral load-bearing capacity is significant. Recommendations for further research should be based on the following facts, taking into account the friction between the timber and the glass:


It is possible to formulate this phenomenon with a common equation, which is needed for the definition of the future mathematical model of the tested type of structural hybrid panel components.

Energy dissipation is possible through friction and ductility of the timber frame angle joints. Ductility of the joints in timber structures is a prerequisite, especially in the seismic zones.

#### **8. Conclusions**

Insight into the existing literature and the current state of the art reveals a gap in the study of composite systems with load-bearing glass, especially on loads of horizontal forces of variable amounts and directions that occur during seismic loading. In the range of larger story drifts, the effect of glass-to-timber friction plays a major role in energy dissipation.

During horizontal loading, friction between glass and timber is a factor that affects the behavior of the timber—load-bearing composite system. Coefficients of friction were determined for CLT on glass surfaces; in particular, the effects of different lateral pressure levels were investigated. Friction depends on the way the elements (especially glass) are processed, as well as on the load introduced into the system. The difference between the coefficient of friction at rough and smooth ground edges is negligible. There are differences in the coefficient of friction when insulating glass or glass with wooden slats is installed instead of laminated glass, but it is not significant. The reason lies in the fact that samples with wooden slats have a higher friction surface, and in addition, do not act as a singular system, as is the case of insulated glass.

The investigation provided the necessary data for the development of design procedures and computational model design guidance for the new design codes.

In the future, glass elements with polished edges could be investigated, thus expanding knowledge about the behavior and interaction of these two materials. During load transfer of such a composite system, the contact surface on the wooden element changes and "disappears" over time. Furthermore, future considerations should include how atmospheric factors affect changes in wood surfaces (swelling and shrinkage) and the eventual deterioration of the wood surface, which would cause changes in the contact zone between the two materials, and consequently friction between them. Analysis and research of changes in the coefficient of friction over time and at cyclic loading would be of great importance.

Obtaining realistic values of friction coefficients for different types of glass elements is extremely important for numerical simulations. The use of extreme and theoretical values of friction coefficients in numerical simulations often does not represent a real situation and can lead to wrong conclusions and misinterpretation of results. This research emphasized that the effects of friction should not be neglected. Consequently, neglecting the effects of friction does not unavoidably produce a more conservative design situation by magnifying the stresses. Experimental tests have been confirmed by numerical simulations, but there is the possibility for a more detailed analysis of the system. The numerical analysis should be extended to the whole composite framework and realistic conditions, and thus evaluate all components and factors involved in load transfer

**Author Contributions:** Conceptualization, N.P. and D.D.; methodology, N.P., V.R. and D.D.; software, N.P.; validation, V.R., N.P., D.D. and J.B.; formal analysis, N.P.; investigation, V.R. and N.P.; resources, N.P.; data curation, N.P.; writing—original draft preparation, N.P., V.R., D.D. and J.B.; writing review and editing, N.P., V.R., D.D. and J.B.; visualization, N.P.; supervision, V.R. and D.D.; project administration, V.R.; funding acquisition, V.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** The herein reported research work was funded by The Croatian Science Foundation-Project no.: IP-2016-06-3811 "VETROLIGNUM: Prototype of multipurpose timber-structural glass composite panel".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available on request due to restrictions, e.g., privacy or ethics. The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## **References**


## *Article* **On the Use of Cloud Analysis for Structural Glass Members under Seismic Events**

**Silvana Mattei , Marco Fasan and Chiara Bedon \***

Department of Engineering and Architecture, University of Trieste, 34127 Trieste, Italy; silvana.mattei@phd.units.it (S.M.); mfasan@units.it (M.F.)

**\*** Correspondence: chiara.bedon@dia.units.it; Tel.: +39-040-558-3837

**Abstract:** Current standards for seismic-resistant buildings provide recommendations for various structural systems, but no specific provisions are given for structural glass. As such, the seismic design of joints and members could result in improper sizing and non-efficient solutions, or even non-efficient calculation procedures. An open issue is represented by the lack of reliable and generalized performance limit indicators (or "engineering demand parameters", EDPs) for glass structures, which represent the basic input for seismic analyses or *q*-factor estimates. In this paper, special care is given to the *q*-factor assessment for glass frames under in-plane seismic loads. Major advantage is taken from efficient finite element (FE) numerical simulations to support the local/global analysis of mechanical behaviors. From extensive non-linear dynamic parametric calculations, numerical outcomes are discussed based on three different approaches that are deeply consolidated for ordinary structural systems. Among others, the cloud analysis is characterized by high computational efficiency, but requires the definition of specific EDPs, as well as the choice of reliable input seismic signals. In this regard, a comparative parametric study is carried out with the support of the incremental dynamic analysis (IDA) approach for the herein called "dynamic" (M1) and "mixed" (M2) procedures, towards the linear regression of cloud analysis data (M3). Potential and limits of selected calculation methods are hence discussed, with a focus on sample size, computational cost, estimated mechanical phenomena, and predicted *q*-factor estimates for a case study glass frame.

**Keywords:** seismic design; structural glass; *q*-factor; engineering demand parameters (EDPs); finite element (FE) numerical models; non-linear incremental dynamic analyses (IDA); cloud analysis; linear regression

#### **1. Introduction**

The large use of glass structures in civil engineering applications represents a challenging issue for designers. In addition to intrinsic mechanical features of the involved load-bearing materials [1,2], careful consideration should be paid in earthquake-prone regions to satisfy rigid resistance and displacement demands. This is the case of primary, stand-alone glass structures, but also secondary glass systems belonging to different primary buildings and constructional assemblies [3–8].

According to various literature studies, the seismic capacity of glass structures can benefit from innovative tools and special fasteners [9–11]. At the component level, refined calculation approaches and investigations of literature have been dedicated to both the pre- and post-cracked analysis of laminated glass (LG) elements [12–14], including considerations of their residual strength [15]. In any case, glass structures are still a rather new domain for several professional designers, and certainly require dedicated methods of analysis [16]. Among others, an open issue is represented by the seismic design of glass structures. Most of the available technical documents do not provide specific recommendations for glass [17,18], but suggest the use of "reliable calculation methods" to verify the seismic resistance/displacement capacity of glass components and restraints. Such a

**Citation:** Mattei, S.; Fasan, M.; Bedon, C. On the Use of Cloud Analysis for Structural Glass Members under Seismic Events. *Sustainability* **2021**, *13*, 9291. https:// doi.org/10.3390/su13169291

Academic Editor: Syed Minhaj Saleem Kazmi

Received: 14 June 2021 Accepted: 16 August 2021 Published: 18 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

technical difficulty is further enforced by the need for a realistic calibration of the expected *q*-factor [19].

The main goal of present study, in this regard, is to assess the sensitivity of *q*-factor for glass structures based on simplified or more advanced calculation approaches. As shown in Section 2, consolidated strategies are common for conventional constructional typologies/materials. Moreover, the *q*-factor itself (with *q* > 1) is known to represent the intrinsic dissipation capacity of the structure/material to verify. At the same time, established performance indicators (or "engineering demand parameters", EDPs) in support of seismic analysis and design are available in literature for structural members and systems made of steel, reinforced concrete, timber, or masonry, while such a calibration is missing for glass. To summarize the present discussion, the numerical analysis is focused on a case study glass frame that was earlier investigated in [19]. Differing from [19], however, the attention was given to the seismic performance and capacity of the full-size frame, rather than its key base connections only. To this aim, an original finite element (FE) numerical model was developed and optimized to support the local/global analysis of the frame as a whole. Extended sets of non-linear dynamic analyses were in fact carried out for the frame under in-plane seismic lateral loads. In doing so, three selected methods of analysis that are deeply consolidated for ordinary constructions (M1 to M3 in Section 3) were adapted to the examined structural glass frame and assessed for the *q*-factor prediction. Basic comparative calculations were first carried out with the support of the incremental dynamic analysis (IDA) approach for the herein called "dynamic" (M1) and "mixed" (M2) procedures. The cloud analysis procedure (M3), as shown, is characterized by high efficiency compared to M1 and M2 methods, but requires the calibration of specific EDPs for glass, as well as an accurate selection of input signals for the structural system to verify. FE comparative results are thus discussed in Sections 4–7, showing the potential and limits of selected M1 to M3 calculation methods, in support of a realistic and computationally efficient estimation of seismic behavioral trends for glass structures, thus resulting in their optimized structural design. lished performance indicators (or " arameters", EDPs) in support of nalysis (IDA) approach for the herein called "dynamic" (M1) and "mixed" (M2) –

#### **2. State-of-Art and Literature Review on** *q***-factor Methods**

Following EC8 [17], the seismic design of buildings is today conducted by using the so-called force-based design (FBD) method. The design base shear is conventionally obtained as the ratio between the elastic base shear and the *q*-factor of the structure to verify (Figure 1). The *q*-factor introduction, as such, simplifies its complex energy dissipation capacity (by means of plastic deformations) to a linear elastic model. Due to its strategic role, the *q*-factor definition is thus a topic which has been deeply discussed in the seismic engineering field as a primary focus of several studies since the 1950s.

**Figure 1.** Examples of design spectra calculated for two different *q*-factor values.

A first simple formulation was proposed for the *q*-factor in the 1980s [20]. Further, it was first recognized by the modern design strategy that structures able to resist severe earthquakes are expected to experience permanent damage. As a matter of fact, design

seismic actions are scaled by taking advantage of an intrinsic plastic capacity that is correlated to irreversible deformations.

Typical *q*-factor values by EC8 are known to span in the range of 1.5 (inverted pendulum systems), 2.0 (torsionally flexible systems), and 3.0 (frame systems), or even higher. Due to this, these constant EC8 values should be generally treated as an upper bound, and thus moving the decision on the acceptable level of damage becomes a designer responsibility. This decision is directly affected by structural features, details, and material properties. In such a general discussion, the choice of standards to adopt constant *q*-factor values looks very conservative. Several literature studies proved that the dissipative capacity of a structure is generally greater than the recommended limit values. For example, as concerns steel moment resisting frames (MRFs), the EC8 prescribes different *q*-factors for medium (DCM, local plastic deformations) or high (DCH, global plastic deformations) ductility classes. Macedo et al. [21] evaluated the consequences of adopting the EC8 recommended *q*-factor and presented a more rational selection methodology based on the specific structure and the site seismic hazard. Costanzo et al. [22] discussed existing design provision for both MRFs and chevron concentrically braced frames (CCBFs), giving evidence of a large lateral overstrength due to the codified design requirements. Also for reinforced concrete frames, studies by Kappos [23], Borzi and Elnashai [24], and Chryssanthopoulos et al. [25] assessed the reliability of *q*-factor values by EC8, and emphasized their high conservativity. Although the cited results from [23–25] looked conflicting, the joint EC8 conservatism was jointly justified with either structural overstrength or ductility supply, or both the aspects. In this context, it is thus recognized that the primary goal of standardizing committees is to simplify, on the safe side, the computational burden for designers. Such a strategy makes it possible to avoid performing complex non-linear analyses, but at the same time can severely limit the actual structural plastic capacity of the examined building systems.

For glass structures, to date, legislative and research efforts have not provided a general recommendation about realistic *q*-factor values that could be adopted in design. Furthermore, it is already required to satisfy global and local verifications for resistance and displacement capacities in seismic conditions [3]. As such, the typical effect often takes the form of fully elastic design (*q* = 1). The present study aimed to investigate further the expected structural behavior trends of seismically loaded glass members, based on the observations of a case study frame. Calibrated parameters are presented to possibly support the adaptation of consolidated general procedures to glass structures. The potentials/issues of available methodologies are assessed towards the *q*-factor calculation for similar structural typologies.

#### **3.** *q***-factor and Selected Calculation Methods**

Different approaches can be used from literature to analytically or numerically predict the *q*-factor of a given structural system [26]. As far as the computational effort and accuracy of a method increase, and the reference EDPs are well defined, moreover, the *q*-factor estimation is progressively more robust and reliable. Figure 2 shows a typical push-over (PO) analysis result for a reinforced concrete building, in which the EDPs are qualitatively pointed out, depending on various performance levels and limit states. As usual, the most common EDPs are represented:


**Figure 2.** Expected building response and damage under seismic events. In evidence, the reference limit states and EDPs for design.

– It is worth noting that recommended EDPs are available for traditional constructional materials and systems, but these parameters cannot be directly transferred to glass structures.

It follows that secondary glass members that take place in a primary building must necessarily accommodate the seismic performance and capacity of the building itself (and thus satisfy the corresponding EDPs). For primary/stand-alone glass structures, otherwise, no recommended parameters are available, and thus the present study tries to provide some research developments in this direction.

The above issue arises for novel structural systems and/or innovative materials (glass included), for which no or indicators are provided by design standards for earthquake resistant buildings. Relevant examples of literature can be found in [30–35]. –

From a practical point of view, the flowchart in Figure 3 can be adapted to general constructions/materials, once standardized procedures are established and a primary calculation method is chosen. In case of structural glass (as well as other innovative solutions), the critical step takes place in #2, as a direct/major effect of the analysis method choice and its input basic assumptions (first of all, the reference EDPs).

**Figure 3.** Reference flowchart for the seismic assessment of novel structural systems/materials.

In the present paper, such an issue is further discussed with a focus on the structural glass frame described in Figure 4 and [19]. Three calculation methods (M1, M2, and M3 from Figure 3) are compared in terms of predicted *q*-factor, computational efficiency, accuracy, and sufficiency of results. In doing so, special care is taken for the detection of reliable EDPs that could be used for design, especially with regard to the key configurations of yielding and collapse. The so-called M1, M2, and M3 methods herein explored find inspiration from literature, but in the current study are specifically adapted to glass frames. Examples for traditional structures can be found in [36], as regards the M1 (Section 3.1) and M2 (Section 3.2) methods whilst, for the M3 one (Section 3.3), the procedure in use for the construction of fragility curves is adapted to glass. As such, the M3 *q*-factor is derived from linear regression on a cloud of points that is obtained from non-linear time-history analyses.

**Figure 4.** Reference structural glass frame (adapted from [19]) for the *q*-factor estimation, based on (**a**) non-linear incremental dynamic analyses (IDA) or (**b**) push-over (PO) numerical procedures.

#### *3.1. Dynamic (or PGA) Method (M1)*

"collapse" or "first yielding" respectively: The dynamic (or PGA) method (M1) conventionally defines the *q*-factor as the ratio between PGA*<sup>u</sup>* and PGA*y*, that is the peak ground acceleration values corresponding to "collapse" or "first yielding" respectively:

$$q = \frac{\text{PGA}\_{\mu}}{\text{PGA}\_{y}} \tag{1}$$

In accordance with Equation (1), IDA were thus carried out in this paper. Based on Figure 4a and a set of input accelerograms, sequential non-linear time-history numerical analyses were performed to estimate the PGA*<sup>u</sup>* and PGA*<sup>y</sup>* values of interest. A minimum set of 7 input signals was taken into account [17].

#### *3.2. Mixed Method (M2)*

The mixed (M2) method examined in this paper still takes advantage from efficient FE simulations. The *q*-factor estimation was based in this case on two different contributions, that is:

$$q = \frac{\text{PGA}\_u}{\text{PGA}\_y} \cdot \frac{V\_y}{V\_d} \tag{2}$$

PGA In Equation (2), PGA*<sup>u</sup>* and PGA*<sup>y</sup>* values agree with the definition in Section 3.1, and can be derived from the non-linear IDA for the structural system object of analysis.

At the same time, *V*<sup>y</sup> and *V*<sup>d</sup> in Equation (2) denote the base shear load corresponding to "first significant yield strength" and "allowable design strength". These base shear values were conventionally derived from non-linear PO simulations according to Figure 4b.

#### *3.3. Cloud Analysis (M3) with Linear Regression*

The *q*-factor of the examined frame was finally estimated in this paper by using the inelastic response spectrum, with the support of the so-called cloud analysis and the spectral acceleration (*S*a) definitions [36]. Successful cloud analysis applications can be found in [37–39] for various structural typologies and materials.

Differing from Sections 3.1 and 3.2 (IDA procedure), the cloud analysis is carried out with the support of a set of unscaled accelerograms. The set of input signals (60 in the present study, Table A1 [40]) must be established to ensure an appropriate distribution of cloud data [41,42]. It is in fact known that the major issue of IDA method may consist of a significant computational cost, and most often a marked scaling of original records to various intensity levels, before the desired EDPs could be achieved. This effort is not required in cloud analysis. Moreover, the use of unscaled natural accelerograms in cloud analysis allows to keep all the information related to the event, also known as "record-to-record variability". From unscaled signals and non-linear dynamic analyses, the correlation is established between selected EDPs and some intensity measure (IM) values of the imposed signals by taking advantage of linear regression [37]. As in case of IDA, however, the unscaled signals require an accurate definition of reference EDPs and are also expected to cover a useful range of values for identifying the required limit states. These signals are thus sensitive to the fundamental vibration period *T<sup>1</sup>* (to predict) and the characteristics of the structure to verify (material properties, damage mechanisms, etc.). Moreover, the input signals should be selected to be representative of the seismic hazard of the site under investigation. When appropriate signals are not available, site-specific ground motion modeling techniques can also be used [43–45]. Based on EC8, the *q*-factor can be finally calculated as:

$$q = \frac{\mathcal{S}\_a(T\_1)\_{\
u}}{\mathcal{S}\_a(T\_1)\_{\
y}}\tag{3}$$

where *S*a(*T*1) is the spectral ordinate corresponding to the characteristic period of the design spectrum. The subscripts "*u*" and "*y*" in Equation (3) refer to the "collapse" and "first yielding" configurations.

#### **4. Case-Study Glass Frame**

#### *4.1. Geometrical and Mechanical Properties*

The current research study follows and extends the analytical and numerical investigation reported in [19]. As such, some geometrical and mechanical features are summarized herein for the system in Figure 5a. Each glass frames followed the layout in Figure 5b, with *H* = 6 m and *L* = 8 m. Both the beam and column sections were composed of heatstrengthened (HS) LG members, with uniform size (*h* = 600 mm high × *ttot* = 66 mm thick) given by 5 × *t<sup>g</sup>* = 12 mm glass layers and *tint* = 1.52 mm thick ionoplast foils. The mechanical connection at each beam–column interception took the form of an ideal pin, see Figure 5b. Possible out-of-plane deformations of the frame were restrained, and the related mechanisms (including lateral-torsional buckling for beams [46,47], or coupled bending-compressive buckling for columns [48]) can be preliminarily disregarded. For the base restraints of columns, stainless steel pins pass through two holes in the glass (*ϕ<sup>g</sup>* = 32 mm in diameter, with *ϕ* = 24 mm the nominal diameter of bolts and *D* = 500 mm their distance). Four mild steel brackets (S235 steel) fix the columns to the foundation (*t<sup>s</sup>* = 15 mm, *B<sup>s</sup>* = 200 mm, *b<sup>s</sup>* = 165 mm, *H<sup>s</sup>* = 300 mm and *L<sup>s</sup>* = 200 mm). The restraint was finally locked by *n<sup>b</sup>* anchoring bolts (Figure 5c).

**Figure 5.** Schematic representation of the case study frame: (**a**) design concept (axonometry), with (**b**) static scheme of the structural glass frame object of analysis and (**c**) detail view of a typical push-pull moment connection at the base of the columns (adapted from [19]).

#### *4.2. Preliminary Elastic Seismic Design of the Frame*

For calculation purposes, the glazed assembly of Figure 5 is located in a high seismicity region of Italy. Based on [17,49,50], its strength and stiffness should be verified to resist the most unfavorable expected seismic combination of actions *E*d, that is:

$$E\_d \le R\_d \tag{4}$$

where *R*<sup>d</sup> is the structural capacity.

— — A simple design of glass members should properly verify that LG columns and beams are not subjected—due to the imposed in-plane seismic loads—to relevant stress peaks and

> ∑

=

premature fracture. As for regular structures in plan and elevation, the input seismic force *Fi* the frame should resist is given by:

$$F\_i = F\_b \frac{z\_i m\_i}{\sum\_j^n z\_j m\_j} \tag{5}$$

with *i* = 1, . . . *n*, *m<sup>i</sup>* , *m<sup>j</sup>* the story masses and *z<sup>i</sup>* , *zj* their height from the foundation, while the base shear *F<sup>b</sup>* is:

$$F\_b = \frac{S\_d(T\_1) \ W \lambda}{q} \tag{6}$$

and:


Given the lack of more appropriate recommendations, the conventional design of the case study columns suggests the assumption that *q* = 1. Disregarding the vertical loads that the LG members must sustain (as a part of the framed system of Figure 5), the in-plane lateral force affects the region of glass holes at the base connections, that is:

$$
\sigma\_{t,\max} = \mathbf{K}\_t \cdot \sigma\_t = 2.71 \cdot 77 \approx 208.9 \,\text{MPa} \tag{7}
$$

with:

$$K\_l = 2 + \left(1 - \frac{\phi\_\mathcal{S}}{h/2}\right)^3 = 2.71\tag{8}$$

the magnification factor for stresses [51–53], while the tensile stress *σ<sup>t</sup>* in glass is given by:

$$
\sigma\_l = \frac{F\_l}{\left(\frac{h}{2} - \phi\_\mathcal{g}\right) \cdot t\_{tot}} \approx 77 MPa\tag{9}
$$

with:

$$F\_t = \frac{(0.5 \cdot F\_b) \cdot H}{D} = 1236 kN \tag{10}$$

The resistance verification of the LG columns in seismic conditions requires that:

$$
\sigma\_{\text{t},\text{max}} \le f\_{\text{g;d}} \tag{11}
$$

with:

$$f\_{\mathbb{S}^{\mathcal{A}}} = \frac{k\_{mod}k\_{ed}k\_{sf}\lambda\_{gA}\lambda\_{gl}f\_{\mathbb{S}^{\mathcal{A}}}}{R\_{M}\gamma\_{M}} + \frac{k\_{ed}'k\_{\mathbb{D}}\left(f\_{\mathbb{B}:k} - f\_{\mathbb{S}^{\mathcal{A}}}\right)}{R\_{M;\boldsymbol{\upsilon}}\gamma\_{M;\boldsymbol{\upsilon}}}\tag{12}$$

the design resistance [49]. Among the coefficients in Equation (12), the short-term duration of seismic events (conventionally set in 30 s [3]) suggests *k*mod = 0.78. Given that the columns are composed of HS glass, Equation (12) results in *f* g;d ≈ 75 MPa, that is ≈ 1/3rd the maximum stress from Equation (7), due to the seismic shear from Equation (6).

To avoid the improper sizing of load-bearing glass members, the design (with given input parameters) would require the exploitation of a minimum *q*min ≈ 3. In other words, the ratio of Equation (7) to Equation (12) and combination with Equation (6) can be used for simple analytical estimates of (minimum) required plastic capacities of the frame, towards the seismic demand. In this regard, it is also worth noting that the analytical model developed in [19] for the ductility estimation of base angle brackets (and properly combined with the stress analysis in Equations (7)–(12)), would result in *q* = 4.58 (collapse governed by stress peaks in the region of glass holes). However, such an analytical

prediction is not able to account for complex mechanisms in the frame as a whole, under dynamic seismic accelerations.

#### **5. Finite Element Numerical Investigation**

#### *5.1. Numerical Model*

The reference numerical analysis was carried out in ABAQUS/Standard [54] on a full three-dimensional model representative of the glass frame object of analysis, inclusive of LG members, and reproducing the geometrical details for base connections (Figure 6). For symmetry, 1/4th the geometry was taken into account.

**Figure 6.** Numerical model of the case study structural glass frame under in-plane seismic loads (detail of the base region and beam/column connection, ABAQUS).

Differing from [19], the seismic response of the frame as a whole was explored for the purpose of this study. To this aim, a novel optimized FE model was developed to maximize its computational efficiency.

Solid brick elements (C3D8R type from ABAQUS library) were used for rigid base support, angular members, frictionless foils, and bolts. For the glass plate, a mix of brick solid elements and shell elements was used to preserve the accuracy of stress distributions in the regions of holes. Finally, the column in elevation (and top beam) were described as pinned rigid beams. The overall symmetry assumption resulted in 18,000 solid/shell elements and 75,000 DOFs for the frame model in Figure 6. After preliminary validation, such a solution was used to replace the FE assembly from [19], in which a total of 45,000 solid elements and 170,000 DOFs were used for half geometry of the base connection only.

#### *5.2. Materials and Contact Interactions*

Key mechanical assumptions for glass and steel members were derived from [19]. An elastic-perfectly plastic law was used for mild steel, with *E<sup>s</sup>* = 210 GPa the modulus of elasticity, *ν<sup>s</sup>* = 0.3 the Poisson' ratio, and *σ*s,y = *σ*s,u = 235 MPa the yielding/failure strength, with corresponding strain values equal to *ε<sup>y</sup>* = 0.112% and *ε<sup>u</sup>* = 25%. The ductile damage material option was also accounted for in FE analyses to detect the possible initiation of ductile failure mechanism in angle brackets.

An elasto-plastic law was also used for steel bolts, with *σb,y* = *σb,u* = 1000 MPa the yielding/ultimate resistance (8.8 resistance class). Finally, the rubber layers were described in the form of an equivalent elastic-perfectly plastic material, with *E<sup>r</sup>* = 30 GPa and *ν<sup>r</sup>* = 0.3. The yielding/ultimate stress was conventionally set at 2.4 MPa [19].

The tensile brittleness of glass was included with the concrete damaged plasticity (CDP) model from the ABAQUS library [19]. While the CDP model was primarily developed for concrete, literature studies show that the same model can be efficiently used for structural glass members (under specific loading/boundary conditions), as in the present application. From the post-processing of FE results, fracture initiation in glass was in fact assumed as a reference for "failure". This means that the overall post-cracked stage was disregarded but the simulation was prevented from additional uncertainties that are typical of the post-cracked response of glass under cyclic loads. In doing so, the nominal mechanical properties for HS were taken into account (*E<sup>g</sup>* = 70 GPa, *ν<sup>g</sup>* = 0.23 and *σtk* =70 MPa). Further, the characteristic compressive strength was set to *σck* = 300 MPa (350–500 MPa [19] the reference strength).

A set of surface-to-surface contacts at the interface of adjacent FE components allowed to reproduce the in-plane lateral response of the frame under seismic loads (with "penalty" tangential characteristic (friction *µ* = 0.3) and "hard" normal features). "Tie" mechanical constraints were also used to rigidly connect some FE components (i.e., the head of each bolt and the corresponding angle bracket, or the frictionless layer and the adjacent angle bracket).

#### *5.3. Loading Strategy*

The frame was investigated by taking into account the presence of in-plane seismic loads and dead loads due to constructional members, plus a vertical accidental load Q<sup>k</sup> = 3 kN/m<sup>2</sup> (with *i* = 1 m). The FE assembly of Figure 6 was used for both the required non-linear static PO and time-history dynamic analyses. As such, two different solving procedures were taken into account, based on two separate steps representative of:


In case of IDA, the main seismic input consisted in the selected accelerogram in Figure 7 (acceleration-time history at the base of the frame for M1 and M2 methods). The used earthquake records were derived from [55], that is considering a PGA of 0.35 g, with type A soil (rock soil), topographic category T1, and a reference nominal life of 50 years. A maximum lower and upper tolerance of 10% was also considered.

For the PO analyses (M2 approach), otherwise, the FE system of Figure 6 was subjected to a linear increasing in-plane shear force in accordance to Figure 7 (base connection rigidly fixed to ground). Finally, in case of cloud analysis (M3), a total of 60 unscaled signals (Table A1 [40]) was taken into account to replace the 7 input signals of M1 and M2 procedures.

**Figure 7.** Reference set of time-acceleration histories, as derived from REXEL [55] for the non-linear incremental dynamic analysis (IDA) of the case study glass frame.

#### **6. Discussion of M1 and M2 Results**

*6.1. Detection of First Yielding, Collapse, and Allowable Strength Parameters*

Given the general definitions of calculation approaches in Sections 4.1 and 4.2 and the structural system object of study, special care must be taken for the definition of the key design configurations. Major design challenges derive from the lack of explicit recommendations for glass structures in current design standards for seismic resistant buildings. On the other hand, the in-plane seismic response of the frame is strongly affected by the intrinsic mechanical properties of its basic components, namely:


In other words, the "first yielding" condition of *σ* δ [19]). Regarding the "collapse" damage In other words, the "first yielding" condition of the system of Figure 6 was defined in this project as the first plastic deformation of angle brackets in tension (with *σ*s,y = 235 MPa the reference strength and δ<sup>y</sup> = 0.677 mm the corresponding vertical deformation [19]). Regarding the "collapse" damage state for the frame, maximum drift amplitudes (or column rotations, or even vertical deformations of the steel angle brackets) should be checked. Globally, for the parametric investigation herein summarized, the control of local and global critical conditions for the frame was primarily based on local stress and displacement controls, namely representing a potential:


• In addition, for comparative purposes, conventional deformation limits available in design standards were also taken into account. For the case study frame, as far as the glass holes are properly protected from potential local damage, the seismic analysis could take advantage of the intrinsic flexibility and dissipative capacity of angle brackets. In this sense, the reliability of limit values of Table 1 from FEMA 356 [56], Vision 2000 [57], UBC 1997 [58], EC8 [17], and NTC2018 [50] documents and their applicability to the examined frame were taken into account in this study.

**Table 1.** Recommended limit configurations for the collapse prevention of steel structures, according to selected international design standards.


Finally, the "allowable strength" condition required by the M2 approach should also be defined. As far as the brackets are assumed responsible of the overall in-plane seismic performance of the frame, the yielding stress of steel suggests that: Finally, the "allowable strength" condition required by the M2 approach should

$$V\_d = f(\sigma\_{s,adm}) = \frac{\sigma\_{s,y}}{\gamma\_M} = \frac{235}{1.05} = 223 \text{ MPa} \tag{13}$$

1.05

with γ<sup>M</sup> the partial safety factor, thus a minimum: γ

$$\frac{V\_y}{V\_d} = \gamma\_M = 1.05\tag{14}$$

to account in Equation (2).

Figure 8 shows the base shear-lateral deformation of the frame from PO analysis. It is clear that the above assumption can strongly penalize the frame response, and the critical glass members prove to offer a safety factor in the order of ≈2.1 against potential tensile cracks. Moreover, it is possible to see that the compressive limit in glass holes is not achieved, neither under large displacements. This results from the high deformation capacity of the frame, thanks to detailing of base connections explored in [19]. ical glass members prove to offer a safety factor in the order of ≈2.1

**Figure 8.** PO analysis of the frame, with evidence of relevant EDPs (ABAQUS).

#### *6.2. Seismic Performance Assessment*

The seismic response of the frame was found to agree with Figure 9, where the typical IDA deformed shape (detail) is proposed for brackets under large in-plane lateral displacements.

**Figure 9.** Example of IDA results (ABAQUS): (**a**) deformed shape (extruded axonometric detail, scale factor = 300, at a total time of 15 s); (**b**) stress envelope in the glass holes; (**c**) Von Mises stress in the angle brackets; and (**d**) in-plane lateral deformation. Results for the seismic record a#7 (×1, PGA = 5 m/s 2 ).

– A total of 140 non-linear analyses was carried out with the imposed scaled accelerograms from Figure 7 (with an average of ≈20 differently scaled simulations for each accelerogram). The IDA results still confirmed the close correlation with PO results in Figure 8, with a qualitative agreement of damage phenomena and maximum effects due to the imposed design loads. Figure 9b, in this regard, presents the evolution of maximum stress peaks in the region of holes, while Figure 9c,d focus on the bracket and frame responses, respectively.

from Figure 7 (with an average of ≈20 differently scaled simulations for each accelerogram).

A more detailed analysis of IDA results can be found in Figures 10 and 11, in terms of relevant EDPs, as a function of the imposed PGA for each one of the input scaled signals. It is worth noting that the parametric analysis was carried out in the ideal PGA range of 0–50 m/s 2 to address the performance of structural components. In this regard, typical PGA values can be seen as associated to limited stress levels in the structure, as is expected due to the limited structural mass and high flexibility of the system. Key benefits derive also from gaps in the region of glass holes to prevent premature stress peaks at the edges. Such an approach is also in line with other studies on the seismic performance of structural systems with flexible joints (see for example [59]).

der of ≈

**Figure 10.** Selected IDA numerical results, as a function of the maximum imposed PGA (ABAQUS): (**a**) tensile and (**b**) compressive stress in glass holes, with (**c**) Von Mises stresses and (**d**) plastic strain in the steel angle brackets.

δ θ he order of ≈1 As shown in Figure 10a,b, the compressive stress peaks were mostly observed to double the corresponding tensile peaks in glass, due to a combination of in-plane lateral and vertical loads. Otherwise, it is also interesting to notice that the LG members can sustain relatively strong earthquake motions, before glass could fracture. A relevant aspect is hence represented, in both figures, by the non-linear evolution of stress peaks with the imposed PGA. A first linear trend of the charts can be observed for PGA up to ≈6 m/s 2 , and such a slope change coincides with first yielding (and progressive plastic deformation) of angle brackets. This limit condition was generally achieved for PGA in the order of ≈4 m/s 2 (Figure 10c,d).

≈

Compared to the stress evolution in the holes region, similar trends can also be observed for the deformations of the frame in Figure 11.

The vertical displacement δ, in-plane lateral drift *u* / *H* and base rotation θ are proposed, as obtained from IDA and maximum envelopes of selected EDPs. Under the input assumptions of this study, the collapse condition is never achieved on the side of angle brackets (Figure 11a). The limit drift of 2% or 5% is exceeded for PGA in the order of ≈10 m/s <sup>2</sup> and 20 m/s 2 (average value), see Figure 11b. The 2% drift, finally, is mostly in line with the 0.03 rad rotation of the frame (Figure 11c).

**Figure 11.** Selected IDA numerical results, as a function of the maximum imposed PGA (ABAQUS): (**a**) vertical displacement; (**b**) lateral drift, and (**c**) base rotation of the frame.

#### *6.3. q-factor Estimates*

The analysis was first focused on the so-called M1 method. For the M2 case, the M1 value was adapted with the magnification factor in Equation (14). The so-calculated IDA results are proposed in Figure 12.

"glass"), while in two cases lateral drift of 2% ("Drift 2%"). The corresponding – Note that the attention was focused on the most unfavorable collapse mechanism for the frame as a whole. This was generally observed to coincide with tensile glass cracking ("glass"), while in two cases, only the deformability of the base joint allowed to reach a lateral drift of 2% ("Drift 2%"). The corresponding *q*-factor values are presented for the M1 method (Equation (1)), in the range from 1 (a#1) to 5 (a#6). The average value of *q* = 2.59 is compared with the corresponding M2 estimate from Equation (2), *q* = 2.72. The preliminary analytical requirement (Equations (7)–(12) combined with Equation (6)) is also highlighted (*qmin* = 3), while the analytical value based on local analysis (*q* = 4.58 from [19]) gives evidence of intrinsic limits due to simple predictions carried out for the base connection only (with collapse governed by stress peaks in the region of glass holes).

**Figure 12.** M1 and M2 calculated *q*-factor from IDA (ABAQUS).

#### **7. Cloud Analysis**

*7.1. Input Records and EDPs*

The maximum inter-story drift (IDR) was chosen as reference EDP for the frame (*T*<sup>1</sup> = 0.3 s), with PGA and pseudo-spectral acceleration *S*a(*T*1) being selected as IM parameters. While the PGA value is only related to the seismic ground motion, the *S*<sup>a</sup> value depends on the dynamic behavior of the structure; therefore, the performance to different IMs depends on the type of structure and the governing failure mechanism.

*σ* The preliminary analysis was focused on the distribution of tensile stress peaks in a glass column (with nominal height *H*) deprived of the angle brackets (rigid base connection). Given the characteristic tensile resistance of HS glass (*σtk* = 70 MPa), an inter-story displacement *u*u= 0.048 m was calculated as in Figure 13. The reference EDP corresponds to contour plots in Figure 13a,b), while Figure 13c shows the data trend obtained at the column base and in terms of maximum envelope.

≈ 0.007 and is It is worthy of interest that the so-calculated value corresponds to *u*/*H* ≈ 0.007 and is in close correlation with consolidated limit values for constructional materials characterized by typical brittle behavior in tension, such as, for example, masonry [50]. At the same time, the calculated value significantly minimizes the expected seismic capacity of the frame, thus confirming the key role of its base connections.

The total set of 60 unscaled ground motion records in Table A1 were chosen from [40], depending on the possible collapse mechanism of the frame. According to procedures for general buildings, special attention was paid to cover a wide range of spectral accelerations, but also to respect the consistency between the characteristics of selected records and the supposed classification for the site of interest. As a result, the selected accelerograms were characterized by a moment magnitude (*M*w) between 5.6 and 7.6, an epicentral distances (*R*) ranging between 3.5 km and 62.9 km, and a soil class type A or B (EC8 classification).

**Figure 13.** Definition of the reference inter-story drift, based on the distribution of tensile stress peaks in glass (ABAQUS): (**a**) stress analysis at collapse (legend values in Pa) and (**b**) corresponding in-plane lateral deformation (legend values in m), with (**c**) calculated trends at the column base or from maximum envelope data.

#### *7.2. Analysis of M3 Results with Linear Regression*

The results of the cloud analysis method related to PGA and Sa(*T*1) are separately collected in Figure 14, with attention to the measured inter-story displacement.

Differing from IDA, one of the potential intrinsic limits of the M3 method can manifest in the availability of natural seismic records that possess sufficiently high accelerations to reach the desired EDPs.

In this regard, Figure 15 gives evidence of the typical observed response for the case study frame. As shown, the imposed records are able to lead the angle brackets to yielding (Figure 15a), but still relatively smooth stress peaks are achieved in glass (holes), with tensile and compressive stress peaks in Figure 15b,c. In the same way, the measured in-plane deformations of the frame are still lower than the % limit values earlier discussed for IDA.

**Figure 14.** Cloud analysis results (M3) in the form of inter-story lateral displacement, as a function of (**a**) PGA and (**b**) Sa(*T*<sup>1</sup> ).

**Figure 15.** Cloud analysis results (M3) in the form of (**a**) Von Mises stress in the angle brackets, (**b**) tensile stress, and (**c**) compressive stress in glass (hole region), as a function of Sa(*T*<sup>1</sup> ).

By considering the entire set of available data from the cloud analysis of the frame, an ordinary least-square linear regression was thus performed. The analysis was carried out in the logarithmic space, given that the use of logarithm of variables improves the fit of the model by transforming the distribution of the features to a more normally shaped bell curve. In order to control the skew and counter problems in heteroskedasticity, both the dependent variable (IM) and the independent variable (EDP) were log-transformed. The final result is proposed in Figure 16, where the linear fits of cloud data are obtained from the least squares method.

**Figure 16.** Cloud analysis results (M3) in the form of lateral displacement of the frame, as a function of (**a**) PGA and (**b**) Sa(*T*<sup>1</sup> ).

The thresholds of "yielding" (EDP ) and "collapse" (EDP as "A" and "B") The thresholds of "yielding" (EDPY,50 = 0.037 m) and "collapse" (EDPCP,50) performance levels are indicated in Figure 16 by vertical lines. Two different thresholds (noted as "A" and "B") were used to identify the collapse prevention limit and to quantify further the influence of the base steel connection in the seismic response of the frame, namely:


"B" collapse value of displacement even if it is outside the available data cloud. Once the In this regard, it should be noted that the regression line was assumed to be valid for "B" collapse value of displacement even if it is outside the available data cloud. Once the regression line is found, the IM characterizing yielding and collapse prevention were obtained in Figure 16 using the following relations:

$$IM\_{\rm CP,50} = \exp(a + b \ln(EDP\_{\rm CP,50})) \tag{15}$$

$$IM\_{Y,50} = \exp(a + b \ln(EDP\_{Y,50}))\tag{16}$$

The *q*-factor estimation can thus be based on Figure 16, for PGA and Sa(*T*1), respectively. Certainly, the obtained results are affected by base steel joints and thus by EDPCP,50. As such, the above outcome should be taken into account as a general approach for basic design considerations of similar structures, given that the resistance and stiffness of joints are strictly responsible for the final ductility of the frame, and thus for the possible fracture initiation in glass.

#### *7.3. Comparative q-factor Predictions*

–

In conclusion, Figure 17 shows the M3 calculated *q*-factor and a comparison of selected methods (average). As expected, *q* significantly decreases as far as the M3 approach at collapse disregards the beneficial effect of brackets in the post-yielded stage (collapse "A"). At the same time, as far as the real ultimate inter-story displacement is considered for the tensile fracture of the column (collapse "B"), Figure 17 proves a stable *q*-factor estimation from M3 or M1–M2 methods.

ure of the column (collapse "B"), Figure 17 proves a stable

≈

"A"). At the same time, as far as the real ultimate inter

–

**Figure 17.** Calculated *q*-factor for the examined frame, based on M1 to M3 methods and EDPs.

The positive outcome is confirmation of intrinsic ductility and post-yielding capacity for the frame as a whole. This is in line with preliminary results from [19]. Moreover, the local analysis of angle bracket ductility and stress peak estimation in the region of glass holes was quantified in [19] up to *q* = 4.58 for the frame (collapse governed by tensile fracture of glass). The present study, consequently, confirms the need for full-size structural analyses for special structures and joint details.

Comparative data in Figure 17 are also a confirmation of simple analytical expectations about the minimum plastic capacity of the frame from Section 4.2 (with *q*min ≈ 3 from Equations (7)–(12) combined with Equation (6)), so as to preserve the glass columns from fracture. At the same time, it is necessary to highlight the relatively stable trend for the *q*-factor numerical estimates from the M1 to M3 selected approaches. Most importantly, this finding seems to confirm the potential of linear regression method based on cloud analysis, thanks to the computational efficiency of the M3 method. The M1 and M2 procedures, while limited in number of signals, are univocal in EDPs detection but could require major calculation efforts compared to M3 (140 scaled simulations, in the present study). Furthermore, the IDA calculated average *q*-factor can be highly sensitive to input signals (7 minimum). While the present investigation suggests a very good correlation of M1 to M3 average *q*-factor predictions, this could not be the case of different structural members, thus requiring even more pronounced calculation efforts from IDA (M1 or M2). Finally, compared to simple analytical estimates that are not able to account for complex mechanical phenomena of the frame as a whole (i.e., *q* = 4.58 from [19]), all the numerical estimates in Figure 17 are on the conservative side, thus confirming the need for refined models and non-linear dynamic procedures in support of seismic design.

#### **8. Conclusions**

Available design standards for seismic-resistant buildings provide various recommendations in support of analysis and safe design of several structures subjected to earthquakes, but no specific details are given for glass systems. Among others, major uncertainties derive from the reliable calculation of the seismic performance and dissipation capacity of glass structures, thus their *q*-factor.

In this paper, attention was focused on the local/global seismic analysis of a structural glass frames under in-plane lateral seismic loads. Careful consideration was paid for the development of efficient finite element (FE) numerical models in support of extended parametric non-linear dynamic analyses that could be used to adapt/assess for glass some consolidated procedures in use for structural systems composed of ordinary materials. For most traditional materials and systems, reference engineering demand parameters (EDPs) are recommended by standards or literature documents. On the other hand, reliable EDPs are still lacking for the methods' adaptation to glass structures.

Three numerical calculation methods were taken into account for *q*-factor estimates, based on the parametric incremental dynamic analysis (IDA; dynamic "M1" and mixed "M2" methods), and the cloud analysis based on linear regression ("M3"). Numerical calculations were also compared to simple analytical estimates.

From the FE parametric outcomes, more in detail, it was proven that the metal joints in use for structural glass applications were the major source of possible critical failure mechanisms, but also a key source of enhanced ductility performances for glass members. Such a finding was confirmed in line with ductility and flexibility capacities discussed in [19], based on local analysis of the base connection of the frame. In addition, the present study also confirmed the need of full-size FE models and non-linear dynamic procedures.

In terms of calculated *q*-factor values, more in detail, it was shown that:


At the same time, the adaptation of M3 method with linear regression to structural glass frames:


but on the conservative side compared to simple analytical predictions from [19], based on the local analysis of bracket ductility and stress peaks in the region of glass holes. Such a finding also confirms the need for complex numerical models able to capture dynamic mechanical phenomena in similar systems, as a more detailed investigation to combine with simplified analytical procedures.

**Author Contributions:** Conceptualization, C.B.; methodology, S.M., M.F. and C.B.; validation, C.B.; formal analysis, S.M. and M.F.; data curation, S.M. and M.F.; writing-original, S.M., M.F. and C.B.; writing-review, S.M., M.F. and C.B.; supervision C.B.; project administration, C.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data will be made available upon request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Reference parameters for the selected ground motions records.



#### **Table A1.** *Cont.*

## **References**


## *Article* **Post-Earthquake Damage Assessment—Case Study of the Educational Building after the Zagreb Earthquake**

**Luka Luli´c, Karlo Oži´c, Tomislav Kišiˇcek , Ivan Hafner and Mislav Stepinac \***

Faculty of Civil Engineering, University of Zagreb, 10000 Zagreb, Croatia; luka.lulic@grad.unizg.hr (L.L.); karlo.ozic@grad.unizg.hr (K.O.); tomislav.kisicek@grad.unizg.hr (T.K.); ivan.hafner@grad.unizg.hr (I.H.) **\*** Correspondence: mislav.stepinac@grad.unizg.hr

**Abstract:** In the wake of recent strong earthquakes in Croatia, there is a need for a detailed and more comprehensive post-earthquake damage assessment. Given that masonry structures are highly vulnerable to horizontal actions caused by earthquakes and a majority of the Croatian building stock is made of masonry, this field is particularly important for Croatia. In this paper, a complete assessment of an educational building in Zagreb Lower Town is reported. An extensive program of visual inspection and geometrical surveys has been planned and performed. Additionally, an in situ shear strength test is presented. After extensive fieldwork, collected data and results were input in *3Muri* software for structural modeling. Moreover, a non-linear static (pushover) analysis was performed to individuate the possible failure mechanisms and to compare real-life damage to software results.

**Keywords:** assessment; earthquake; Zagreb; case study; cultural heritage

**Citation:** Luli´c, L.; Oži´c, K.; Kišiˇcek, T.; Hafner, I.; Stepinac, M. Post-Earthquake Damage Assessment—Case Study of the Educational Building after the Zagreb Earthquake. *Sustainability* **2021**, *13*, 6353. https://doi.org/10.3390/ su13116353

Academic Editor: Giacomo Salvadori

Received: 29 April 2021 Accepted: 27 May 2021 Published: 3 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

On 22 March 2020, at 6 h 22 min, Zagreb Metropolitan area was hit by an earthquake of medium magnitude M<sup>L</sup> = 5.5, and intensity of VII, in the epicenter, according to the EMS-98 scale [1]. At 7 h 1 min followed the strongest subsequent earthquake of magnitude M<sup>L</sup> = 5.0 and intensity of VI. The main earthquake damaged most of the buildings in the Lower Town, including residential buildings, universities, schools, kindergartens, hospitals and public buildings. The vast majority of buildings built after the first mandatory earthquake regulations in former Yugoslavia (1964) [2,3] either remained intact or suffered small damage. Nonetheless, the larger part of the city's historical center (Upper and Lower Town) was severely damaged because the buildings in the center were built before any seismic regulations. The damage to historical buildings is enormous. Numerous museums, churches and university buildings have been severely damaged (Figure 1). At the end of the year, Croatia was hit by another devastating earthquake with an epicenter in Petrinja, located approx. 50 km from Zagreb (M<sup>L</sup> = 6.3). The quake caused subsequent damage to already damaged buildings, but to a lesser extent.

As well as most parts of the European region, many existing buildings in Croatia are built in masonry. Given that most of the so-called "strategic" buildings of cultural significance and high historical importance are built using masonry, such a condition is suggesting that the assessment and rehabilitation of existing masonry structures must be conducted on a very high level [4–8]. An important part of the structural assessment is numerical analysis. When it comes to existing buildings, a more refined non-linear analysis should be adopted. Non-linear static analysis or pushover analysis is important and is recommended in Eurocode 8-3 as a reference method for such situations.

**Figure 1.** Typical damage to educational buildings after the Zagreb earthquake (photo credit: M. Stepinac).

called "strategic" buildings of cultural signif- After the earthquake, the first to respond to a disaster were civil engineers who led and coordinated the entire organization of building assessment and damage detection. Various similar post-earthquake assessment procedures are used worldwide [9–11].

– – In the first week, a large number of buildings were inspected, with a rapid postearthquake assessment. The most endangered buildings in the city's center were the ones under cultural heritage protection. The aim of a rapid assessment of buildings is to determine the degree of damage to buildings concerning the protection of life and property, that is, to determine if the buildings are usable, temporarily unusable or unusable. Emerging technological advances allow the usage of artificial intelligence in the postearthquake assessment process in the form of machine learning methods for more efficient and precise results [6,12–15].

quake assessment. The most endangered buildings in the city's – Zagreb's historic urban complex is a protected area regulated by the Law on Zagreb's historic urban complex is a protected area regulated by the Law on the Protection and Preservation of Cultural Heritage. The area is divided into two zones, zone A and zone B (Figure 2) [1]. Zone A includes the oldest and most architecturally valuable parts of Zagreb and is characterized by densely-built blocks of buildings made of stone, brick or a combination of materials. Most buildings consist of massive longitudinal and orthogonal walls and masonry ceiling vaults or wooden ceiling beams and wooden roofs (Figure 3) [1]. Many hospitals, schools, business premises, residential and government buildings, cultural institutions, monuments, churches and chapels are located in zone A and are protected either as part of a historic urban complex or as individual heritage buildings per se. A total of 72% [1] of buildings in zone A suffered major damage due to the earthquake; to compare the damage suffered by this area is almost proportional to the

value of its cultural heritage. Zone B consists of a variety of urban patterns and a large number of immensely valuable buildings [16]. According to the World Bank report, in the educational sector, 106 buildings intended for preschool education, 214 primary and secondary school buildings and education centers, and 12 pupils' dorms were damaged. In the higher education subsector, the damage was reported to 152 buildings. In addition, the buildings of 29 research institutes were also affected. The total value of damage and losses to the education sector is estimated at EUR 1.8 billion at pre-disaster prices, with 97.9% affecting the City of Zagreb [1]. ondary school buildings and education centers, and 12 pupils' dorms were damaged. In ondary school buildings and education centers, and 12 pupils' dorms were damaged. In

**Figure 2.** Protected zones A and B with the location of the case study inside the Lower Town of Zagreb (yellow dashed line).

**Figure 3.** Typical Lower Town masonry building with timber floors and timber roof.

– – In 2020, there were destructive earthquakes all around the world, causing loss of lives, building collapses, and severe structural and non-structural damage and economic losses (e.g., M7.0 Aegean Sea (Turkey–Greece) [17], M6.7 Elazig (East Turkey) [18], M5.5 and M6.4 Croatia [1]). Identification of vulnerability characteristics and earthquake performance assessment of existing structures are essential steps in reducing earthquake losses, and the topic of seismic assessment of existing masonry structures is actual worldwide. Based on the available state-of-the-art literature on assessment and rehabilitation of existing masonry structures (e.g., [19–23]), this paper presents the Croatian perspective and shows it on an actual case study.

#### **2. The Case Study**

The subject of this paper is a building (Figures 4 and 5) located at Vlaška Street 38 as an attached building inside a block. Today's building was built in 1895 by adapting and upgrading two one-story buildings that were built in the early 19th century. The building was upgraded in 1906, while the building's current shape established complete reconstruction in 1997. The building was retrofitted in 1997 for educational purposes and seismic strengthening was not implemented. The building has a rectangular ground plan with the main orientation, which is the longer side of the building in the east-west direction. The building's external dimensions are 12 × 53 m, with two wings: one, 4.4 × 7.6 m and the other, 4.2 × 5.4 m located at the south side of the building. The total floor area of the building is approximately 685 m<sup>2</sup> . The building consists of a basement, first, second, third floor and attic. According to the Croatian seismic hazard map [24], the building is located in the area of peak ground acceleration intensity of agR = 0.255 g for a return period of 475 years. The building serves as an educational-scientific institution. The condition of the building before the earthquake, regarding the vertical loads, was satisfactory, and the building was regularly maintained. located at Vlaška Street 38 as . Today's building was built in upgraded in 1906, while the building's current shape The building's external dimensions are 12 × 53 m, w

–

— **Figure 4.** Aerial view of the building—north façade (photo credit: M. Stepinac).

— **Figure 5.** Aerial view of the building—south façade (photo credit: M. Stepinac).

(unusual for today's standards) used in the

The original drawings show foundations that are approximately 1.0 m wide and 1.15 m deep in relation to the surrounding terrain. They were probably built in brick or stone, which is in line with the construction technology of the time. The building is built of a solid brick of a standard format 30 × 15 × 6.5 (unusual for today's standards) used in the late 19th century. The load-bearing wall thicknesses vary throughout the building, reducing with height, and are 51, 43, 28 cm (Figures 6 and 7). Plaster thickness also varies throughout the building from 3 to 6 cm. The ceiling structures before the reconstruction in 1997 were wooden beams, except for the basement ceiling and the first floor where the masonry vaults are located. After the 1997 reconstruction, the attic and 2nd floor ceilings (Figures 8 and 9) were replaced with reinforced concrete slabs, 12 cm and 16 cm, respectively. The 1st floor ceiling is a semi-precast masonry/concrete floor system (Fert ceiling) inserted between the existing wooden beams. In contrast, the ceilings on the ground floor and basement remained masonry vaults. The building also has two auxiliary staircases at the ends of the building made during the 1997 reconstruction and are made of reinforced concrete, and the main central staircase is older and is made of prefabricated stone stairs supported by a wall on one and a beam on the other edge. (unusual for today's standards) used in the

—

**Figure 6.** Ground floor plan with load-bearing walls in red.

**Figure 7.** 1st floor plan with load-bearing walls in red. The floor plan of the 2nd floor is identical to the 1st floor plan.

**Figure 8.** Longitudinal building section.

**Figure 9.** Transversal building sections.

#### **3. Methodology**

#### *3.1. Assessment Procedure*

The first step in a complete post-earthquake building assessment is a rapid, preliminary assessment of the usability [25,26] of all buildings damaged in the earthquake. Additionally, it is of great importance to preserve the three-dimensional data of the facades of culturally-protected goods in the form of point clouds obtained by laser scanning or drone imaging. The mentioned data can also be used to assess existing structures for the creation of a 3D numerical model. Similar technology was used in the following papers [6,27,28]. In cases where it is needed, detailed assessment and available Non-Destructive Testing (NDT) assessment methods [29] are favorably used. A rapid preliminary assessment is conducted as early as possible after the earthquake, bearing in mind the safety of civil engineers in the field. In Croatia, this type of assessment consisted of a quick visual inspection of individual elements of the load-bearing structure, stating the appropriate degree of damage and deciding on the classification of the building into one of six possible categories (Figure 10): U1 Usable without limitations (Green label), U2 Usable with recommendations

(Green label), PN1 Temporary unusable—detailed inspection needed (Yellow label), PN2 Temporary unusable—emergency interventions needed (Yellow label), N1 Unusable due to external impacts (Red label) and N2 Unusable due to damage (Red label). —

—

**Figure 10.** Six categories of usability divided into three original labels (in Croatian) [30]. •

## *3.2. Rapid Preliminary Assessment Results*

A rapid assessment of the building in question was conducted on 23 March 2020. After a quick visual inspection of individual elements of the load-bearing structure, a decision was made to classify the building as temporarily unusable (Yellow label) with a recommendation for a detailed assessment (PN1). Basic conclusions of the preliminary assessment are: • • •


**Figure 11.** Cracks on the 1st floor: wall, lintel (**a**) and ceiling (**b**) (photo credit: I. Hafner).

The second floor and attic suffered minor damage, while the most severe damage is found on the eastern (Figure 12a,b) and central staircase wings (Figure 13a,b). Recommendations were given that the building can be used with a restriction in the zones where there is a danger of plaster falling. Additionally, the eastern and western staircase can be used with a restriction in the number of people, while the central staircase is not to be used until a detailed assessment is conducted.

**Figure 12.** Cracks on the eastern staircase: exterior (**a**) and interior (**b**) (photo credit: I. Hafner).

**Figure 13.** Diagonal cracks on the central staircase: exterior (**a**) and interior (**b**) (photo credit: I. Hafner).

—

#### *3.3. Detailed Assessment Results*

– — — – — According to the current standards for the design of structures—a series of Eurocodes, HRN EN 1990–1998 and the relevant national annexes, the building that is the subject of this study is in the range of peak ground acceleration of 0.255 g; that is, the expected earthquake intensity is IX according to EMS-98 scale for a return period of 475 years. No geotechnical tests have been performed for the site in question for this article, but based on empirical data, a category B foundation soil (deposits of very compacted sand, gravel or hard clay, at least several tens of meters deep) or category C (deep deposits of compacted or medium—compacted sand, gravel or hard clay with a thickness of several tens of meters to several hundred meters) can be assumed. Moreover, based on the latest findings obtained from the research of the Croatian Geological Institute in cooperation with the University of Zagreb, a seismic microzonation map was prepared according to Eurocode 8 standards for the Zagreb area [31]. According to the mentioned seismic microzonation (2017–2019), the

soil in the immediate vicinity of the assessed building belongs to the category of soil type C. Soil type C causes a certain amplification of the soil shaking, which must be taken into account when assessing the condition of the structure.

All damage, structural and non-structural, is photographed and described. They are plotted in the floor plans of the building (Figures 14–16). The building was inspected from the air by an unmanned aerial vehicle, and no damage was observed to the main load-bearing structure or the building's roof structure. Decorative crosses, statues and reliefs were also inspected. For the purposes of digital preservation, the 3D model of the building was made on the basis of photogrammetric images.

— — **Figure 15.** 1st floor of the building—damage scheme and shear strength testing positions. —

— — **Figure 16.** 2nd floor of the building—damage scheme and shear strength testing positions.

" "

"

" "

" "

"

" "

A detailed inspection of the building revealed the following damage: at the ground floor, the damage is visible in the form of cracks on the wall coverings, arches, vaults and ceilings, as well as separation and local decay of plaster. Cracks in the barrel vaults are mostly parallel to the supporting joints, probably due to lateral movements during the earthquake. They are the result of the occurrence of tensile stresses perpendicular to the supporting joint. Such cracks can cause hinge formation and consequent loss of stability if they propagate deep enough. Fortunately, the cracks in the assessed building are narrow and mostly found in the plaster. Minor local damage to structural elements (walls, columns, arches) is also visible. In the central core of the building where the main staircase is located and in the eastern part of the building, diagonal cracks are visible on the load-bearing walls, which can also be seen on the north side of the building. Damage is visible on all floors in the form of cracks and falling plaster on the walls. Minor local damage to the walls on the 1st floor is also visible, and cracks at the joints of partition walls and ceilings are locally visible. In the central wing of the building where the main staircase is located and in the eastern wing of the building, diagonal cracks are visible on the load-bearing walls, which can also be seen on the north side of the building.

The 2nd floor and attic suffered minor damage. Particular attention should be paid to the central part of the building, occupying the wing with the staircase. The formation of cracks on the central wing transverse walls is clearly visible, indicating a possible failure mechanism out-of-plane of the entire wing. A wedge was made, and the cracks were interconnected and propagated inside the building (they also appear in the stair beams). There was no displacement of the walls out-of-plane, but the preconditions for its failure were met. The central wing needs to be strengthened as soon as possible as part of the entire building's renovation.

The east wing was also damaged at the ground floor and 1st floor level. The cracks that appeared propagated were through the entire wall of the south facade of the wing. It is unfavorable that the cracks are joined and continue to the transversely-joined walls and lintels. The cause of such cracks can be the slight contribution of a torsional response of the building as a whole, where the boundary elements are the most loaded ones, and their failure occurs. Additionally, that part with the building is connected to the neighboring building. Although this can generally have a positive effect on the whole building, in the case of walls on the east wing, such a boundary condition can cause additional forces. If the walls are not well connected to the diaphragms by a tensile compression connection, this can cause them to fail. Since there has been no displacement of the wall out-of-plane, it is not in danger of collapsing, but it should be strengthened soon, and further damage propagation should be prevented. Minor damage in the form of cracks on the walls can be seen on the west and east staircases.

The building can be used in its entirety except for the main staircase. Depending on the possibilities in the future, a static and seismic analysis of the existing condition of the building should be made, and regardless of whether the entire building will be reinforced, the main staircase and other walls with cracks along the entire width of the walls must be repaired and reinforced. Before that, it is necessary to do all the research work to determine the characteristics of the masonry and other necessary data for the structural analysis.

#### *3.4. In Situ Masonry Shear Strength Tests*

In order to assess the condition of the structure after the earthquake and corresponding analysis of the existing condition of the load-bearing structure, in situ tests were carried out. Determination of shear strength (mortar in the composition of load-bearing masonry) of solid brick masonry [32,33] was performed "in situ" using a small hydraulic press "Holmatro" with a load capacity of 200 kN. The mortar was moved horizontally in the vicinity of one brick in order to determine the shear strength. At the same time, the structure of the existing wall was minimally damaged. A total of eleven positions on the load-bearing walls were selected for testing the shear strength of the masonry (Figures 14–16): five positions on the ground floor (PS-PR-1 to PS-PR-5), three positions on the 1st floor (PS-1-1

to PS-1-3) and three positions on the 2nd floor (PS-2-1 to PS-2-3). The shear strength test of the masonry was carried out in nine places. After removing the plaster, it was found that due to the method of masonry (no bricks were found in the longitudinal direction—the inner part of the wall is built of bricks "on edge"), conditions to perform the test were not met for positions PS-PR-3 and PS-PR-4. – — " " – — " "

Results of the Shear Strength Tests

The shear strength of the masonry was obtained based on the registered horizontal force Humax acting on one brick at the time of reaching the shear strength in that brick and the corresponding mortar area on both sides of the shear is transmitted (Ag + Ad). The test method can be seen in Figure 17 and photographs (Figures 18 and 19).

**Figure 17.** Masonry shear strength test method.

credit: L. Lulić) credit: L. Lulić) **Figure 18.** Masonry shear strength tests at measuring positions PS-PR-1 (**a**) and PS-PR-2 (**b**) on the ground floor (photo credit: L. Luli´c).

Lulić)

Lulić)

**Figure 19.** Masonry shear strength tests at measuring positions PS-1-1 (**a**) and PS-1-3 (**b**) on the 1st floor (photo credit: L. Luli´c).

The masonry shear strength test results in the load-bearing walls of the ground floor, 1st floor and 2nd floor are shown in the following Table 1. The values from testing site positions PS-PR-5, PS-1-1, PS-1-2 and PS-2-3 were disregarded because of significant deviations from other results (Figure 20).


**Table 1.** Masonry shear strength test results in load-bearing walls.

— **Figure 20.** fv(shear strength)—σ<sup>0</sup> (vertical stress) diagram.

It is the shear strength f<sup>v</sup> with the contribution of σ0—vertical stress. During the test, each measuring position is precisely located to calculate the vertical load (G0), that is, vertical stress (σ0) from the numerical model. Therefore, for each test site, in addition to the floor plan position, the height of the measuring position from the upper edge of the ceiling structure (h) is recorded. When analyzing the shear strength of mortar, the vertical constant load is taken into account, that is, vertical stresses at a particular test position. Shear strength according to Mohr–Coulomb law is calculated by Equation (1). — 

$$\mathbf{f}\_{\rm vm} = \boldsymbol{\mu} \cdot \boldsymbol{\sigma}\_0 + \mathbf{f}\_{\rm vm0} \tag{1}$$

It can be seen from the diagram that the shear strength without vertical pressure is 0.316 MPa and that coefficient of friction is 1.303. Due to the high value in comparison to EC standard recommendation, the friction coefficient was taken as µ = 0.40 (according to [34]). The obtained shear strength of masonry without vertical pressure, i.e., cohesion, is higher than the one provided by the regulations for the case when there are no tests (fvm0 = 0.32 MPa > 0.10 MPa). The results show that the quality of masonry is good in contrast to similar buildings from that period. – ∙ 

#### *3.5. Numerical Modeling*

The 3D numerical model of the assessed building is obtained in 3Muri software. The macro-element approach is adopted due to computational efficiency and high precision [35]. Its versatility in modeling (implementing elements of various materials, realistic floor stiffnesses, strengthening and many more) makes it highly valuable in a region where a vast majority of building stock is made of masonry. Similar case studies in 3Muri software were used as a base for our research [36–38]. –

Macro-element approach implies equivalent-frame method which uses non-linear beam elements. Macro-elements (or non-linear beam elements) are divided into three categories which are piers, spandrels and rigid nodes. In piers and spandrels, all deformation is concentrated, and they are connected with rigid nodes. Figure 21 shows an equivalent frame model made of mentioned macro-elements.

**Figure 21.** 3D model and 3D equivalent frame in 3Muri.

' Non-linear static pushover analysis [39,40] allows us to check the overstrength ratio used in linear analysis and it gives us more detailed insight into critical elements, possible failure mechanisms, and global behavior of the building as a whole. Pushover analysis is performed with constant gravity loads and monotonically-increasing horizontal loads. Two different distributions of the horizontal loads along the structure's height are used for the pushover analysis. The first distribution has a uniform pattern where the horizontal load is proportional to the mass of the building. The second distribution has a modal pattern where the horizontal load is distributed along with the building's height proportionally to the first vibration mode shape of the building determined through elastic analysis (Figures 22 and 23). These horizontal loads are applied at the location of the masses in the model, i.e., at each floor level in the center of masses. Moreover, accidental eccentricity is taken into account to cover uncertainty in the calculation of the center of masses of the

building. The 5% of the building's length perpendicular to seismic load direction is taken into account on each side for both x (longitudinal) and y (transverse) directions. '

'

**Figure 22.** Mode shape used for pushover in y-direction, T = 0.37 s.

**Figure 23.** Mode shape used for pushover in x-direction, T = 0.24 s.

– – – – – – In each incremental step, internal forces are redistributed according to the element equilibrium, and stiffness is degraded in plastic range. Additionally, ductility is controlled with maximum drift which is 0.004 for shear failure and 0.008 for bending failure [34,41]. Generally, masonry walls have three main failure modes as described in [41–44]. Turnsek– Cacovic constitutive law is used as diagonal cracking is usually the dominant failure mode for existing unreinforced masonry structures [45,46]. In the work of [47], diagonal failure strength is correlated with shear strength used for Turnsek–Cacovic constitutive law. Shear strength obtained by in situ testing is compared with shear strength approximation from the visual MQI method which is explained in more detail in [48] and further developed in [49,50]. Compared shear strengths are very close for this case study which implies good precision of the MQI method. Cracked stiffness of vertical elements is used in a model as recommended in [34], so that cracking that occurs during lifetime because of expected earthquakes of a smaller magnitude is taken into account. Shear and flexural stiffness are taken as half value of initial stiffness.

In a 3D model, floors are modeled as horizontally rigid diaphragms, which is precise enough due to the real in-plane stiffness of the horizontal floor structures. Axial in-plane stiffness of rigid diaphragms in software is infinite and the mass of the real slab is taken into account. Many similar old masonry buildings have an unfavorable distribution of seismic forces due to traditional flexible timber floors [22]. In seismic analysis, the roof is excluded from the load-bearing structure because it does not significantly affect the response of the structure and does not contribute to the global resistance of the structure. Although it was left out of the structural part, its contribution in the form of load on the structure itself was not neglected.

The mean values of material characteristics used in the numerical model (Table 2) are a combination of the literature review [34,51] and on-site testing. Regarding experimental in situ tests and detailed inspection of the structure knowledge level 2 (normal knowledge) can be defined. Based on the achieved knowledge level, confidence factor was taken as 1.2.


**Table 2.** Masonry material characteristics.

According to the work in [52], the building is classified as regular in height but irregular in floor plan, requiring 3D modeling. The building is classified as a torsional stiff system.

First, static analysis is performed according to [53]. Next, the seismic analysis is done. The educational building belongs to importance class III because its seismic resistance is of great importance given the consequences associated with a collapse. Hence, importance factor is γ<sup>I</sup> = 1.2. Three PGA values are used for two limit states. According to the new law "Law on the Reconstruction of Earthquake-Damaged Buildings in the City of Zagreb, Krapina-Zagorje County and Zagreb County (NN 102/2020)" [54], ultimate limit state return period can be different depending on the level of strengthening for old masonry buildings damaged in the recent earthquakes. Limit state of significant damage with a return period of 475 and limit state of damage limitation with a return period of 95 years were checked [55]. In the new law [54], the return period of 225 years which corresponds to a probability of exceedance of 20% in 50 years is introduced for a limit state of significant damage.

Elastic response spectrums for acceleration are calculated for all three return periods taking into account parameters for soil type C, which is found on the location of the building. Altogether, 24 pushover analyses are performed; for x- and y-direction in both orientation, with two load distributions, without and with −/+ 5% of accidental eccentricity.

The result of the performed seismic analysis is a capacity curve that shows the ratio of the shear force in the base of the structure and the displacement of the control node. The control node was selected in the immediate vicinity of the center of mass and is located on the top floor of the building. Obtained capacity curves for all 24 analyses can be seen in Figure 24. Bilinearized pushover curves for the x- and y-direction are shown in Figures 25 and 26. Total base shear in kN is plotted on the y-axis and the displacement of the control nodes in mm is plotted on the x-axis.

**Figure 24.** Pushover curves for the x- (blue) and y- (red) direction.

**Figure 25.** The most relevant pushover curve for the x-direction.

**Figure 26.** The most relevant pushover curve for the y-direction.

Figures 27 and 28 show the damaged state at the last step of the pushover curves for the x- and y-direction where yellow color represents shear damage, and the red color represents bending damage.

**Figure 27.** Damage at maximum displacement capacity for pushover in the x-direction.

a parameter α here α is

eration on type A ground. Parameter α is given for all limit states.

Γ

**Figure 28.** Damage at maximum displacement capacity for pushover in the y-direction.

After the response of the structure, the capacity of the structure is obtained and checks are carried out according to the basic requirements relating to the state of structural damage, defined by limit states. Parameters for equivalent SDOF systems from Figures 27 and 28 are shown in Table 3. These parameters are obtained during bilinearization based on the equivalent energy principal and are used for target displacement determination.


**Table 3.** SDOF parameters for pushover in x- and y-direction.

a parameter α here α is eration on type A ground. Parameter α is given for all limit states. Results are also given in the form of a parameter α (Table 4), where α is the ratio between the limit capacity acceleration of the building and reference peak ground acceleration on type A ground. Parameter α is given for all limit states. A problem arises with old masonry buildings that often cannot be strengthened to such an extent that today's building codes regarding seismic resistance are satisfied. Hence, a new legal document [54] that followed recent earthquakes allows for different levels of earthquake resistance after reconstruction.

**Table 4.** α values.


For the return period of 475 years, none of the 24 analyses satisfies the limit state of significant damage, while for the return period of 225 years, only 8 of the 24 analyses satisfies the same limit state. For a return period of 95 years, none of the 24 analyses satisfies the limit state of limited damage.

The target displacement is determined in accordance with Eurocode 8 (Appendix B) [52]. Since the building's natural period of vibration is less than the period of TC, the target

displacement is determined by the procedure for short periods shown in Figure 29. Buildings with short natural periods of vibration do not comply with the equal displacement rule [56] as buildings with medium-long and long natural periods of vibration.

, Since the building's natural period of vibration is less than the period of T

**α (x α (y**

extent that today's

**Figure 29.** Target displacement determination for short periods.

α values.

Since 3Muri does not consider the bending failure of walls out-of-plane during the global non-linear static analysis, it is necessary to subsequently perform an analysis of the bending of the walls out-of-plane. This takes into account the seismic action perpendicular to the walls. The analysis is performed for the boundary condition near collapse. A return period of 475 years is used. Walls are considered to have passed the out-of-plane bending check if their MRd/MEd ratio is greater than or equal to 1.0. In Figure 30, walls that did not satisfy the check are colored red.

**Figure 30.** Bending out-of-plane results.

The 3Muri program does not consider out-of-plane loss of stability of local mechanisms in the global analysis either. It assumes that proper connection is established between walls and between walls and diaphragms. That way, out-of-plane local mechanisms are prevented so that the global in-plane response of the building can be evaluated [43]. Therefore, the resistance to local losses of stability is checked in a special program module; certain parts of a single wall are checked and the interaction of parts of several individual walls that can together form different local mechanisms. Figure 31 shows an example of a local mechanism where the kinematic block is colored red.

Local mechanisms are arbitrarily defined according to the structure's shape, common failure mechanisms and earthquake damage. The emergence of local mechanisms is often due to the poor interconnection of walls and walls with floor structures. Linear kinematic analysis is used. Defining a local mechanism consists of three steps. To begin with, it is necessary to define a kinematic block that is a part of a wall that is considered absolutely

'

Parameter α repre-

rigid, and that is subject to movement or tilting relative to another block or the rest of the wall. Then, the boundary conditions are defined and finally, the load needs to be set. Some of the possible local failure mechanisms of the observed structure and its resistance to the shown forms of failure are presented below (Figure 32 and Table 5). Parameter α represents the ratio between the spectral acceleration of the activation of the mechanisms and spectral acceleration of the seismic demand.

**Figure 31.** Example of a local mechanism.

**Figure 32.** Local mechanism LM1 (east wing), LM2 and LM3 (central wing).

' **Table 5.** Results for local mechanisms' analysis.


– ' ' After the performed analyses, the real damage was compared with the damage distribution previewed in 3Muri (Figures 33–36). Description of building's capacity through displacement rather than forces allows us to better understand and accurately predict a building's response in the form of a damage initiation and propagation throughout all phases of analysis all the way until the formation of failure mechanism and collapse. That is possible, thanks to incremental non-linear static pushover analysis that shows us damage patterns for every macro-element in every step of the analysis. Compared results are very similar to real damage, which implies good accuracy of the used analysis and software.

(photo credit: T. Kišiček)

–

'

**α**

'

**Figure 33.** First comparison (photo credit: T. Kišiˇcek). (photo credit: T. Kišiček)

'

č **Figure 34.** Second comparison (photo credit: T. Kišiˇcek). (photo credit: T. Kišiček)

č **Figure 35.** Third comparison (photo credit: T. Kišiˇcek). (photo credit: T. Kišiček)

**Figure 36.** Fourth comparison (photo credit: T. Kišiˇcek).

Another important piece of information is that the studied building is a part of a masonry building aggregate, which is not rare in Zagreb and other European cities. The behavior of buildings in the aggregate is extremely complex, where many parameters affect their response during an earthquake. Moreover, buildings at the ends of aggregates are often more severely damaged than those inside aggregates during earthquakes [57,58]. Modeling of the entire aggregate or, at least, adjacent buildings is recommended. Adjacent buildings can be modeled as mutually independent, mutually interconnected by the same walls, or mutually interconnected by various connections that accurately simulate real coherence. Those connections transfer axial and/or shear load between adjacent buildings in aggregate and they can be modeled as linear or non-linear elements [58–60]. Due to the lack of data on neighboring buildings, modeling of only the observed building is often resorted to. This is on the safe side in the case of a unit inside the aggregate when deformations are observed, but there is a risk of misinterpretation of the failure mechanism [61].

č

č

The observed building was modeled as an isolated unit due to insufficient data on neighboring buildings. In order to roughly assess the effect of neighboring buildings on the observed building, the vulnerability index proposed in the articles [62,63] is used. The index is based on five categories depending on the height of the surrounding buildings, the position of the building in the unit, the mismatch of floor heights, differences in material characteristics between buildings and differences in the areas of openings on the facade walls. This method was developed as an upgrade to the vulnerability index based on the ten parameters described in [64]. According to the mentioned five parameters, the buildings in the aggregate have a minimal impact on the observed building, and it can be assumed that the modeling of the building is as isolated as on the safe side.

#### **4. Discussion and Conclusions**

Based on similar case studies from the recent earthquakes in Italy [37,38,65–69], the assessment of the existing educational building is made in order to determine its seismic resistance. The building was damaged in the recent earthquakes in Croatia. The building was not constructed according to the principles of seismic design, but the reconstruction in 1997 partially improved the condition of the structure by adding transverse walls and replacing old wooden beams with reinforced concrete floors. Such rigid diaphragms enable good connection of all walls and thus, better behavior of the building in an earthquake. That is why increasing the stiffness of the traditional timber floors in old masonry buildings

is usually one of the first measures in seismic retrofitting. On the other hand, new research suggests that replacing traditional wooden floors with rigid diaphragms, i.e., RC floors, can induce some unwanted consequences such as cracks on the edges of the two materials or, in the worst scenario, disintegration and collapse of the masonry walls. However, for earthquakes with expected smaller magnitudes as the ones in Zagreb, this measure is convenient and serves to reinforce the existing structure to horizontal actions.

The assessment of the existing unconfined masonry structure was performed using 3Muri with an equivalent frame method. Rather than linear analysis, non-linear static (pushover) seismic analysis is used due to its various benefits. Limit state checks were carried out for the return periods of 95, 225 and 475 years. The results are in line with the expected behavior with respect to the existing wall distribution in the structure. The structure is less rigid and has greater displacements in the y-direction. Additionally, the capacity of the structure is smaller for the y-direction. There is also a small eccentricity between the center of rigidity and the center of mass that causes a slight but negative impact of torsion on the global behavior of the structure. The critical elements are the walls of the central stairwell and the transverse walls on the west side of the building. Additionally, bending failure out-of-plane was checked. Local failure mechanisms were also analyzed using linear kinematic analysis. A comparison of the actual damage with the damage obtained as a result of the conducted non-linear static seismic analysis in the 3Muri program was also conducted.

After assessing the building's behavior in an earthquake, it is obvious that strengthening is needed to raise the seismic resistance level of the building. Damage caused by an earthquake must be repaired to prevent the progression of damage and the possible threat to the global resistance and stability of the building.

As Croatia was hit by two severe earthquakes last year, the knowledge of 'build back better' is fully appreciated. That means that sustainable materials and innovative concepts [70–73] should be used and energy efficiency ensured [74,75]. Various methods are used for strengthening of masonry buildings [76–78] and materials such as FRP and TRM will probably be adopted due to their compatibility and reversibility.

This study provides a detailed insight into the behavior of the building during the earthquake, as well as the extent and distribution of damage and critical elements for earthquakes of different intensities.

**Author Contributions:** Conceptualization: M.S., L.L. and T.K.; methodology: M.S. and T.K.; software: L.L.; validation: M.S. and T.K.; investigation: I.H., L.L. and K.O.; writing—original draft preparation: M.S., K.O., L.L. and I.H.; writing—review and editing: L.L., M.S. and T.K.; visualization: L.L., K.O., I.H. and M.S.; supervision: M.S. and T.K.; project administration: M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Croatian Science Foundation, grant number UIP-2019-04- 3749 (ARES project—Assessment and rehabilitation of existing structures—development of contemporary methods for masonry and timber structures), project leader: Mislav Stepinac.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy reasons.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References and Note**

