**Determination of Young Elasticity Modulus in Bored Piles Through the Global Strain Extensometer Sensors and Real-Time Monitoring Data**

### **Hossein Moayedi 1,2,\*, Bahareh Kalantar 3, Mu'azu Mohammed Abdullahi 4, Ahmad Safuan A. Rashid 5, Ramli Nazir <sup>5</sup> and Hoang Nguyen <sup>6</sup>**


Received: 16 June 2019; Accepted: 12 July 2019; Published: 29 July 2019

**Abstract:** For friction piles depending on the friction resistance, accurate prediction of unit skin friction around the pile shaft is the dominant resistance to measure the final bearing capacity of a bored-pile. The present study measures the stress–strain transferring in two instrumented bored-piles (BP #1 & BP# 2) embedded within the soil layer in Kuala Lumpur by real-time monitoring global strain extensometer (GSE) sensors. Two bored-piles (i.e., having 1.80 m and 1.0 m diameters, as well as 36 m and 32 m lengths) have been loaded with two times to their design working loads. Extensive data are analyzed to measure the changes in stress–strain in the bored-pile. The effect of loading and unloading stages on the pile's head and base settlement has been monitored, indicating that Young modulus of elasticity in concrete bored-pile (*Ec*), average strain, and unit skin friction changed along the bored-pile based on the ground site conditions and stress registered. One example of two case studies with great real-time monitoring data has been provided to further design.

**Keywords:** bored-pile; global strain extensometer; pile friction resistance; real-time monitoring

#### **1. Introduction**

In recent years, calculation of pile bearing capacity data in-situ test has been broadly applied by geotechnical engineers and building foundations, because these data are more accurate and reliable than small-scale laboratory tests. In fact, bored-piles have been considered empirically more as an art-work than science [1], and are formed using appropriate machines (capacity-type) to fill the holes with applicable concrete and reinforcement. Their usual sizes are 400 mm to 3,000 mm diameter, with a capacity that reaches up 45,000 kN based on the pile size and geological profile close to the pile, so an excellent pile capacity has reduced the pile cap size and pile numbers in the group [2,3]. Therefore, bored-pile designing in most countries has relied on the results of the conventional standard penetration test (SPT). According to the literature, the procedure of bored-pile design has consisted of stages such as: (1) the calculation of the end bearing capacity (*fb*) in bored-pile; (2) calculation of the shaft bearing capacity in bored-pile (*fs*), note that, the sum of both values is the ultimate bearing

capacity of an individual pile; and (3) pile working load which has been presumed from the ultimate bearing capacity by using a safety factor that permits the piles' interaction within the group [4,5].

Regarding the empirical approach of *fs* in relation to *Ks* × SPT, *fb* has been related to *Kb* × SPT widely applied in pile designing. To evaluate *Ks* (i.e., skin friction ratio, *fs*/*SPT*) and *Kb* (i.e., end bearing ratio, *fb*/*SPT*), the value including local soil condition has required vibrating wire strain gauges (VWSG) and a mechanical tell-tale rod installed within the pile to measure the axial stress–strain relation and movement in different levels down to the pile's toe and shaft. Bored-pile of length (L) and diameter (Ds) supporting a vertical head load (*Ph*) by the mobilized shaft and base resistance (*Ps* − *Pb*) is illustrated in Figure 1a. The vertical displacement of the pile's base and head are defined as Δ*<sup>b</sup>* − Δ*h*, followed by the pile compressing (shortening) as *ep* ignoring the pile weight as it is in Equation (1) [6]:

$$P\_h = P\_s + P\_b \tag{1}$$

According to the cross-sectional area of concrete *A*, the stress along with the pile σ*<sup>h</sup>* is as Equation (2) [6]:

$$
\sigma\_{\rm li} = \frac{P\_{\rm li}}{A} \tag{2}
$$

$$
\Delta\_h = \varepsilon\_p + \Delta\_b \tag{3}
$$

The length δ*<sup>z</sup>* has been located at depth *z* below the pile head level as shown in Figure 1a–c. Shaft resistance δ*Ps* mobilized along the segment δ*<sup>z</sup>* equals the axial force changing in axial force δ*Pz*. The unit shaft resistance *fs*(*z*) mobilized in δ*<sup>z</sup>* has been related to *P*(*z*) as shown in Equation (4):

$$\frac{dP(z)}{dz} = -\pi D\_s f\_s(z) \tag{4}$$

The negative sign shows that *Pz* decreased as *z* increased (Figure 1c). Omer et al. [7] mention that if the pile's elastic feature has been denoted (*Ep*), ignoring all vertical soil moving and pile displacement *w*(*z*) (e.g., pile movement at depth *z* below the pile head level) at depth *z* have been offered as Equation (5):

$$\frac{dw(z)}{dz} = \frac{-4P(z)}{\pi D\_s^{-2} E\_p} \tag{5}$$

Equation (5) is different with *z* and *dP*(*z*)/*dz*, replaced by using Equations (1)–(4) as Equation (6):

$$\frac{d^2 w(z)}{dz^2} = \frac{4}{D\_s E\_p} f\_s(z) \tag{6}$$

All the mentioned equations according to force equilibrium and displacement correspondents are valid irrespective of pile type and soil grouping [7].

*Appl. Sci.* **2019** , *9*, 3060

with length of δ*z*.

Few more studies have been performed on longitude strain, which have measured the instrumented piles [8–13], transferring of load in rapid pile axial loading [14], static, dynamic, seismic, and cyclic lateral load of pile classifications [15–17], rigid and flexible pile behaviors in diverse soft soils [18–20], and skin friction resistance measurement in piles [21]. Some studies have been conducted based on numerical simulations, however, others have been performed on field monitoring (Ng and Sritharan [22]; Hung et al. [23]; Tafreshi et al. [24]; Mascarucci et al. [25]; and Lee et al. [26]). Ochiai et al. [27] propose a reliable designing model for bored-piles following in-situ tests by SPT. Poulos [28] has also introduced an appropriate design for piled rafts, comprising three stages: (1) assessing the feasibility of piled raft application accompanied by the required pile number, (2) evaluating where piles are needed and its general features, and (3) obtaining the optimum number, location, and configuration of a pile and computing the settlement distributions, bending moment and shear in the raft, and the pile loads and moments.

Sego et al. [29] studied the effect of an enlarged base on the total and end bearing resistance of a pile for use in ice-rich permafrost. Therefore, the total capacity of a bored-pile has mainly inclined by belled-pile usage. Shariatmadari et al. [30] studied the bearing capacity of driven piles in sands following SPT by applying 60 previous cases. Data included 43 full-scale, 17 dynamic tests, and static pile load testing analyzed by control and provisioning of wireless access points (CAPWAP). Note that, SPT data have been used close to the pile locations. Another model, as standard penetration test N(SPT-N), has been conducted and proved to have less scatter with higher accuracy. Zhang et al. [31] applied an elastic–plastic model showing the load–settlement relationship, and provided a simple method to analyze the behavior of a pile group and/or a single pile embedded in multilayered soils using two methods like an approach to enable a quick estimation of pile group settlement and/or a single pile embedded in multilayered soils providing time and cost saving. Ruan and Zuo [32] explained the relations between the ultimate vertical bearing capacity and SPT for the jacked pile. These case studies confirm the accuracy and reliability of the formula, mainly in silent pipe pile bearing capacity calculation. Assessing the friction capacity of driven piles in clay has been performed through the use of multivariate adaptive regression spline (MARS) by Samui [33], which led to D (diameter), L (length), undrained shear strength, and effective vertical stress as MARS input and friction capacity as output compared using artificial neural network (ANN). MARS, as a robust model, has been applied to predict the friction capacity in driven piles installed in clay.

Meanwhile, another research carrying the same approach has been conducted on cohesionless soil by Samui [34]. Chae et al. [35] investigated the uplift capacity of belled-piles in weathered sandstones of the Persian Gulf coast. Accordingly, few full-scale pullout load tests on belled tension piles in Abu Dhabi have also been performed. Comparing the results from the field, 3D finite element (FE), and theoretical methods have overestimated the ultimate pullout resistance of belled-pile without bell shape considerations [35]. Sakr [36] compared the results of high strain dynamic and static load tests of single helical-screw piles in cohesive soils. High strain dynamic pile load experiments have been performed on both driven steel open-pipe and helical piles. Case Pile Wave Analysis Program (CAPWAP) and full-scale static load test have confirmed high strain dynamic testing (HSDT) as a reliable tool for assessing static helical pile capacity. Zhang et al. [37] stated that in a typical design, the skin friction and the base resistance are majorly independent. Furthermore, the ultimate bearing capacity of a single pile has been composed in ultimate end resistance and confining skin friction.

Yao and Chen [38] promoted a flexible plastic solution for the uplift belled-pile. In the comparison of the provided solution outcome to the theoretical calculation one, the theoretical method results revealed that the elastic–plastic analytical solutions is a good method. Aksoy et al. [39] developed a new chart to estimate the friction angle between the pile and soil materials. Accordingly, in the current research, soil including different internal friction angles (φ) has been initially provided, then the skin friction angles (delta) of the mentioned soils with fiber-reinforced plastic (FRP) like a composite material, wood (pine), and steel (st37) have been defined by undertaking the interface shear test to provide a pile design diagram to determine skin friction angles of the soils and pile materials.

Wang et al. [40] investigated the controlling effectiveness and settlement behavior of two types of rigid pile structure embankments (PSE) constructed on collapsible loess soils beneath a high-speed railway. The results have shown that this type of PSE has mainly reduced embankment settlement so that embankments have to be maintained on collapsible loess. Meanwhile, pile spacing has significant efficiency in settlement reduction. Therefore, the current study has focused on stress–strain transferring throughout the instrumented bored-pile within the layered soil to measure the parameters of soil-pile interaction comprising Young elasticity module in concrete (*Ec*), average strain, and unit skin friction changed along with the pile.

#### **2. Material and Methods**

#### *2.1. Testing*

Maintained load test (MLT), is known as static loading test (i.e., the load will remain constant until the settlement ended to small values), has been followed by ASTM Standard (D1143/D1143M-07). An apparent distance between the reaction pile and test pile should not be less than five times to diameter (D) of the immense pile. Considering setup, the piles are loaded by applying hydraulic jacks toward the main beam operated by an electric pump. Then, the applied load is accurately calibrated by vibrating wire load cells (VWLC) (Table 1). Both instrumental piles are located in Jalan Ampang and Kuchai Lama, Kuala Lumpur, Malaysia.


**Table 1.** Instrumental bored pile load test summary.

#### *2.2. Loading Procedure*

Poulos [28] states that applied loading is crucial to the bending moment and differential settlements, but less critical to load-sharing or maximum settlement between the piles and raft. In bored-piles, time for loading test has been determined by the piles' concrete strength (Zhang et al. [37]). In addition, Tomlinson and Woodward [1] suggest that at the testing time, concrete should be in its seven days (at least) with at least doubled strength of applied stress, moreover, in the current pile testing, the load cycles have started 28 days after pile construction. A schematic view of the MLT is shown in Figure 2. Instrumented piles are tested by MLT per two loading-cycle calculated by calibrated VWLCs in every 10 min for one hour (Figure 3).

**Figure 2.** Schematic view of the maintained load test (MLT)

**Figure 3.** *Cont*.

**Figure 3.** Variation of the total applied static stress on pile versus time, (**a**) Bored pill (BP)# 1 and (**b**) BP# 2.

#### *2.3. Instruments Monitoring System*

The influence of geologically weak zones through multilayer site conditions has dramatically changed the designing parameters when the majority of these parameters change by depth. In other words, new global strain extensometer (GSE) sensors have recently measured the change of design parameters through the depth, so the pile top settlement has been monitored using: (1) four linear variation displacement transducers (LVDTs) mounted on the reference beams with their plungers vertically placed against glass plates fixed on the pile top, and (2) vertical scale rules attached to the pile top and sighted by a precise level instrument. Vertical scales have also been shown on the reference beams to check any movement in the loading test. Indeed, VWLC, global strain extensometer (GSE) sensors, VWSG, and LVDTs have been automatically logged across the use of Micro-10x data logger and Multi logger at 5-minute spaces for precise controlling during loading/unloading stages, adding that only accurate level readings have been taken manually. Figure 4 shows the cross section of BP# 1 (D = 1,800 mm) and BP# 2 (D = 1,000 mm) and the sensors' placement along the main steel rods.

**Figure 4.** Cross-section of BP# 1 (D = 1,800 mm) and BP# 2 (D = 1,000 mm).

Accordingly, the bored-pile has been tested using MLT through the reaction pile. All instruments using Micro-10 data logger and multilevel software have been automatically logged. The conventional instrumented method used VWSG and mechanical tell-tales. While by being attached to the steel cage of bored-pile, VWSG and mechanical tell-tales have permanently been embedded in the concrete. GSE is the second instrument applied in axial load measurement and settlement distribution in the bored-pile.

On static load testing, loaded pile deformation has produced a related moving between every two anchored intervals changed in strain gauge wire tension, in addition to a corresponding change in its resonant vibration-frequency measured by plucking GSE sensors/transducers through a signal cable to readout box/data logger to measure the frequency and display of the shortening and strain reading. Considering the installation set up, GSE has measured shortening and strains on all test pile sections in every load step of a static pile load test, so it has integrated the strains on a larger and more representative sample. Therefore, using a defined instrumental scheme, data derived from instrumented load testing have provided reliable information. The results of GSE have been compared to the conventional instrumentation bored pile results. Subsequently, regarding the test piles, including BP# 2, the Geokon VWSG and tell-tale extensometers have been installed internally to monitor the strain developing and shortening of pile behavior on the test. BP# 1, with the instrumentation of GSE, has been placed with seven levels depending on the pile length and vertical varieties of sub-soil cases in sonic logging tubes (Figure 5).

Due to calibration of the applied axial load and the measured average strain, a calibrated GSE sensor has been installed near the pile's head, in which there is no interaction between the soil friction and pile shaft. The GSE sensor measured the strain of other levels to determine the axial load transferring in all pile shaft sections, also to measure the loads contributed from the toe or/to end bearing resistance. VW Extensometer is installed on the anchored interval at eight levels (Figure 6). Pile deformation in loading has produced a related motion between every two anchored intervals producing strain gauge wire tension alteration in VW transducers, in addition to a correspondent variation in its frequent resonant vibration measured by plucking GSE sensors to the readout box/data logger. Therefore, the process measured the frequency displayed by the shortening reading and strain reading. VWSG for BP# 2 is also installed at five levels (A to E) in four numbers per level (Figure 6). The connected VWSG to the steel cage and electrical lead wires from the sensors coming to the top, have been illustrated in Figures 7 and 8.

The pile head displacement has also been analyzed using dial gauges and LVDTs, resulting in a reading with 0.01 mm accuracy mounted on stable reference beams and protected from direct sunlight and disturbance of personnel in the whole system. Settlement measuring through the use of proper leveling techniques has also been implemented as a useful backup to check the reference beams' movement. Vibrating wire load cells, strain gauges, retrievable extensometers, and LVDTs have been automatically logged by applying a Micro-10x Datalogger and Multilogger software at the 3-minute spacing for precise control during load/unload stages. Only accurate level readings are manually considered.

Briefly, the pile movement monitoring system, such as pile top and bottom settlements, has been monitored using the vertical scale rule fixed to the pile top observed by tunnel boring machine (TBM) for correction purposes. On the other hand, LVDTs and GSEs mounted to the reference beam accompanied by a plunger pressed vertically against a glass plate fixed to the pile's top. The vertical scales have been provided to monitor any movement in the load test. Indeed, strain gauges are manufactured in Geokon, USA. VWSG principles in strain measurement were presented due to the frequent natural vibration of taut wire restrained at both ends and varied with the square root of wire tension. Indeed, this change of tension in the wire has shown all the strain alterations in the structural steel member where the gauge is mounted. Moreover, all strain gauge has been mounted to two end blocks as arc-welded to the main reinforced bar at a proper level. The signal cables of the picked-up sensors fixed to the strain gauges have been tightly tied in the reinforced bars to the top of the piles terminated in a multiplexer box and observed using a Micro-10 Datalogger.

**Figure 5.** Various instrumental installments at different global strain extensometer sensors levels in BP# 1.

**Figure 6.** Vibrating wire strain gauges and tell-tales' extensometer arrangement in BP# 2.

**Figure 7.** VWSG attached to the steel cage.

**Figure 8.** The electrical lead wires from the sensors come to the top.

#### *2.4. Site Condition*

In the current paper, two full-scale maintained static load experiments on bored-piles are conducted. The first experiment (BP# 1) was performed at Cadangan Pembangunan, Lorong Stonor, Kuala Lumpur, Malaysia. The test pile was initially loaded up to 2 times to pile structural capacity, therefore, regarding BP# 1 with the structural capacity of 22,200 kN, the nominal diameter of 1,800 mm with a penetration depth of 36.95 m from the current piling platform level is RL 36.25 m. The pile was initially examined by up to 44,400 kN (2 x working load) in two loading cycles. The second test (BP# 2) was applied

at Utama Lodge, Jalan Senangria, Kuala Lumpur, Malaysia. The summary of soil types, besides the SPT-N values measured near the pile location, are illustrated in Table 2. Noting that the soil stratum was classified according to the unified soil classification system.



\* Note: SPT-N is the number of SPT driven into the ground (e.g., at the bottom of a borehole) by blows from a slide hammer with a mass of 63.5 kg.

The gauges were investigated prior and after installation, after cage placement in the borehole, and after concreting. The strain gauges' signal cables were reserved for testing approximately after 28 days, allowing for the concrete to achieve the design strength. Therefore, on test day, the strain gauges' cable was linked to the switch box connected to the data logger to ensure the functional sequences. Regarding the rod extensometer, galvanized iron (GI) pipes were tied to the main reinforced cage with steel wires at each terminating depth (Figure 9). A mild steel rod (10 mm) was inserted untill it touched the bottom of the pipe. In addition, a steel plate was welded to the rod's end for the plunger to sit on along the experiment (Figure 10).

**Figure 9.** global strain extensometer pipe for tell-tales extensometer is pre-installed at VWSG.

**Figure 10.** Existing platform in full-scale maintenance of the pile load test.

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

#### *3.1. Stress–Strain Variation in the Piles*

Stress distributions along the piles (BP# 1 and BP# 2) in two continuous loading/unloading cycles have been illustrated in Figures 11–14, showing that in a specified load, normal stress used in the pile's surface area is reduced by depth. Therefore, the reduction rate varied along with the pile—it is low in soft soil layers (having low SPT-N) but sharp in stiff layers (with high SPT-N). These results indicate great skin friction capacity of soil where rapid change between two different depths.

Tosini et al. [41] have declared that the forecast of deep foundation settlements in a layered soil profile is not always straightforward due to the problems in value defining of mechanical parameters affecting them. Distribution of average change and their cycles, in both strains (BP# 1 & BP# 2), in addition to the distribution of back-calculated modulus of concrete as *Ec* (kN/mm2) in both strains (BP# 1 & BP# 2) are shown in Figures 15 and 16. In BP# 1, *Ec* for the entire length remained constant at 25 kN/mm2 except at depths below 5 m within stiff sandy silt with little gravel when the higher load showed slightly lower *Ec* (Figure 15a). Therefore, the results obtained for *Ec* (in BP# 2) changed based on the depth and applied load. The lower applied stress through the first cycle provided minor changes in *Ec* distribution along the pile length. However, higher applied load amounts caused a higher change in *Ec* with depth (Figure 15). The highest and lowest *Ec* for the depth = 17.816 m (in BP# 2) were 38.2 and 33.9 kN/mm2, resepectively, for the applied stress of 1,750 kPa and 12,904 kPa. The reverse relation of *Ec* with applied stress occurred when higher applied stress on the pile cap showed lower *Ec* (Figure 16a). Consequently, the highest variation in *Ec* was observed within the silty and/or sandy silt layers, while the lowest change in *Ec* was measured at the weathered sandstone zone at depths above 20 m (Figure 16b).

**Figure 11.** Stress distribution measured along BP# 1 in (**a**) 1st cycle and (**b**) 2nd cycle.

**Figure 12.** Stress distribution measured along BP# 2 in (**a**) 1st cycle and (**b**) 2nd cycle.

**Figure 13.** Distribution of average change in strain measured along BP# 1 in (**a**) 1st cycle and (**b**) 2nd cycle.

**Figure 14.** Distribution of average change in strain measured along BP# 2 in (**a**) 1st cycle and (**b**) 2nd cycle.

**Figure 15.** Distribution of back-calculated elastic modulus of concrete, *Ec* (kN/mm2) for BP# 1; (**a**) 1st cycle and (**b**) 2nd cycle.

%DFNFDOFXODWHGPRGXOXVRIFRQFUHWH(FN1PP

**Figure 16.** Distribution of back-calculated elastic modulus of concrete, *Ec* (kN/mm2) for BP# 2; (**a**) 1st cycle and (**b**) 2nd cycle.

Average SPT-N and unit skin friction in BP# 1 & BP# 2 presented in Figure 17, indicate the range of ultimate skin friction with value change of SPT-N. Furthermore, in multilayer site conditions, the least SPT-N sum for a soil layer provided the least unit skin friction and vice versa. The outcome derived from GSE sensors was computed according to the displacement amounts recorded by the sensors. A higher alteration for the recorded axial force within two continuous levels provided a larger unit skin friction for a specified soil layer (Figure 17).

**Figure 17.** Average number of standard penetration test and the unit skin friction in BP# 1 and BP# 2.

#### *3.2. Pile Movement Monitoring*

Zhang et al. [42] state that the stress load-settlement curves reflect (1) the pile–soil interaction law, (2) the load transfer law, and (3) the pile load destruction mode. According to the previous explanations, the pile top and base settlements have been monitored for each load increment by applying the dial and strain gauges. Applied stress alteration versus total pile top and base settlement for two well-instrumented field tests (BP# 1, and BP# 2) in multilayered soils are presented in Figures 18 and 19, respectively. In the 1st cycle, the highest sighted pile top settlement at loading 22,418 kN was 9.60 mm. During unloading to zero, the pile was rebounded to a residual settlement of 0.36 mm. On the contrary, in the 2nd cycle, the maximum sighted pile top settlement at the peak load of 44,036 kN was 24.63 mm, so during unloading to zero, the pile was rebounded to a residual settlement of 5.34 mm. Similarly, residual settlements of 2.4 mm and 6.58 mm were recorded for the BP# 1 and BP# 2, respectively. According to Omer et al. [7], the variation of *fs*(*z*) with depth (*z*) is affected by different parameters including pile–soil properties, such as(1) pile–soil interface geometry and slip properties, (2) stress performance on the pile–soil interface, (3) pile installment technique, and (4) pile load method and ratio.

**Figure 18.** Different applied stresses versus total pile top settlement in (**a**) BP# 1 and (**b**) BP# 2.

(**b**)

**Figure 19.** Different applied stress versus total pile base settlement in (**a**) BP# 1 and (**b**) BP# 2.

Different total pile settlement versus time during BP# 1 and BP# 2 is presented in Figure 20. The settlement rate in the pile's head is almost linear in the loading steps. However, when unloading begins, the settlement rate of unloading steps significantly increased depending on the loading/unloading sequences, showing that the loading time for the second cycle of BP# 1 and BP# 2 were 400 and 600 min. The measured settlement on the pile's head was rebounded to permanent vertical deformation of 5.34 mm and 8.55 mm for both tests (BP# 1 and BP# 2) after unloading to zero. Correspondingly, the highest settlements of 24.63 mm and 19.54 mm were recorded for the pile's head vertical deformation.

**Figure 20.** Variation of the total pile settlement versus time during pile test in (**a**) BP# 1 and (**b**) BP# 2.

Some of the critical factors that were considered to be constant during pile load tests (e.g., BP# 1 and BP# 2.) were (i) piling technology, (ii) concrete maturity, (iii) the location of the groundwater, and most importantly, and (iv) the soil parameters (e.g., both mechanical and physical properties changes). Such a problem may affect *Ec* measurement during the pile servicing period. For instance, it is established that piling technology will influence soil–pile interactions [43–45]. Piling techniques

that lead to changes in soil properties can affect the axial force, load-displacement response, and tip resistance of each model pile. In most cases, to assess such influences, a large number of small-scale experimental works, or real-scale numerical and analytical studies are helpful.

On the one hand, the other factor that can influence the results of such investigations is concrete maturity [46,47]. Concrete maturity factor reveals the amount of concrete strength gain during the curing period, which is typically challenging to be taken into consideration as a separate variable on the full-scale experimental programme. It is important to note that the lack of such information may cause a significant change in the load-settlement behaviors of the pile during its working lifetime. Factors such as soil properties (e.g., shear strength parameters such as soil internal friction angle, cohesion, etc.) as well as groundwater levels, are primary terms that can remarkably alter the soil–pile responses to heavy external loadings. As an example, saturating the soil can cause soil shear strength reduction, which will influence the pile settlement as well as reducing the pile bearing capacities [48]. Another critical issue that can increase the complexity of soil–pile reactions, as well as load-settlement responses, as highlighted by Chisari et al. [49]. In the study provided by Chisari et al. [49], the influence of dynamic and static loading conditions (e.g., for identification of the primary material properties of a base-isolated bridge) are investigated. Their results showed that static identification is much less complicated compared to dynamic analysis. Although the current study covers the static load test, future work could evaluate the effect of a dynamic loading test in real-time monitoring of *Yc*.

#### **4. Conclusions**

The main objective of this study was to find a reliable estimation of the *Ec* in the installed bored pile. Two full-scale maintained a static loads test on instrumented bored-piles had been conducted in Kuala Lumpur, Malaysia to obtain a reliable range for ultimate skin friction with SPT-N value (i.e., blow counts) alteration. The effects of geologically weak zones through the layered soil ground conditions in crucial parameter-design changing such as elastic concrete modulus and strain–stress along the piles have also been researched. The details of the conclusion are as follows.


**Author Contributions:** Data collection and experimental works: H.M., B.K., H.N., Writing, discussion, analysis: H.M., A.S.A.R., and R.N., M.M.A.

**Funding:** This research has received funding from Ton Duc Thang University.

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

#### **List of Abbreviations**


#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Modeling of the Snowdrift in Cold Regions: Introduction and Evaluation of a New Approach**

**Mengmeng Liu 1,2, Qingwen Zhang 1,2,\*, Feng Fan 1,2 and Shizhao Shen 1,2**


Received: 17 July 2019; Accepted: 14 August 2019; Published: 17 August 2019

**Featured Application: A new approach for investigating snowdrifts around buildings and snow loads on building roofs based on a new facility named the "snow–wind combined experimental facility" is proposed in this paper. First, the method proposed in this paper is convenient and economical, and also the experimental results prove that the method is feasible and reliable. In this paper, only experiments on snowdrift around a cube were chosen for the verification for the sake of generality, but this method could also be adopted for the prediction of snowdrift around complex building environments and snow loads on various shapes of building roofs.**

**Abstract:** Unbalanced, or non-uniform, snow loads caused by snow drifting or sliding in cold regions with heavy snowfalls, can be a serious problem for the building industry. However, the methods for predicting snow distribution still need to be improved. Field observation is the most direct and reliable method to study snow distribution, but because the natural environment is uncontrollable and varies dramatically, sometimes conclusions may be confused under the influence of the many variables in the investigation. This paper proposes a snowing experiment approach using an outdoor snow–wind combined experiment facility for the study of snow distribution. The facility can produce a stable and controllable wind field and snowfall environment. Experiments which focused on snowdrift around a building were conducted during the winter to make an evaluation of the repeatability and reliability of the new approach. Finally, from the analysis of results, it was demonstrated that the experimental facility was stable and that the similarity criterion adopted for the snowing pattern was reliable. Especially, the minimum value of the friction speed ratio was suggested to ensure the test accuracy.

**Keywords:** snow–wind combined experiment facility; snowdrift; field observation; scale experiments; similarity criterion

#### **1. Introduction**

Since ancient times, buildings have been constructed to protect people from the natural environment. The environment dictates what type of building should be constructed, and snow loads are one of the dominant live loads that should be considered in the design of buildings in cold regions with a severe winter climate. Investigations show that snow loads are usually unbalanced, due to the action of the wind. And snowdrifts around buildings may cause problems for the vehicular traffic and pedestrians.

Currently, field observations, numerical simulations, and wind tunnel tests are the three main methods to estimate snow distribution, and much work has been undertaken to explore how the wind affects snow distribution.

Field observation is the most direct and reliable approach for obtaining information on snow distributions. Investigations of snowdrifts around an isolated building or group of buildings have been carried out in several previous studies [1,2]. Thiis and Gjessing [3] investigated the snowdrifts around three different model buildings in a valley in Spitsbergen, Norway, and the results showed that a change of roof shape could lead to large differences in snow distribution around the building. Beyers et al. [4] carried out research to investigate the snowdrift around a 2.0 <sup>×</sup> 2.0 <sup>×</sup> 2.0 m<sup>3</sup> cube at the SANAE IV research station, Antarctica, during the summer of January 2002 and verified the accuracy of a numerical simulation model using the results of field observations. Høibø [5,6] measured snow loads on nearly 200 building roofs under different wind and sun/shade environments from 1966 to 1986 in southern Norway. This method has strict requirements on the environment of the test site, such as wind velocity, temperature, humidity, and snowfall, but these important parameters in the observation were unstable and uncontrollable.

In 1991, Uematsu [7] first used numerical simulation to analyze snowdrifts, and now this method is widely used by researchers. Tominaga et al. [8] analyzed the snowdrift around an actual apartment building using the revised k-ε model. Beyers et al. [4] predicted transient snowdrift around a cubic structure using the standard k-ε model. And a large amount of numerical work has been done recently due to its significant advantage in efficiency and convenience [9–12]. However, only a few attempts have been made to apply CFD (Computational Fluid Dynamics) to roof loading so far [13]. Thiis et al. [14,15] predicted the snow distribution on a curved roof of a sports hall located in Oslo and compared the results with measurements. Some researchers are confused by the selection of the turbulence models, hence there is no approved and unified model.

Wind tunnel testing [16–18] is the most effective method for investigating snow distribution under the action of airflow. Kind [19,20] identified the important similarity criteria by analyzing the saltation process and discussed the necessary compromises in the modeling procedures. Delpech et al. [21] used artificial snow to simulate snow drifting around buildings for the Antarctic Concordia research station. Okaze et al. [22] investigated the development of drifting snow in a boundary layer wind tunnel. LÜ et al. [23] conducted a series of experiments in a wind tunnel to investigate the motion of natural snow and found that the threshold wind speeds for fresh and old snow were 6.2 m/s and 6.8 m/s (at 10 m height). Zhou et al. [24,25] investigated the snow distribution on a stepped flat roof using different granular materials in an open straight-circuit wind tunnel and also conducted a series of experiments on snow loads on flat roofs using high-density silica in the open test section of a boundary-layer wind tunnel. However, the snowdrift simulations in a boundary-layer wind tunnel struggle with scaling issues and similarity criteria. Some previous experiments on snowdrifts using wind tunnels have used bran, or other particles, as a substitute for natural snow, but these cannot completely replicate the distribution of natural snow. Tests have been carried out in advanced climate wind tunnels [26], where the influencing factors that may affect the results, such as the wind speed, temperature, and snowfall rate, are fully controllable. However, they are very expensive and not readily available for most researchers.

Snowdrift loads around buildings or on the roofs, which have a great impact on functionality and security of buildings, require particular attention. However, for the existing research methods, numerical simulation still needs to be improved with data from field observations and wind tunnel tests. Field observations are rarely available when considering security and technical issues, even though they should provide the most reliable data. Wind tunnel tests are the main source of data, but the expense of experiments in climate wind tunnels is so high that few experiments are reported, and it has restricted further studies. What is more, the wind speed and direction of field observations vary dramatically, so some conclusions may be confused when the natural environment is problematic. Hence, this paper proposes a new approach to investigate the snowdrift loads with a facility in which the wind speed and direction are controllable and a snowing environment, just like natural snowfall, is available.

The purpose of this facility is to provide a new approach to study snow distribution under the action of wind. In comparison with traditional wind tunnel tests, there will be a 90% or more decrease in the construction cost because the facility need not be equipped with a refrigerating system. Besides that, the operating expense of the facility is much lower than climate wind tunnels.

Transportation of snow particles can be classified into three processes, i.e., creep, saltation, and suspension. Saltation, which is widely accepted as contributing most to the total transfer volume, is a process in which snow particles move with repeated leaping up, or jumping, on the snow surface.

The height of saltation is usually below 0.1 m. Generally, the length of saltation *L* [27] is suggested to be

$$L \cong 10 \text{ h.}\tag{1}$$

Yet, to ensure that the saltation process could be accurately reproduced on the roofs or around the building, it is suggested that the characteristic length of the prototype model is larger than the length of saltation *L*.

In this study, the design parameters of the facility, experiment method, and the adopted similarity criteria are introduced and discussed. Finally, a prototype model with external dimensions of 1.0 <sup>×</sup> 1.0 <sup>×</sup> 1.0 m3 was selected, taking both economy and theory into consideration, to study the surrounding snow distribution and evaluate the repeatability and reliability of the new approach.

#### **2. Design of the Snow–Wind Combined Experiment Facility**

#### *2.1. Brief Introduction of the Facility*

In the study of snow distribution, numerical simulations and wind tunnel tests should be verified by field observations. For this purpose, a snow–wind combined experiment facility as shown in Figure 1a has been constructed at the Harbin Institute of Technology (HIT) to produce stable wind fields to investigate snow distribution [28] The climatic characteristics of the winter months in Harbin provide convenient conditions for the experiments. Harbin has a long winter, which can last for almost six months, with snowy conditions and a temperature range from −30 ◦C to 0 ◦C. The location of the facility is shown in Figure 1b. Because the facility was to be sited outdoors, a narrow strip between two buildings was chosen to reduce the interference of the external environment. The dimension of available ground is 12 m × 20 m, and the heights of the surrounding buildings are all greater than 10 m.

**Figure 1.** (**a**) Photo of the snow–wind combined experiment facility. (**b**) Location of the facility.

#### *2.2. Introduction of the Main Component*

The facility, as shown in Figure 2, consists of a power section, a flow conditioning section, an experimental section, and a snowfall simulator.

**Figure 2.** Schematic view of the snow–wind combined experiment facility.

The power station is the 3 × 2 fan-matrix, and the area affected by the model during a snowdrift is thought to be about three times larger than the characteristic length of the model [27]. Therefore, the cross-section of the facility was designed to be 4.5 <sup>×</sup> 3.0 m2. Considering the fact that there are usually light winds when it is snowing, the facility was designed to be able to produce a stable wind field for a specific duration with a range of speeds from 0.5 to 11.5 m/s. Wind pressure is the main indicator of the output of the fan-matrix and is defined as

$$P = \frac{\rho \text{k}lI\_0}{2}'\tag{2}$$

where *P* is the wind pressure; ρ is the density of air, and ρ = 1.395 kg/m3 when the temperature is −20 ◦C; k is the resistance coefficient of pressure of the flow conditioning section, which is defined as the ratio of the outlet area of the facility to the total outlet area of the six airflow fans, and k = 1.734; *U*<sup>0</sup> is the wind speed at the outlet or immediate output from the fan matrix, and *Umax*, which is measured at the section 4 m away from the outlet, is regarded as the available maximal wind speed of the facility and *Umax* = 11.5 m/s for design purposes.

Because the facility is open to the natural environment, a numerical simulation was undertaken to explore the relationship between *U*<sup>0</sup> and *U*, where *U* is the wind speed at the experiment section. Based on the environment of the experimental site, a calculation model of the facility, with a computational domain of 10 <sup>×</sup> 10 <sup>×</sup> 20 m3, was formulated. The result showed that when the initial velocity of the wind export of the facility was 15 m/s, the uniform wind speed in the range of 4–20 m was about 11.5 m/s, hence, the value of *U*<sup>0</sup> is approximately equal to 15 m/s.

Accordingly, the minimum wind pressure is 272.1 Pa by calculation, and a fan with 286 Pa wind pressure was chosen. The maximum speed measured in the field was 11.5 m/s, which reasonably meets the requirements of the experiments.

The flow conditioning section, which stabilizes and smooths the airflow produced by the fan, is a combination of the draft tube, the honeycomb, and the damping net. Two pieces of damping net were stuck on both import and export sections of the honeycomb, and the draft tube served as connection between the power section and the honeycomb. The flow conditioning section was designed to be short to reduce the loss of wind power because of the open experimental environment. Considering the square pattern of the cross-section and the design experience from the wind erosion tunnel designed by Saxton [29], a square cross-section with the dimension of 30 <sup>×</sup> 30 mm2 was chosen as the single unit of the honeycomb with a length of 300 mm. The experimental section is a platform with a circular dial, which is used to fix the model at the required angle.

Snowfall suitable for the experiments was not always available, so to accumulate a large amount of data on snow distribution, it was necessary to find a way to simulate snowfall using natural or artificial snow particles.

In this work, a snowfall simulator situated 1 m away from the wind export was designed to simulate snowfall during the experiments. Before the experiments started, snow, which had been collected in heat-insulation boxes, was placed in a metal box with wire gauze on the bottom. Then, shaking the metal box with an eccentric motor to generate friction between the snow and wire gauze, made snow particles fall from the mesh into the stable wind field caused by the fan-matrix so that a snowfall environment was formed in the test section. Different snowfall conditions are available using different kinds of wire gauze with varied vibration frequency. There are three different kinds of the aperture of wire gauze, 2 mm, 3 mm, and 5 mm, and the vibration frequency ranges from 1 to 5 Hz. The metal box, which is 5 m wide, can be raised and lowered freely in the range of 0.5–3 m.

#### *2.3. A Brief Introduction to the Experimental Procedure*

Figure 3 shows a sketch of the experimental procedure. The facility can provide a controlled snowfall environment in the test section so that experiments examining snow distribution around the building or on roofs can be undertaken.

**Figure 3.** A sketch of the experimental procedure.

When an experiment is going to be conducted, there are four steps: First of all, enough snow should be put into the box of the snow simulator, next the model should be fixed according to the experiment design along with an anemograph for monitoring the wind speed at the reference point; then the fan-matrix is started and adjusted until the wind speed reaches a given value, and then the snow simulator can be turned on for the required duration; finally, the snow load is measured after all of the equipment has been turned off.

#### **3. Verification of the New Approach**

#### *3.1. Field Measurements of Wind Field Produced by the Facility*

Obviously, the wind field has a significant impact on the results of the experiments, so field measurements of the actual wind field were also made to investigate it and determine the best location of the experimental platform.

Figure 4 shows the method for measuring the wind field. Because the cross-section is symmetrical, the measuring points are arranged symmetrically in the full section. The wind speed was measured by a hot-wire Thermo-Anemometer and each measuring point was measured at least 30 s with a sampling frequency of 1 Hz. The vertical and horizontal measurement points in each cross-section formed a matrix of 13 × 9. Before the measurements, a design speed U was defined. If the measured wind speeds at the 1.5 m point on the first line of the five sections were all between 0.9~1.1 U, this would

be recognized as the design speed U. Five cross-sections were measured with a design speed (U) of 3.0 m/s and the distance between them and the export was 1 m, 3 m, 4 m, 5 m, and 7 m.

**Figure 4.** A sketch of the measuring method.

Figure 5a shows that the mean wind speed on the first line of five cross-sections varies with height; the Y-axis represents the height of the cross-section and the X-axis represents the ratio between *U\** (the measured wind speed) and *U*. Figure 5b shows the turbulence intensity of wind on each measuring point. It is apparent that the turbulence intensity is higher in the section before 3 m and the wind is more uniform in the range 3–7 m; this was also shown in the numerical simulation. However, the wind speed decreased faster with increasing distance from the export, which is slightly different from the numerical simulation. To ensure a stable and uniform wind field during the experiments, the effective range of the experimental section was 3–7 m away from the wind export.

**Figure 5.** Wind field on the first line of the five sections: (**a**) mean wind and (**b**) turbulence intensity.

Figure 6a shows how the mean wind speed varies with width in the cross-section 4 m away from the export, and Figure 6b shows the turbulence intensity of wind on each measuring point. And from the results, it is apparent that the wind speed decreased rapidly close to the edge of the cross-section. To ensure the veracity of the experiments, the effective width of the experimental section was 4 m.

From a comparison of the measured wind speed with the normalized wind speed in the atmospheric boundary layer, the discrepancy above 1.5 m is non-negligible, but it is relatively acceptable for the given experiment. The turbulence intensity was small except at some points close to the wind export and at the edges of the cross-section, which indicates that the quality of the flow field in the open test section was reasonable and acceptable.

During the experiments, temperature, the natural wind field, and humidity were measured on a PC-4 automatic weather station, which is located on the side of the facility, to rectify the experimental wind field. This has three anemometers, at heights of 0.5 m, 1.5 m, and 2.0 m, and the sampling interval was 1 min.

**Figure 6.** Wind field on the section which is 4 m away from the outlet: (**a**) mean wind and (**b**) turbulence intensity.

#### *3.2. Physical Properties of Snow*

Properties of snow particles have a great influence on the result of the test [23] so that the physical properties of the stored natural snow particles were measured before the experiments. The diameter and shape of the snow particles were observed by an optical microscope and the angle of repose was decided from accumulated experimental observations. The terminal velocity of the snow particles was decided by measuring the time that the snow particles freely fell from 2 m to the ground level. A comparison of the physical parameters of the snow particles is given in Table 1.


**Table 1.** Comparison of the physical parameters of the snow particles.

#### *3.3. Experiments for the Verification of Repeatability*

To confirm the repeatability of the new approach, 20 repetitive experiments using stored snow on the snowdrift surroundinga1m cube were undertaken. Environmental conditions, namely temperature, humidity, and natural wind field were measured by an automatic weather station that was sited nearby the test section and results are shown in Table 2. The reference point for the experimental wind speed was set at 1 m in front of the model and 0.5 m above the ground. The snowing rate was defined as the accumulated snow depth on the ground in unit time and it was determined by the kind of wire gauze and vibration frequency. For the 20 experiments, the snowing rate was controlled to be the same. Figure 7 shows some photos of the experiment results.

Figure 8 gives the distribution of average snow depth and the standard deviation of each measured point along lateral and streamwise lines crossing the model. It is clear that the distribution of snow depths in the 20 experiments is relatively stable. The Pearson Correlation Coefficient for any two experimental results along the measuring lines shown in Figure 8 was above 0.90, which indicates that there is a strong correlation between the results from the 20 experiments.


**Table 2.** Environmental conditions of the 20 repetitive experiments.

**Figure 8.** *Cont.*

**Figure 8.** Results of the 20 experiments: (**a**) streamwise direction (center section) and (**b**) lateral section (center section).

From the above analysis, it can be concluded that the results of the 20 experiments, with the same test conditions, show good correspondence with each other. This indicates that the method provides good stability and repeatability.

#### *3.4. Experiments for the Verification of Reliability*

#### 3.4.1. A Field Observation on the 1 m Cube Prototype Model

To confirm the reliability of the new approach, the comparison of field observations and experiments should be done. So, the field observation on a 1 m cube prototype model was conducted on the HIT campus on 11 September 2015 from 0:00 to 6:00. To obtain a reasonably stable wind direction, the test site was located on a narrow strip between two buildings, as shown in Figure 1b. Figure 9 shows the method of observation; the wind field and environment conditions were measured by a PC-4 automatic weather station with a sampling interval of 1 min and the box-type snow flux trapper draws on the design experience summarized by Kimura [31].

**Figure 9.** Abridged general view of field observation.

The weather station and snow flux instrument were all mounted on the centerline of the model. The wind field during the observation is given in Figure 10. The average snowfall was 29 mm, which was measured on an open and flat ground 20 m away from the model. Figure 11 shows the photo of the snow distribution of the field observation.

**Figure 10.** Wind field during the observation.

**Figure 11.** Photo of snow distribution for the field observation.

#### 3.4.2. Test Similarity Criteria Based on Snowfall Mode

A reliable scale in experimental modeling depends on the right similarity criterion. In order to reproduce the non-uniform snowdrift phenomena caused by structures with scale models, the widely accepted fact is that it is necessary to meet geometric, kinematic, and dynamic similitude requirements.

Since the beginning of snowdrift research" many different similarity criteria have been proposed, validated, and applied by some famous scholars, such as Strom [32], Odar [33], Calkins [34], Kind [19], Iversen [16], Anno [35], Isyumov [18], Naaim [36], Delpech [21], Beyers [2], and so on.

However, just as it is not possible to satisfy the Froude numbers and Reynolds numbers simultaneously [37], it has been proved that some of the similarity criteria reveal incompatibility, therefore only the most important and widely accepted similarity parameters are selected and discussed.

In the present study, considering the fact that the snow particles falling into the snow bed are provided by the sowing pattern, the experimental modeling is mainly focused on the reliable simulation of the Froude number and drifting volume. The original form of the Froude number is defined as

$$\frac{\mathcal{U}^2}{\mathcal{g}^{L'}}\tag{3}$$

where *U* is the reference wind velocity, *g* is gravitational acceleration (m/s2), and *L* is the reference length (m). But considering the discrepancy of particle density between the prototype and the test, Odar [33] and Calkins [34] proposed the densimetric Froude number:

$$\frac{\rho}{\rho\_P} \frac{\mathcal{U}^2}{\mathcal{g}^L} \tag{4}$$

where ρ is the air density, and ρ*<sup>p</sup>* is the particle density. To confirm the reliable simulation of drifting volume, Anno [35] suggested that the test should satisfy the following time scaling parameter:

$$\frac{TQ\eta}{\rho\_P L^2} \tag{5}$$

where *T* is experiment time (s), *Q* is the transport rate of snow (kg/m·s), and η is the snow collection efficiency. Except for the above similarity, a number of similarity requirements, like the flow field, ejection process, particle trajectory, and deposition pattern, must be satisfied in the experiments on snow drifting on roofs to ensure a reliable simulation.

First, the near ground turbulent wind velocity profile for stable atmospheric conditions is

$$
\Delta U = \frac{\mathfrak{u}\_\star}{\kappa} \ln(\frac{z}{z\_0}),
\tag{6}
$$

where κ is the Von-Karman constant, *u*<sup>∗</sup> is the friction velocity, *z* is the height above the surface, and *z*<sup>0</sup> is the roughness height parameter. Geometric similarity should be enforced through the relationship

$$D\_p/L.\tag{7}$$

Except for the densimetric Froude number expressed as Equation (4), the following Froude number based on the threshold friction velocity should also be satisfied to simulate the shear stress of particles near the ground [16]:

$$\frac{\rho}{\rho\_p - \rho} \frac{u\_{\ast \mathbf{f}}^2}{\mathcal{g}D\_p} \tag{8}$$

where *DP* is the particle diameter, *<sup>u</sup>*∗*<sup>t</sup>* is the threshold friction velocity, and <sup>υ</sup> is the kinematic viscosity.

To ensure a similar movement trajectory of a particle, Kind [24] indicated that the following equations must be satisfied:

$$\frac{\rho}{\rho\_p - \rho} \frac{\mu\_{\ast t}^2}{\mathcal{g}L} \tag{9}$$

$$w\_f/\mathcal{U}\_\prime \tag{10}$$

where *wf* is the settling velocity of snow particles (m/s). The deposition or erosion mechanisms could be modeled satisfactorily through the similarity of the particle ejection process [19,38]:

$$\mathfrak{u}\*/\mathfrak{u}\*.\tag{11}$$

In the present study, because the simulation of the Froude number is considered to be more important, the limit of the Reynold number is released. However, to ensure that the inertial forces are dominating the flow, Kind [19] suggested that the following relationship should be satisfied:

$$\frac{\mu\_{\rm \ast f}^3}{2gv} > 30.\tag{12}$$

In summary, the parameter (4) and (5) were adopted for determining the wind velocity and experiment time, and several similarity parameters were also chosen for evaluating the results when using a sowing pattern to study snowdrift. In order to verify the reliability of this method, four scale models were selected for the verification test. The experimental conditions are listed in Table 3 and the prototype is the field observation presented above. The reference point for the experimental wind speed was set at 1 m in front of the model along the central line and the height above the ground was equal to the side length of the model. The major similarity parameters for the prototype and scale model are given in Table 4.

Snow flux in the vertical plane was measured by additional experiments with only the snow flux trapper on the testing ground. The snow load coefficient is the ratio of the depth of snow on the roof to the average snowfall. For the experiments, the average snowfall on the ground was determined by an additional experiment conducted for the same duration but without any obstacles in the test section. After the experiments, the snowdrift geometry was measured with the laser total station for about 200~300 points around the model and the measured snow depth data contributed to the contour maps of the snowdrift geometry.


**Table 3.** Environmental conditions of the scale experiments.


**Table 4.** Similarity parameters of prototype and scale models.

#### 3.4.3. Results and Discussion

There are two key factors to evaluate the accuracy of snowdrift given by the experiments. One is the similarity of the shape of snowdrift; another is the similarity of the snow distribution coefficient. The distribution coefficient is the ratio of measured snow depth and average snowfall at any measuring point.

Figure 12 shows the contour maps of the snowdrift geometry around the model of field observation and scale experiments. The photos of scale experiments are shown in Figure 13. All of the contour maps showed a similar shape which was formed like a horseshoe, in which snow mainly accumulated at the front of the windward side and beginning at the two corners of the windward side, there were two distinct erosion areas near the lateral sides along the downwind direction. From the comparison, if just focusing on the overall distribution shape of the snowdrift, it seemed that the result when the scale was 1/10 was more similar to the field observation and the main erosion areas (near the lateral sides) were larger than field observations for the other scale experiments. The basis for explaining this phenomenon needs an objective and appropriate understanding of the result of the field observation. The main erosion is more distinct in the experiment due to the stable wind speed. However, sometimes there were slight or no wind during the field observation, and snowfall contributes more to decrease the erosion effect in this condition, which would decrease the non-uniform snow distribution to some extent.

However, considering that two distinct snowdrifts were formed in front of the windward side, the results of the 1/1, 1/2, and 1/4 scale experiments had a more similar characteristic of snow distribution. When the scale was 1/10, the snow only accumulated at the feet of the windward side, indicating that wind velocity was excessively reduced; the near to the ground reverse flow in front of the windward side was weak and could not transport snow particles away from the windward side, so that only one distinct snowdrift was formed.

**Figure 12.** Contour maps of the snowdrift geometry around the model: (**a**) field observation; (**b**) experiment with a 1 m Cube model; (**c**) experiment with a 0.5 m Cube model; (**d**) experiment with a 0.25 m Cube model; and (**e**) experiment with a 0.1 m Cube model.

Considering that the effect of fluctuant wind velocity in the field observation could not be ignored, the 1/1 model may be treated as the self-test prototype. In this situation, according to the comparison of results, the 1/2 and 1/4 scale experiment results would have relatively acceptable accuracy of the reproduction of the snowdrift shape of the field observation, but the 1/10 scale experiment result would have a lower precision reproduction.

In order to compare the snow distribution in a more precise method, Figure 14 compares the snow distributions along lateral and streamwise lines crossing the model.

**Figure 13.** Photos of scale experiments.

**Figure 14.** Comparison between normalized snow depths obtained from field observations and scale experiments: (**a**) streamwise direction (center section) and (**b**) lateral direction (center section).

In Figure 14a, it could be seen that the snow distribution of the field observation has obvious characteristics: Near the windward side, the snow formed a compressed "N" shape, and near the leeward side, there formed a compressed "U" shape. When the scale was 1/1 and 1/2, the distribution results showed a good agreement with the field observation, even though the position where the extreme value occurred moved towards to the windward side as the scale ratio decreased. But, as the scale decreased, especially when the scale was 1/10, the snow tended to only accumulate at the foot of the windward side, and the snowdrift away from the windward side disappeared. As for the snow distribution near the leeward side, when the scale was 1/1 and 1/2, the experiments could precisely reproduce the snow distribution, even though the overall value of the coefficient decreased. When the scale was 1/4 and 1/10, the results gave a lower estimate of the snow accumulation, but for the snow distribution away from the leeward side, the results showed better agreement with the field observation both in value and shape.

Snow distribution along the lateral direction is shown in Figure 14b. Except when the scale was 1/10, the snow distribution showed a U-shaped distribution for both field observation and scale experiments. For the 1/10 scale experiment, a similar distribution was found, but the location where the minimum value occurred was much closer to the lateral side. This phenomenon indicates that the erosion effect was excessively decreased.

In summary, considering both the shape of the snowdrift and snow distribution coefficient, except when the scale was 1/10, the results of scale experiments reproduce the field observation accurately, suggesting that the method is reliable.

Especially, it is known that the relationship between friction velocity *u*<sup>∗</sup> and threshold friction velocity *u*∗*<sup>t</sup>* is of great significance in the study of the snowdrift. From the analysis of the similarity parameters in Table 4, it could be derived that a too-small value of the parameter *u*∗/*u*∗*<sup>t</sup>* may adversely affect the accuracy of reproduction. Therefore, it is recommended to limit the minimum value of the parameter *u*∗/*u*∗*<sup>t</sup>* for scale experiments.

Figure 15 shows the Pearson correlation coefficients of the distribution coefficient between the prototype and scale experiments along the measured line defined as the front, rear, left, and right.

**Figure 15.** Pearson correlation coefficients of the distribution coefficient between the prototype and scale experiments.

It could be concluded from Figure 15 that the results are of poor quality when the value of the parameter *u*∗/*u*∗*<sup>t</sup>* is lower. Assuming that the experiment results are accepted when the Pearson correlation coefficient is higher than 0.6, the lower limit values of the parameter *u*∗/*u*∗*<sup>t</sup>* are 0.353 and 0.287, respectively, when the field observation or 1/1 scale experiment is chosen as the prototype. At present, the result of the field observation is affected by the adverse effects of the unstable wind field and the quantitative assessment of the effect is still unclear. So, when giving a suggestion, the analysis based on treating the 1/1 scale experiment as a prototype should be more respected and, finally, 0.3 is suggested to be the lower limit value of parameter *u*∗/*u*∗*<sup>t</sup>* in the present study.

#### **4. Conclusions**

In the present study, a new approach based on a snow–wind combined experiment facility, for investigating snow distribution, is introduced, including its concept, design, operation, and verification. According to the above work, the following conclusions could be drawn:


In this paper, only experiments on snowdrift around a cube were chosen for the verification for the sake of generality, but this method could also be adopted for the prediction of snowdrift around complex building environments and snow loads on various shapes of building roofs. It could provide a basis for site selection planning for buildings and the safety design of roof structures in snowy areas.

**Author Contributions:** Conceptualization, Q.Z.; data curation, M.L.; formal analysis, M.L. and Q.Z.; funding acquisition, F.F.; investigation, M.L.; methodology, M.L. and Q.Z.; project administration, F.F.; resources, F.F.; software, M.L.; supervision, Q.Z. and S.S.; validation, M.L., Q.Z., F.F. and S.S.; visualization, Q.Z.; writing—original draft, M.L.; writing—review and editing, Q.Z.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number 51508133, 51478147, and the National Outstanding Youth Science Fund Project of the National Natural Science Foundation of China, grant number 51525802.

**Acknowledgments:** This work was financially supported by the National Natural Science Foundation of China (Grant no. 51508133, 51478147), the National Outstanding Youth Science Fund Project of National Natural Science Foundation of China (51525802), and Foundation of Key Laboratory of Structures Dynamic Behavior and Control (Ministry of Education) in Harbin Institute of Technology (HITCE201704). The authors are grateful to the members of the Space Structures Research Center in Harbin Institute of Technology, for providing invaluable information and advice in this study.

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **Numerical Analysis and Dynamic Response of Optimized Composite Cross Elliptical Pressure Hull Subject to Non-Contact Underwater Blast Loading**

**Mahmoud Helal 1,2,3,4, Huinan Huang 1,2, Elsayed Fathallah 1,2,5,\*, Defu Wang 1,2,\*, Mohamed Mokbel ElShafey <sup>6</sup> and Mohamed A. E. M. Ali <sup>5</sup>**


Received: 1 August 2019; Accepted: 19 August 2019; Published: 23 August 2019

**Abstract:** Among the most important problems confronted by designers of submarines is to minimize the weight, increase the payload, and enhance the strength of pressure hull in order to sustain the hydrostatic pressure and underwater explosions (UNDEX). In this study, a Multiple Intersecting Cross Elliptical Pressure Hull (MICEPH) subjected to hydrostatic pressure was first optimized to increase the payload according to the design requirements. Thereafter, according to the optimum design results, a numerical analysis for the fluid structure interaction (FSI) phenomena and UNDEX were implemented using nonlinear finite element code ABAQUS/Explicit. The propagation of shock waves through the MICEPH was analyzed and the response modes (breathing, accordion and whipping) were discussed. Furthermore, the acceleration, displacement and failure index time histories at different locations were presented. The results showed that the greatest acceleration occurred in the athwart direction, followed by the vertical and longitudinal directions. Additionally, the first bubble pulse has a major effect on athwart acceleration. Moreover, the analysis can be effectively used to predict and calculate the failure indices of pressure hull. Additionally, it provides an efficient method that reasonably captures the dynamic response of a pressure hull subjected to UNDEX.

**Keywords:** underwater explosion; composite pressure hull; whipping; breathing; failure index

#### **1. Introduction**

Significant research work has been presented in order to simulate the behavior of underwater vehicles under severe loading conditions. For instance, Reddy [1] illustrated the effect of shock pressure loading on a ring-stiffened submersible hull using finite element analysis. The failure analysis indicated that fibers failed in tension while matrix failed in shear when the explosion charge exceeded 25 kg TNT. Furthermore, Jen [2] worked on minimizing the weight of pressure hull by enhancing the pressure hull strength taking into consideration both hydrostatical pressure and underwater explosion. Additionally, the results indicated that the dynamic motion of the pressure hull has an accordion mode, a whipping mode and a breathing mode. Also, Cho et al. [3] derived an empirical formula for predicting the collapse strength of composite cylindrical structures under hydrostatic pressure as a function of important design parameters such as the geometric dimensions and the layered angle. In addition, Chen et al. [4] experimentally investigated the dynamic performance of ship body coated

by rubber. The results demonstrated that the coating was not effective at reducing the low-frequency whipping motion excited by the bubble pulse. However, it was able to moderate the high-frequency response excited by shock wave. Moreover, Ramajeyathilagam et al. [5] numerically and experimentally investigated the effect of underwater explosions (UNDEX) on rectangular plates. It was revealed that the underwater explosion failures can be predicted using the strain rate effects. Also, Liu et al. [6] illustrated the global responses of ship subjected to UNDEX. The results demonstrated that the UNDEX waves changed the added masses on the ship and effectively affected the global responses of its body. On the other hand, Liang et al. [7] examined the transient dynamic responses of submarine pressure hull exposed to hydrostatic pressure and shock loading. It was concluded that the collapse depth (maximum diving) of the submarine pressure hull was about 700 m. In addition, they observed that the loading condition depends not only on the hydrostatic pressure but also on the shock loading. Likewise, Kwon and Fox [8] studied experimentally and numerically the dynamic response of a cylinder subjected to UNDEX. The results illustrated that the largest strains on a cylinder subjected to a far-field side-on UNDEX occurred near the two ends on the near side of the cylinder to a far field explosive charge. Additionally, the damage occurred at the center on the opposite side of the cylinder. Similarly, McCoy and Sun [9] combined FSI code and finite element modeling techniques to investigate the dynamic response of a thick-section hollow composite cylinder. The results showed that the fluid-structure coupling has a significant effect on stress distributions within the structure.

Furthermore, Shin and Hooker [10] predicted numerically the damage response of submerged imperfect cylindrical structures exposed to UNDEX. Based on these results, the introduction of initial imperfections greatly affected the response of the cylinder when compared with the response of a perfect cylinder. On the other hand, Qiankun and Gangyi [11] demonstrated the shock response of a ship section to non-contact UNDEX using the finite element software package ABAQUS. It was revealed that the fluid thickness and size of fluid mesh effectively affects and improve the modeling accuracy. Additionally, Adamczyk and Cichocki [12,13] performed a numerical study to obtain the shock response of an underwater hybrid structure subjected to UNDEX. Zhao et al. [14] predicted the damage features of RC slabs subjected to air and UNDEX. The shock wave propagation and damage mechanisms from contact explosions in air and water were compared. The dynamic response of the RC slab is highly localized in the air contact explosion. Furthermore, the crater failure is observed at the top surface of the RC slab due to the direct impact of the air contact blast loading and the spalling failure occurs at the bottom of the RC slab. On the other hand, when the RC slab subjected to underwater contact explosion the top surface of the reinforced concrete RC slab almost completely destroyed. Therefore, underwater contact explosion can cause significantly more damage to the RC slab than the same amount of explosive in air. Rajendran and Narasimhan [15] investigate the damage of clamped circular plates subjected to contact UNDEX. The deformation contours were a spherical viewing the maximum absorption of energy for the depth of bulge attained. Also, Jacinto et al. [16] applied a linear dynamic analysis of plate models under explosions. The element size has a great effect on the results. Likewise, Kumar et al. [17] experimentally studied the blast effect on carbon composite panels the results showed that. There were two types of dominant failure mechanisms observed, fiber breakage and inter-layer delamination. Furthermore, Guo et al. [18] presented a new shock factor of twin hull water plane (catamaran) subjected to UNDEX. The shock factor parameter is used to describe the response of ships exposed to this loading condition based on shock wave energy. Similarly, Wang et al. [19] proposed a new method involving an analytical technique connected to elastic dynamic response of laminated plates exposed to UNDEX. The method was validated by comparing its results with those achieved by a semi analytical method and the experiment results. Zhang et al. [20] predicted the dynamic bending moment of a UNDEX bubble acting on a hull. The predicted numerical results showed that UNDEX bubble propagated a longitudinal bending which caused sagging and hogging damage for the ship. Furthermore, Gong and Khoo [21] presented a transient response of an UNDEX bubble on a glass/epoxy composite deep-submersible pressure hull. It was observed that the UNDEX bubble produced a huge deformation around the stand-off point in the pressure hull immediately after the collapse of the bubble, and the minimum volume was observed beneath the composite hull. On the other hand, Wang et al. [22] investigated numerically and experimentally the failure mode and dynamic response of a ship structure subjected to shock wave and bubble pulse. Jun et al. [23] investigated the impact environment characteristics of floating shock platform subject to UNDEX. Young and Leonard [24] modeled and simulate a surface ship shock to UNDEX, the results showed that the cavitation effect must be taken into account in the ship shock simulation, and that cavitation volume must be large enough. Gannon [25] investigated the response of a submerged stiffened cylinder to UNDEX using a coupled Eulerian Lagrangian model and experimental approaches.

This present study developing a procedure and describes a numerical modelling methodology for calculating the dynamic response of optimized Multiple Intersecting Cross Elliptical Pressure Hull (MICEPH) exposed to non-contact UNDEX. First, a submarine with pressure hull in the form of MICEPH is optimized using non-linear finite element analysis software ANSYS. Thereafter, according to optimization, the finite element model is built using non-linear finite element code ABAQUS/Explicit to examine the dynamic response of the pressure hull exposed to non-contact UNDEX.

#### **2. Analysis of Underwater Shock Loading and Bubble Pulse**

The most important element of any submersible body is the pressure hull. It contributes about one-fourth to one-half of the total underwater vehicle weight. Figure 1 presents various pressure hull configurations used in submersible bodies [26–28]. Exposing the hull structure to shock wave and bubble pulsation upon UNDEX events leads to great damage to the hull structure [2]. Figure 2 presents the different configurations that occur during UNDEX events [29,30]. The loading mechanisms resulting from UNDEX include incident wave, free surface reflection, shockwave, bottom reflection wave, gas bubble oscillation, bubble-pulse loading and bulk cavitation [31]. The compressive load on the structure of the hull is increased due the reflections from the bottom of the ocean which propagates from the shock wave, while the reflection of the shock wave from the ocean's free surface causes a reduction in the pressure [32,33]. Equation (1) presents the evaluation technique of pressure time history (*Pin*(*t*)) for the pressure profile, as follows [2,34–38]:

$$P\_{\rm in}(t) = P\_{\rm max} e^{-(\frac{t-t\_1}{\Theta})} \quad \text{ (Mpa)} \\ \text{?} \quad t \ge t\_1 \tag{1}$$

where (*t*) denotes the time elapsed after detonation of charge (ms), (*Pmax*) denotes the peak pressure (Mpa), (*t*1) denotes the arrival time of shock wave to the target after the detonation of the charge (ms), and (θ) denotes the time decay constant that describes the exponential decay (ms). The peak pressure and decay constant depends on the size of the explosion and the stand-off distance from the charge at which the pressure is measured. The peak pressure (*Pmax*), decay constant of the wave (θ), impulse (*I*), bubble oscillation period (*T*), the maximum radius of the first bubble of explosive gas (*Rmax*), and the energy flux density/energy per unit volume (*E*) can be expressed as per Equations (2)–(7) [39–41].

$$P\_{\max} = K\_1 \left(\frac{W^{\frac{1}{\lambda}}}{R}\right)^{A\_1} \quad (\text{Mpa}) \tag{2}$$

$$\theta = K\_2 W^{1/3} \left(\frac{W^{\frac{1}{3}}}{R}\right)^{A\_2} \quad \text{ (ms)}\tag{3}$$

$$I = K \mathfrak{z} \mathcal{W}^{1/3} \left(\frac{W^{\frac{1}{3}}}{R}\right)^{A\_3} \quad (\text{Mpa-sec})\tag{4}$$

$$T = K\_{\mathbb{S}} \frac{W^{\frac{1}{2}}}{\left(D + 9.8\right)^{3/6}} \quad \text{(sec)}\tag{5}$$

$$R\_{\text{max}} = K\_6 \left( \frac{W}{(D + 9.8)} \right)^{1/3} \quad (\text{m}) \tag{6}$$

$$E = K\_4 \mathcal{W}^{\frac{1}{3}} \left( \frac{\mathcal{W}^{\frac{1}{3}}}{(\mathcal{R})} \right)^{A\_4} \quad (\text{m} \cdot \text{kPa}) \tag{7}$$

where (*K*1, *K*2, *K*3, *K*4, *K*5, *K*6) and (*A*1, *A*2, *A*3, *A*4) are constants that depend on the explosive charge type with values as in [39,42]. These input constants are illustrated in Figure 2. Phenomena of the underwater explosions (UNDEX): shock wave, high pressure and bubble motion.

**Figure 1.** Various wall structures utilized for pressure hulls. (**a**) Reproduced with permission from [26], Copyright Elsevier, 2011. (**b**,**d**–**f**) reproduced with permission from [27], Copyright Elsevier, 2003. (**c**) reproduced with permission from [28], Copyright Elsevier, 2016.

**Figure 2.** Phenomena of the underwater explosions (UNDEX): shock wave, high pressure and bubble motion.

(*W)* denotes the mass of the explosive charge in (Kg), (*D*) denotes the charge depth in *m*, and (*R*) denotes the distance between the explosive charge and target in (m). In UNDEX events, there are two types of subsequent cavitations that occur: the bulk and the local cavitation. The latter is caused by the reflection of the shock wave at the free surface and must be taken into consideration during the analysis of the surface ship. When the shock wave impinges upon the structure, then the total pressure *Ptot*(*t*) at the fluid structure interface can be expressed according to Equation (8) [43]:

$$P\_{tot}(t) = P\_{in}(t) + P\_{\mathcal{C}}(t) + P\_{st} \tag{8}$$

where *Pin*(*t*) denotes the incident shock wave, *Pc*(*t*) denotes the scattered pressure, and *Pst* denotes the hydrostatic pressure. Accordingly, the local cavitation occurs in the water as the pressure drops to vapor pressure (about 0.3 psi) [31]. Thereafter, the cavitation collapses and reload the structure.

#### **3. Composite Failure Criteria**

An appropriate failure criterion is needed in order to find the maximum permissible load before lamina failure. Therefore, it is necessary to develop theories to compare the state of stresses and strains in materials [44–46].

#### *3.1. Maximum Stress Failure Criteria*

The maximum stress criteria are single mode failure criteria. Fracture occurs if stresses at the principal material coordinates are higher than their respective strength. The failure index (*IF*) is presented in Equation (9) [47]:

$$I\_F = \max\begin{cases} \sigma\_{11}/X\_t & \text{if } \sigma\_{11} > 0 \text{ or } -\sigma\_{11}/X\_t \text{ if } \sigma\_{11} < 0\\ \sigma\_{22}/Y\_t & \text{if } \sigma\_{22} > 0 \text{ or } -\sigma\_{22}/Y\_t \text{ if } \sigma\_{22} < 0\\ & |\tau\_{12}|/S \end{cases} \tag{9}$$

where *Xt*, *Xc*, *Yt*, *Yc* and *S* denote the ultimate longitudinal, transversal and shear strength constants, respectively. Also, (σ11, σ<sup>22</sup> and τ12) denote the applied longitudinal, transversal and shear stress components, respectively, and can be calculated using Equation (10) [48]:

$$
\left\{ \begin{array}{c} \sigma\_{11} \\ \sigma\_{22} \\ \tau\_{12} \end{array} \right\} = \left\{ \begin{array}{cc} Q\_{11} & Q\_{12} & 0 \\ Q\_{12} & Q\_{22} & 0 \\ 0 & 0 & Q\_{66} \end{array} \right\} \left\{ \begin{array}{c} \varepsilon\_{11} \\ \varepsilon\_{22} \\ \gamma\_{12} \end{array} \right\} \tag{10}
$$

$$\text{where } Q\_{11} = \left(1 - \nu\_{12}\nu\_{21}\right)^{-1} \text{E}\_{1\prime} \text{ Q}\_{12} = \left(1 - \nu\_{12}\nu\_{21}\right)^{-1} \text{E}\_{1}\nu\_{21\prime} \text{ Q}\_{22} = \left(1 - \nu\_{12}\nu\_{21}\right)^{-1} \text{E}\_{2\prime} \text{ Q}\_{66} = \text{G}\_{12}.$$

#### *3.2. Tsai-Hill Failure Criteria*

Tsai-Hill failure criterion assume that there is an interaction between longitudinal, transversal and shear strength in the damage progress. The Tsai-Hill failure criterion can be expressed as in Equation (11) [49]:

$$
\sigma\_{\rm 11}^2 / \mathcal{X}^2 + \sigma\_{\rm 22}^2 / \mathcal{Y}^2 - \sigma\_{\rm 11} \sigma\_{\rm 22} / \mathcal{X}^2 + \tau\_{\rm 2}^2 / \mathcal{S}^2 = 1 \tag{11}
$$

where σ11, σ<sup>22</sup> and τ<sup>12</sup> are the applied longitudinal, transversal and shear stress components, respectively, and can be calculated using Equation (10). Meanwhile, *X* and *Y* are the longitude and traverse strength, respectively, whether in tension or compression, which depends on the stress status in the laminates. *S* is the shear strength constant.

#### *3.3. Tsai-Wu Failure Criteria*

The Tsai-Wu failure criterion is the most generalized criterion that distinguishes between compressive and tensile strength. The criterion can be expressed as shown in Equation (12) [50,51]:

$$FI = \sigma\_{11} \left(\frac{1}{X\_t} - \frac{1}{X\_c}\right) + \sigma\_{22} \left(\frac{1}{Y\_t} - \frac{1}{Y\_c}\right) - \frac{\sigma\_{11}^2}{X\_t \times X\_c} - \frac{\sigma\_{22}^2}{Y\_t \times Y\_c} - \frac{\tau\_{12}^2}{S^2} \tag{12}$$

where *FI* is the failure index, *Xt*, *Y*t, *Xc*, *Y*c, σ11, σ22, τ<sup>21</sup> and *S* are as aforementioned. The failure occurs when the calculated stresses reaches the ultimate stresses and the *FI* reaches or exceeds the value 1 [52]. In finite element procedures, failure criteria are presented using a defined failure index and can be presented as per Equation (13):

$$FI = \frac{Stress}{Strength} \tag{13}$$

The results presented in this work are based on maximum stress, Tsai-Hill and Tsai-Wu failure criteria.

#### **4. Optimization and Geometrical Configuration of MICEPH**

In this work, the proposed form of submarine pressure hull is a multiple intersecting cross-elliptical hull, constructed from Carbon/Epoxy composite (USN-150) with stacking sequence [(α/−α)4] *s*. The optimization is performed using the ANSYS V14.5 (ANSYS, Canonsburg, PA, USA) parametric design language (APDL) to maximize buckling load and minimize the buoyancy factor under the constraints of failure strength and deflection (δ*max*) according to the design requirements. The model was built using SHELL99 for the shell and BEAM189 for the ring and long beams [53,54]. The material properties and the strength parameters are presented in Table 1 [55]. Figure 3 shows the proposed MICEPH utilized in this study. Figure 3 shows the multi-objective optimization procedure flow chart. The random design generation method (RDGM), which is a sub type of the sub problem approximation used in (ANSYS) will be considered in this study. On the other hand, Table 2 illustrates the results of the multi-objective optimization for MICEPH which indicates the optimal design point and objective function. The table contains the failure index (*FI*) for the maximum stress (FMAXF) and Tsai-Wu (FTWSR) failure criteria. The optimal objective function (*MOF*) is 0.0487, the buoyancy factor (*B.F*) is 0.197, and the buckling strength factor is 4.056, with a total hull weight of 394.57 kg. The optimal angle-ply orientation (α*)* is 49◦, with a layer thickness (*t)* equal to 1.294 mm. The major diameter (*Dmajor*) and minor diameter (*Dminor*) are 2.0 m and 1.6836 m, respectively. The intersecting angel (θ*)* is 40◦, and the radius (*R)* is equal to 0.50 m.

**Table 1.** Strengths of unidirectional composites and material properties of the sandwich components.

**Figure 3.** Multi-objective optimization procedure flow chart.

**Table 2.** The results of optimal design for pressure hull without core (Carbon/epoxy composite (USN)).


#### **5. Modeling and Simulation of MICEPH**

The fluid structure interaction (FSI) and UNDEX phenomena are implemented in non-linear finite element code ABAQUS V6.13 (Dassault Sys Simulia Corp, Providence, RI, USA)/Explicit. The MICEPH used in this study consists of five cross-ellipses, as demonstrated in Figure 4. The optimized structural hull and water domain are imported from ANSYS Program as (\* iges). Since the water domain does not have internal space for the pressure hull, it therefore needs to be modified to have identical internal shape for the hull. This can be simply done using the merge/cut function in the CAE model part in ABAQUS. The pressure hull was modelled using the shell elements. The stiffeners were added in the radial and longitudinal directions to provide structural integrity. The fourth node shell element (S4R) was used to model the pressure hull and stiffeners. An assemblage of 4-node acoustic tetrahedral elements (AC3D4) was used to represent the external fluid. The total horizontal length of the fluid model with the spherical ends was 11.56 m. The vertical length of the fluid domain was 7 m. The modeled water domain part includes the hull structure. The model and the charge mass were located at depth 100 m below the free surface. In dynamic problems that involve fluid and coupled solid medium, the interface between the two domains must be identical. Also, the water domain must have bulk modulus and density since it is acoustic domain. Figure 5 shows the finite element model and the meshing technique for the MICEPH surrounded by the fluid. The boundaries of the fluid around the MICEPH may cause shock wave reflection or refraction, which may cause a change in its superposition or cancellation by the incident wave [56,57]. To overcome this problem, the boundary condition of the fluid is executed as a non-reflective boundary condition during the analysis. The pressure hull was exposed to UNDEX produced by various amounts of explosives and offset distances (30 kg TNT at 5.5 m, 20 kg TNT at 5.5 m, 15.5 m and 20.5 m). In this study, the stand-off distance is located at the right side of the cross-elliptical submersible pressure hull. The stand-off point represents the location where the incident wave defined and represented by reference point1 (RP1). The location of the charge (source point) defined as reference point 2 (RP2), which represents the position of the charge as illustrated in Figure 6. There are three types of input parameters for UNDEX: physical charge, material and bubble model. Defining the UNDEX bubble in the interaction module was applied by defining the source point, stand off point which specifying the wave properties and charge depth. Taking into consideration that the infinite surface shouldn't reflect the shock wave; thus, acoustic impedance (non-reflecting) has to be applied. The initial boundary condition for the infinite surface of water domain (acoustic wave) was set to zero, which means that the domain is calm water, i.e., no interference [58].

**Figure 4.** Miceph.

(**b**) fluid domain

**Figure 5.** Finite element modeling of: (**a**) cross-elliptical submersible pressure hulls, and (**b**) fluid domain.

**Figure 6.** The location of charge (source point (RP2)) and stand-off point (RP1).

#### **6. Results and Analysis**

#### *6.1. Propagation of Shock Wave*

In this study, the MICEPH is subjected to shock wave and bubble pulse owing to UNDEX. Its performance depends mainly on the strength of the waves and the resilience of its structure. The stand-off point at node (*A*) was first struck by the shock wave. The charge is located at different stand-off distances—5.5 m, 15.5 m and 20.5 m—measured from the stand-off point. In the modelling procedures, two output variables were provided for the acoustic pressure: (i) the acoustic pressure (POR), which represents the total dynamic pressure in the wave formulation analysis including additional pressure induced by the incident and scattered waves; and (ii) the absolute acoustic pressure (PABS), which is the sum of the acoustic pressure and hydrostatic pressure. The acoustic pressure magnitudes are in Pascal [56]. Figure 7 demonstrates the POR at which the total dynamic pressure

instantaneously occurs at the time of explosion (zero time instant) due to total wave formulation. Also, the figure shows the spherical shape of the shock wave and the initial instantaneous shock wave propagation at the time of the explosion due to different explosive weights and stand-off distances. It can be observed that the maximum shock wave pressure determined was 3.471 <sup>×</sup> 107 Pa due to 30 kg TNT at offset distance of 5.5 m. Figure 7b–d illustrates the computed POR field distribution due to 20 kg TNT at offset distances of 5.5 m, 15.5 m and 20.5 m. The maximum shock wave pressure measured was 2.373 <sup>×</sup> 107 Pa, 6.777 <sup>×</sup> 106 Pa and 4.853 <sup>×</sup> 106 Pa at explosive offset distance 5.5 m, 15.5 m and 20.5 m, respectively. The location of the source and standoff point greatly affect the computed POR field distribution, propagation and its magnitude.

**Figure 7.** Computed POR field at the zero time instant due to total wave formulation.

Figures 8 and 9 present the computed POR field distribution at several time instants at the front and back sides of fluid domain. The propagation of the shockwave and the dynamic pressure in the surrounding water were clearly presented in these figures. The figures demonstrated that there are some radial POR waves that hits the structure and consequently leads to some deformations. Furthermore, it was observed that the dynamic pressure is affected by the reflections and the emissions from the pressure hull in addition to the incident field from the source point. When the shock wave hits the cross-pressure hull, it is reflected and generates negative pressure. Also, it is interaction with the incident wave decreases the dynamic pressure in the surrounding fluid. This is attributed to the fact that the water cannot withstand tension. That is why the total dynamic pressure acquired negative values as illustrated in the figures. Furthermore, cavitation occurs immediately after the incident shock wave hits the MICEPH. Additionally, local cavitation around the MICEPH was observed once the acoustic pressure declines to the steam pressure of the fluid. Subsequently, when the local cavitation disappears, the load of the fluid on the pressure hull generates vibrations. Additionally, Figure 9 illustrates that the shock wave propagates symmetrically from the stand-off point to the aft and fore,

while it expands in the perpendicular direction. The aforementioned results are in accordance with the results reported in [2,43,59]. Also, Figure 9 presents the volume changes in gas bubbles. If the oscillating gas bubble is close enough to the pressure hull, the differential pressure will be created. When the bubble decreases in volume (due to resistance to water flow close to the hull), that would result in bubble collapsing onto the hull and producing high speed water jet, which in some instances is capable of destroying or holing the pressure hull.

**Figure 8.** Computed POR field distributions at several time instants at the front side of the fluid domain (on the blast side).

**Figure 9.** Computed POR field at several time instants at the back side of the fluid domain (the side most remote from the charge).

#### *6.2. Responses of Submarine Pressure Hull to UNDEX*

The optimum MICEPH subjected to non-contact UNDEX will show three major response modes: (i) motion in the axial direction that makes the accordion motion, (ii) motion in a direction at right

angle to the fore-and-aft line of the MICEPH that makes the whipping mode parallel to the direction of the shock wave propagation, and (iii) motion in the vertical direction that makes the breathing motion perpendicular to the shock wave direction of the travel. Several locations were chosen in the model (location *A*, *B*, *C*, *D*, *E* and *F*) to demonstrate the responses in the MICEPH as illustrated in Figure 10.

**Figure 10.** Finite element model of Multiple. Intersecting Cross Elliptical Pressure Hull (MICEPH) and the locations of different test points.

#### 6.2.1. The Acceleration Response at Different Locations

First, the response of the MICEPH in the axial direction is illustrated in Figure 11. The figure presents the time history response of *Z*-axial acceleration (*A*3). UNDEX of 20 kg TNT at different offset distances of 5.5 m, 15.5 m and 20.5 m were discussed. At points E and F, located at the center of each end of the MICEPH, it was observed that the acceleration responses (*A*3) were in the opposite directions and occurs at the same instance. The peak value of the acceleration is 26 <sup>×</sup> 103 m/s2 at an offset distance of 5.5 m, and occurs at 4 m/s. This measurement then oscillates and decays after 4 m/s. In addition, while the offset distance increases, the acceleration decreases and acquires 8.3 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2 and 6.65 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2 at offset distances of 15.5 m and 20.5 m, respectively. This is attributed to increasing the standoff distance of the explosive charge. Also, the figure demonstrates that the bubble pulses show a minor impact on nodes *E* and *F*.

Similarly, Figure 12 plots the time history curves of *Z*-axial acceleration (*A*3) for nodes *C* and *D*. The curves illustrate that while the pressure hull stroked by the shock wave, node *C* moves to the left and node *D* moves to the right at the same time instance which causes the accordion motion. This is attributed to the propagated compressive pressure and its subsequent release on the MICEPH in the axial direction. The peak value at node *<sup>D</sup>* of the measured acceleration is 27.4 <sup>×</sup> 103 <sup>m</sup>/s2 in the negative *Z*-direction, owing to a 20 kg TNT charge and an offset distance of 5.5 m. Likewise, at node *C*, the peak value of acceleration is 23.9 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2, which occurs at 2 m/s. Also, as the offset distance increases at node *C*, the acceleration decreases and achieves 4.6 <sup>×</sup> 10<sup>3</sup> m/s2 and 2.8 <sup>×</sup> 103 m/s2 at offset distances of 15.5 m and 20.5 m, respectively. Similarly, the maximum acceleration at node *D* are 5.1 <sup>×</sup> 103 m/s<sup>2</sup> and 3.2 <sup>×</sup> 103 m/s<sup>2</sup> at offset distances of 15.5 m and 20.5 m, respectively. The aforementioned results demonstrate that the peak acceleration at nodes *C* and *D* were achieved before nodes *E* and *F*.

(**a**) Fore and aft Acceleration (*A*3) at offset distance 5.5 m. (**b**) Fore and aft Acceleration (*A*3) at offset distance 15.5 m.

(**c**) Fore and aft Acceleration (*A*3) at offset distance 20.5 m.

**Figure 11.** The acceleration-time history (*A*3) at nodes *E* and *F* due to 20 kg TNT on the pressure hull.

(**c**) Fore and aft acceleration (*A*3) at offset distance 20.5 m.

**Figure 12.** The acceleration-time history (*A*3) at nodes *C* and *D* due to 20 kg TNT on the pressure hull.

On the other hand, Figure 13 illustrates the acceleration response (*A*2) of the MICEPH in the vertical direction of nodes *E* and *F* under UNDEX of 20 kg TNT at offset distances of 5.5 m, 15.5 m and 20.5 m. It is observed that node E moves opposite to node *F* throughout the transient response of the MICEPH. Also, the figure illustrates the breathing motion caused by the expansion and subsequent contraction of the MICEPH. For instance, at node *E* and an offset distance of 5.5 m, the peak acceleration (*A*2) is 9.45 <sup>×</sup> 103 <sup>m</sup>/s2 in the *<sup>Y</sup>*-direction, while at node *<sup>F</sup>*, the peak acceleration is 8.96 <sup>×</sup> 103 <sup>m</sup>/s2 in the negative *Y*-direction and occurs at a time instance of 8 m/s. Furthermore, at higher offset distances of 15.5 m and 20.5 m, the peak measured accelerations are 7.16 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2 and 5.26 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s<sup>2</sup> in negative *<sup>Y</sup>*-direction at node *<sup>E</sup>*. While at node *<sup>F</sup>*, the peak accelerations are 6.43 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2 and 5.15 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2 in the *Y*-direction and occur at a time instance of 2 m/s.

(**c**) Vertical acceleration (*A*2) at offset distance 20.5 m.

**Figure 13.** The acceleration-time history (*A*2) at nodes *E* and *F* due to 20 kg TNT on the pressure hull.

Furthermore, Figure 14 shows the acceleration time histories (*A*2) in the vertical direction for nodes *C* and *D* which located at top and bottom of the MICEPH under the effect of 20 kg TNT charge located at offset distances of 5.5 m, 15.5 m and 20.5 m. It is revealed that the upper point (node *C*) moves in a direction opposite to that of the lower point (node *D*). The peak accelerations at node *C*, measured downward, are <sup>−</sup><sup>40</sup> <sup>×</sup> 103 <sup>m</sup>/s2, <sup>−</sup>10.4 <sup>×</sup> 103 <sup>m</sup>/s2, and <sup>−</sup>7.37 <sup>×</sup> 103 <sup>m</sup>/s2 at offset distances of 5.5 m, 15.5 m, and 20.5 m, respectively. Likewise, the peak accelerations at node *D*, measured upward, are 49.2 <sup>×</sup> <sup>10</sup><sup>3</sup> <sup>m</sup>/s2, 12.6 <sup>×</sup> 103 <sup>m</sup>/s2, and 8.25 <sup>×</sup> 103 <sup>m</sup>/s2 at offset distances of 5.5 m, 15.5 m and 20.5 m, respectively. The peak acceleration time histories (*A*2) at node *D* are higher than its counterpart at node C, while the frequencies of accordion and breathing motion are nearly the same. These results ensure that the accordion and breathing motion are directly correlated as per [2].

(**c**) Vertical acceleration (*A*2) at offset distance 20.5 m.

Similarly, Figure 15 shows the acceleration time histories (*A2*) at nodes *A* and *B*. It is observed that, once the shock wave struck the stand-off point at node *A*, the vertical movements of nodes *A* and *B* follow the same direction. Node *A* moves first, followed by node *B*. The peak acceleration at node *A* is <sup>−</sup>8 <sup>×</sup> 10<sup>3</sup> m/s2, <sup>−</sup>1.8 <sup>×</sup> 103 m/s2 and <sup>−</sup>1.3 <sup>×</sup> 103 m/s2 at offset distance of 5.5 m, 15.5 m and 20.5 m, respectively, and occur downward. Likewise, the peak acceleration at node *<sup>B</sup>* is <sup>−</sup>3.5 <sup>×</sup> 103 <sup>m</sup>/s2, <sup>−</sup>1.1 <sup>×</sup> 10<sup>3</sup> m/s2 and <sup>−</sup>0.81 <sup>×</sup> 10<sup>3</sup> m/s<sup>2</sup> at offset distance of 5.5 m, 15.5 m and 20.5 m, respectively, following the same direction at point *A*. Figure 16 plots the athwart acceleration time histories (*A*1) at nodes *A* and *B* which is the primary direction of shock wave propagation. It is revealed that the shock wave arrives at node *A* first then hits node *B*. The maximum athwart acceleration (*A*1) at node *A* occurs in the horizontal direction and its peak is about 52.4 <sup>×</sup> 103 m/s2 at an offset distance of 5.5 m. This is attributed to the release of the shock wave, occurrence of local cavitation and effect of bubble pulse. On the other hand, the second peak values, measured at offset distances of 5.5 m, 15.5 m and 20.5 m, occur owing to the uploading of the MICEPH. Similarly, the maximum athwart acceleration (*A*1) at node *B* also occurs in the horizontal direction, and its peak is about 29.4 <sup>×</sup> 10<sup>3</sup> m/s2, 10.9 <sup>×</sup> 103 m/s2 and 7.9 <sup>×</sup> 103 <sup>m</sup>/s2 at offset distance of 5.5 m, 15.5 m and 20.5 m, respectively. From the aforementioned analysis, it is concluded that the local cavitation has a major effect on athwart acceleration at standoff point at node *A*. Likewise, the uploading of the structure and the first bubble pulse also have a major effect on athwart acceleration. Additionally, for nodes *A* and *B*, the transverse responses are very severe at these locations, and the greatest acceleration occurs in the athwart direction, which is the main direction of shock wave with values of 52.3 <sup>×</sup> 103 m/s2 and 29.4 <sup>×</sup> 103 m/s2 at time intervals of 2 m/s and 4 m/s, respectively, followed by the vertical and longitudinal directions. While, for nodes *C* and *D*, which located at the top and bottom of MICEPH, the peak acceleration occurs in the vertical direction with values of <sup>−</sup>40 <sup>×</sup> 103 m/s2 and 49.2 <sup>×</sup> 103 m/s2 at time interval of 4 m/s, followed by the athwart and longitudinal directions. Furthermore, at nodes *E* and *F*, which located at the center of each

end of the MICEPH (fore and aft direction), the peak acceleration occurs in the longitudinal direction followed by the vertical and athwart directions.

(**c**) Vertical acceleration (*A*2) at offset distance 20.5 m.

**Figure 15.** The acceleration-time history (*A*2) at nodes *A* and *B* due to 20 kg TNT on the pressure hull.

(**c**) Athwart acceleration (*A*1) at offset distance 20.5 m.

**Figure 16.** The acceleration-time history (*A*1) at nodes *A* and *B* due to 20 kg TNT on the pressure hull.

#### 6.2.2. The Displacement Response at Different Locations

Figure 17 illustrates the deformation response in the translation direction (*U*2) at nodes *A* and *B* owing to 20 kg TNT at various offset distances (5.5 m, 15.5 m and 20.5 m). It is observed that the vertical displacement (*U*2) at nodes *A* and *B* move in opposite directions throughout the transient response and follow the same variation. In addition, the figure shows a rigid body translation of the MICEPH. Figure 17a shows the high-frequency response measured at nodes *A* and *B*, at an offset distance of 5.5 m. It is revealed that at a small offset distance, the structural wavelengths are longer than the acoustic wavelengths. Additionally, the surrounding fluid acts on the structure as a simple damping mechanism and the energy transported away from the structure through the acoustic radiation. Furthermore, Figure 17b illustrates the intermediate-frequency response at nodes *A* and *B*, at an offset distance of 15.5 m. At an intermediate offset distance, the structural wavelengths are nearly similar to the acoustic wavelengths. In addition, the surrounding fluid added mass and radiation damping on the structure of the hull. Moreover, Figure 17c demonstrates the low-frequency response at nodes A and B, at offset distance of 20.5 m. It is revealed that the structural wavelengths are nearly shorter than the acoustic wavelengths and the surrounding fluid added effective mass to the structure of the hull on the wetted interface. Figure 18 presents the translation (*U*2) in the vertical direction of two points located at the top (node *C*) and bottom (node *D*) of the amid-ship section at different offset distances (5.5 m, 15.5 m and 20.5 m) due to 20 kg TNT. It is observed that the upper point (node *C*) is moved in a direction opposite to the lower point (node *D*). The figure also illustrates that the displacement at node *C* is higher than its counterpart at node *D.* In addition, the pressure hull shows an elastic deformations and rigid body motions with a significant heaving and pitching response. These results match those achieved by [60]. Similarly, Figure 19 illustrates the axial displacement response (*U*2) in the vertical direction of points located at the center of each end of the MICEPH (*E* and *F*). It is observed that the two nodes move in opposite direction with almost similar frequencies forming the accordion motion. Additionally, the shape of the displacement response (*U*2) is significantly affected by the offset distance. While the offset distance increases, the frequency decreases. The aforementioned results agree well with that acquired by [40,43].

**Figure 17.** Displacement-time history responses (*U*2) at nodes *A* and *B* on the MICEPH.

(**c**) *U*2 at offset distance 20.5 m.

**Figure 18.** Displacement-time history responses (*U*2) at nodes *C* and *D* on the MICEPH.

**Figure 19.** Displacement-time history responses (*U*2) at nodes *E*, *F* on the MICEPH.

#### 6.2.3. Failure Analysis of MICEPH

To evaluate the maximum allowable load that can be sustained before lamina failure, an appropriate failure criterion is needed. Therefore, the failure of the MICEPH is analyzed in this study using different failure criteria (maximum stress and Tsai-Hill failure criteria). The failure was assessed through the calculation of the failure index. If the failure index achieved unity, then the material had failed. Figures 20 and 21 demonstrate the failure index distribution on the structure of the pressure hull due to 20 kg TNT at different offset distances of 5.5 m and 15.5 m at various time instants.

**Figure 20.** Failure criteria distributions on pressure hull due to 20 kg TNT and offset distance 5.5 m.

**Figure 21.** Failure criteria distributions on pressure hull due to 20 kg TNT and offset distance 15.5 m.

Generally, it is observed that the failure indices decrease with increasing stand-off distance. Furthermore, the maximum failure indices are measured around the ring and long beams. Moreover, it is revealed that failure is firstly initiated near the stand-off point. Figure 21 also shows the distribution of different failure criteria such as Azzi-Tsai-Hill, maximum stress, and Tsai-Hill criteria on the MICEPH. The results demonstrate that the difference among the different failure criteria is relatively too small. Figure 22 presents the maximum stress and Tsai-Hill failure-time histories for elements *A*, *B*, *C*, *D*, *E* and *F* on the structure of the pressure hull. The results illustrate that the maximum failure index occurs at offset distance of 5.5 m. The failure index owing to the maximum stress and Tsai-Hill criteria showed a similar trend with relatively small differences.

**Figure 22.** Maximum stress and Tsai-Hill failure-time history for elements *A*, *B*, *C*, *D*, *E* and *F* on the pressure hull.

#### **7. Conclusions**

In this study, the dynamic response of optimized multiple intersecting MICEPH under the effect of non-contact UNDEX was explored. The simulation technique was followed to avoid the expenses and complexity of physical tests. First, the multiple intersecting MICEPH subjected to hydrostatic pressure was optimized. Thereafter, using the optimum design results, the numerical simulation was carried out. Then, the response modes, breathing, accordion and whipping for the MICEPH subjected to non-contact UNDEX were discussed. Furthermore, the effects of applying various failure criterions such as maximum stress, Azzi-Tsai-Hill, and Tsai-Hill criteria for the damage initiation on failure strengths of pressure hull were studied. Based on the analysis and the simulation results, the following conclusions were drawn:


**Author Contributions:** E.F. and M.H.; methodology, E.F. and M.H.; software, E.F. and M.H.; validation, E.F.; M.E.; M.A.; formal analysis, E.F. and M.H. investigation, H.H. and D.W., M.E., M.A.; resources, E.F. and M.H.; data curation, E.F.; writing—original draft preparation, all, writing—review and editing, E.F.; M.E.; M.A.; E.F.; visualization, E.F. and D.W; supervision, E.F, H.H. and D.W., project administration, D.W.; funding acquisition.

**Funding:** The authors are grateful for the support of this research by the 13th Five Years Key Programs for Science and Technology Development of China (Grant No. 2016YFD0701300), the Chinese Natural Science Foundation (Grant No. 51405076), and Heilongjiang Province Applied Technology Research and Development Program Major Project of China (Grant No. GA16B301). Also, the authors are grateful to Military Technical College (Cairo, Egypt), Taif University (Taif, KSA) and Mansoura University (Mansoura, Egypt) for providing all the required facilities to carry out the present research.

**Acknowledgments:** The authors are grateful to Military Technical College (Cairo, Egypt), Taif University (Taif, KSA) and Mansoura University (Mansoura, Egypt) for providing all the required facilities to carry out the present research.

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

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Issues on the Vibration Analysis of In-Service Laminated Glass Structures: Analytical, Experimental and Numerical Investigations on Delaminated Beams**

#### **Chiara Bedon**

Department of Engineering and Architecture, University of Trieste, Piazzale Europa 1, 34127 Trieste, Italy; chiara.bedon@dia.units.it; Tel.: +39-040-558-3837

Received: 5 August 2019; Accepted: 7 September 2019; Published: 19 September 2019

**Abstract:** Load-bearing laminated glass (LG) elements take the form of simple members in buildings (i.e., columns, beams, and plates) or realize stand-alone assemblies, where glass and other traditional constructional materials can interact. Among several relevant aspects, the dynamic response of LG structures requires dedicated methods of analysis, towards the fulfilment of safe design purposes. A combination of multiple aspects must be taken into account for dynamic calculations of even simple LG elements when compared to static conditions, first of all the sensitivity of common interlayers to the imposed vibration frequency. The challenge is even more complex for the vibration serviceability assessment of in-service LG structures, where the degradation of materials and possible delamination effects could manifest, hence resulting in structural performances that can markedly differ from early-design conditions. Major uncertainties can be associated to the actual mechanical characterization of materials in use (especially the viscoelastic interlayers), as well as the contribution of restraints (as compared to ideal boundaries) and the possible degradation of the bonding layers (i.e., delaminations). All of these aspects are examined in the paper, with the support of extended analytical calculations, on-site experimental measurements, and parametric Finite Element (FE) numerical analyses. When compared to literature efforts accounting for ideal boundaries only, an analytical formulation is proposed to include the effects of flexible restraints in the dynamic performance of general (double) LG beams. Special care is also spent for the presence of possible delaminations, including size and position effects. In the latter case, existing formulations for composite laminates are preliminarily adapted to LG beams. Their reliability and accuracy is assessed with the support of test predictions and parametric FE simulations.

**Keywords:** laminated glass (LG); free vibrations; fundamental frequency; mechanical restraints; field experiments; analytical modelling; Finite Element (FE) numerical modelling

#### **1. Introduction**

Laminated glass (LG) is largely used in several ways, and mostly increasing is its application for load-bearing structural members for buildings and constructions. There, special care is required for the optimal and safe design of structural elements that are able to ensure resistance, robustness, stiffness, redundancy, etc., even under extreme loads [1,2].

From a structural point of view, critical design conditions can include explosive events [3–5] and impacts [6–8], natural hazards and earthquakes [9,10], and dynamic loads in general (including moving loads, in the case of pedestrian systems [11–16]), whose effects need dedicated calculation methods.

While the research community is spending efforts for the refinement or development of analytical/graphical design procedures that are able to capture the limit states of interest for design, the actual dynamic behaviour of even simple laminated glass elements is often described by means of approximate methods/working assumptions [17]. Such a strategy can be beneficial at the preliminary

design stage. However, the same approaches are often not able to properly capture the real dynamic performance of in-service LG structures. This is especially the case of LG systems that are subjected to severe operational conditions, where mechanical properties and boundary conditions might be strongly affected by installation methods, ambient, etc. (i.e., Figure 1).

**Figure 1.** Example of glass structures under severe operational conditions, due to (**a**) temperature variations or (**b**) overcrowding.

In this context, the paper aims at investigating the dynamic behaviour of the in-service LG members in free vibrations. Following earlier studies of literature, special care is focused on the sensitivity analysis of their fundamental frequency to a multitude of aspects of technical interest, including the type of interlayer (and its frequency-dependent shear stiffness), the aspect ratio of LG beams, the presence of flexible restraints that can provide only partial translational/rotational constraints (with respect to ideal boundaries), or the effects of possible delaminations in the bonding interlayers.

To this aim, classical analytical methods are first recalled in Section 2, so as to introduce the literature concept of "adjusted dynamic" effective thickness for LG beams [17]. Hence, preliminary analytical calculations are carried out on a wide set of LG beams, giving evidence of the expected sensitivity of frequency estimates to some variations in the LG beam composition and/or geometrical configuration.

As shown, major uncertainties in frequency estimates can derive from the actual stiffness of restraints, as well as from the dynamic response of the interlayers in use, to the imposed vibration frequency. Even further sensitivity is related—for in-service structures—to the possible presence of damage, like delaminations or material degradations, which should be properly taken into account for safe vibration assessments, with respect to early-stage design calculations. Such an aspect is further explored with the support of on-site vibration experiments (Section 3) that were carried out on selected LG beams belonging to an existing structure. The sensitivity of test predictions to major influencing parameters is discussed in Section 4, with the support of refined Finite Element (FE) numerical models (ABAQUS [18]) and parametric analytical calculations. Based on experimental observations and refined FE numerical models (ABAQUS [18]), practical analytical expressions are proposed in the paper, so that the expected vibration frequencies of a general (double) LG beam with flexible restraints could be rationally calculated. Dedicated analytical methods are also adapted to LG members and proposed to include possible delaminations in the bonding interlayers, thus resulting in more accurate dynamic estimations for in-service glass structures, based on past literature studies on composite laminates.

#### **2. Classical Analytical Formulation for Frequency Calculations**

#### *2.1. Reference System*

The attention of the current study is focused on the vibration performance of double LG beams agreeing with Figure 2. For simplicity, the reference cross-section is symmetrical, inclusive of two glass layers (with thickness *h*<sup>1</sup> = *h*3) and a middle viscoelastic foil (*h*2), providing a certain shear coupling.

The resisting LG member has total width *b*, with *L* the span, and *d* ≥ 0 representing the distance (if any) between each restraints (from the middle axis) and the beam end section. Accordingly, *L*<sup>0</sup> >> 2*d* is the actual bending span. The thickness *hi* of each layer is relatively small when compared to the global dimensions *b* × *L*.

In addition, *E*<sup>1</sup> = *E*<sup>3</sup> = 70 GPa is the nominal modulus of elasticity (MoE) of glass, with ρ<sup>1</sup> = ρ<sup>3</sup> = 2500 kg/m3 the density and ν<sup>1</sup> = ν<sup>3</sup> = 0.23 the Poisson' ratio [19]. Disregarding the interlayer type and composition, the bonding foil has generally a relatively low density ρ<sup>2</sup> as compared to glass (ρ<sup>2</sup> <sup>≈</sup> 1000 kg/m3), while its stiffness can strongly modify with operational conditions, see Section 2.2.

**Figure 2.** Double LG beam in free vibrations: (**a**) transversal and (**b**) longitudinal cross-sections, with (**c**) selected ideal restraints.

#### *2.2. Existing Closed-Form Solutions*

Assuming the ideal restraint configurations that were schematized in Figure 2c for simply supported (S-S) or clamped (C-C) members, classical theories for slender beams can be taken into account for frequency analyses of LG elements, as far as the dynamic mechanical properties of the interlayer in use are properly described.

According to Figure 2, a given LG beam in free vibrations must, in fact, satisfy the well-known Euler–Bernoulli differential equation of motion, that for a given a monolithic (*A* = *b* × *h*) × *L*<sup>0</sup> member is given by [20]:

$$\frac{\partial^2}{\partial \mathbf{x}^2} EI(\mathbf{x}) \left( \frac{\partial^2 \upsilon(\mathbf{x}, t)}{\partial \mathbf{x}^2} \right) + \rho A \frac{\partial^2 \upsilon(\mathbf{x}, t)}{\partial t^2} = 0 \tag{1}$$

Moreover, in Equation (1), *v*(*x*,*t*) is the vertical displacement at the abscissa 0 ≤ *x* ≤ *L*<sup>0</sup> and time instant *t*; *E* and ρ the MoE and density of the material in use; and, *I* the second moment of area. Thus, the fundamental frequency is conventionally given by the compact expression [20]:

$$f\_n = \frac{\omega\_n}{2\pi} = \frac{1}{2\pi} \sqrt{\frac{\beta\_n^4 E}{12\overline{m}}} h^3 \tag{2}$$

with *n* the mode order and β the wavenumber (Table 1).


**Table 1.** Reference wavenumbers β*<sup>n</sup>* and shape coefficients Ψ for beams with simple restraints (Figure 2c) and bending span *L*0.

According to several studies of literature, the challenge for a given LG member lies in the estimation of the actual composite stiffness, being strongly related to the shear coupling effect of the bonding interlayer. Besides the availability of simplified analytical approaches that are based on the use of an equivalent, monolithic glass thickness *hef* = *h* for sandwich sections, according to Figure 2, it was shown in [17] that the "adjusted dynamic" effective thickness (adapted from [21] for modal analysis purposes) is able to offer reliable frequency estimates for double LG beams with ideal boundaries. In particular, such an "adjusted" thickness is given by [17]:

$$h\_{ef} = \sqrt[3]{\frac{1}{\frac{\eta}{h\_1^3 + h\_3^3 + 12I\_\varkappa} + \frac{1-\eta}{h\_1^3 + h\_3^3}}} \tag{3}$$

where η represents the shear coupling of the composite section:

$$10 \le \eta = \frac{1}{1 + \frac{E\_1 h\_2}{C\_2 b} \cdot \frac{I\_1 + I\_3}{I\_{\text{lat}}} \cdot \frac{A\_1 A\_3}{A\_1 + A\_3 \cdot \Psi}} \le 1 \tag{4}$$

and the other relevant terms are given by:

$$I\_i = \frac{bh\_i^3}{12} A\_i = bh\_i \tag{5}$$

$$I\_s = \frac{h\_1 h\_3}{h\_1 + h\_3} \cdot [h\_2 + 0.5(h\_1 + h\_3)]^2 \tag{6}$$

$$I\_{\rm tot} = I\_1 + I\_3 + \frac{A\_1 A\_3}{A\_1 + A\_3} \cdot \left[ h\_2 + 0.5(h\_1 + h\_3) \right]^2 \tag{7}$$

The coefficient Ψ in Equation (4) depends on the normalized shape of deflections, for a given homogeneous beam. For basic boundary conditions and several mode orders *n*, Ψ can be calculated from Table 1.

Finally, η in Equation (3) is strictly affected by the shear modulus *G*<sup>2</sup> = *G*2(ω) of the bonding layer. Given that the common materials for LG applications have a viscoelastic behaviour that depends on the material composition and its vibration frequency and/or ambient conditions, this turns out in an effective thickness *hef* = *hef*(ω), which explicitly reflects the dynamic response of the interlayer itself, as a part of a composite system it belongs. However, according to Equation (2), it is also ω = ω(*hef*), and hence an iterative calculation approach is required for accurate thickness/frequency estimates.

Both the real, frequency-independent term (storage modulus *G*2,0) and the imaginary part (loss modulus, *G*2,ω), are in fact involved in the frequency domain, where:

$$\mathcal{G}\_{\mathsf{Z}}(\omega) = \mathcal{G}\_{\mathsf{Z},0} + \mathcal{G}\_{\mathsf{Z},\omega}(\omega) \tag{8}$$

and their typically high sensitivity to frequency is shown in Figure 3 (selected examples reproduced from [7,17,22]).

In this regard, literature projects have been dedicated to the mechanical characterization of interlayers in use for LG systems, under the assumption of various severe conditions of temperature or time loading (i.e., [23–30]). As far as the interlayer composition and the test method both modify; however, different mechanical properties can be derived for a given interlayer material [31]. The study reported in [32] also emphasized that the mechanical properties of the interlayer samples (i.e., material test), or interlayer foils belonging to small portions of LG sections (i.e., section test), can result in markedly different stiffness estimates, due to variation of the actual boundary conditions. Finally, for in-service glass structures, it is generally recognized that the degradation of interlayers can affect several material features, including the shear stiffness, but also the adhesion properties, and other thermo-mechanical parameters that could indirectly affect the overall structural performance of a given LG section (see [33–35]).

**Figure 3.** Examples of (**a**) storage and (**b**) loss moduli variation with frequency, for Polyvinyl Butyral (PVB) or SentryGlas Plus (SGP) interlayers at 25 ◦C, according to [7,17,22].

In other words, Equation (3) represents a practical tool for design, being able to simplify the original dynamic problem of composite beams with flexible, viscoelastic connection and ideal boundaries. On the other side, the reference *G*<sup>2</sup> = *G*2(ω) value for dynamic estimates should be properly assessed, including possible delaminations (Section 4).

In Figure 4, frequency calculations are proposed for selected LG beams in the S-S or C-C conditions, as a function of *G*2. The collected frequency values are derived from Equation (2), while assuming that *hef* (Equation (3)) modifies with *G*2, and 10<sup>−</sup><sup>4</sup> MPa < *G*<sup>2</sup> < 10<sup>5</sup> MPa. In the figures, both the limit "*layered*" and "*monolithic*" conditions can be easily detected. In addition, the grey regions represent all the possible frequency values that could characterize the dynamic performance of a given LG beam geometry, as far as its end restraints are characterized by a certain translational/rotational stiffness that can be comprised within the limit conditions of ideal simple supports (S-S) or clamps (C-C).

**Figure 4.** Example of analytical frequency estimations for double laminated glass (LG) beams, as a function of the restraint type (S-S or C-C) and shear stiffness of the interlayer. The grey region denotes the sensitivity of analytical frequency estimates to the restraints and interlayer stiffnesses.

As far as higher vibration modes are taken into account in Figure 4, it is possible to notice that the grey region progressively minimizes for *n* = 2 and *n* = 3. As such, major uncertainties for simplified analytical calculations can be expected, especially for *n* = 1, that in most of the cases is a key parameter for design purposes.

#### *2.3. Restraints and Delaminations for In-Service LG Systems*

A relevant influencing parameter for the vibrational analysis of LG beams is certainly represented by the effect of real restraints, with respect to the ideal supports (Figure 2c).

It was shown in [12], for example, that the restraints characterized by a certain flexibility (i.e., axial (*Ks*) and rotational (*Kr*) stiffness due to the presence of soft layers, gaskets, etc., see Figure 5) should be properly taken into account for the dynamic analysis of even monolithic glass members with cantilever or beam behaviour, affecting both frequency and damping calculations.

**Figure 5.** Example of LG beams with end mechanical restraints: (**a**) real system and (**b**) corresponding mathematical model.

The presence of delaminations and their effects on the flexural stiffness of the resisting sandwich section can represent another crucial aspect to account for dynamic analysis of LG members. Critical regions for delaminations are commonly represented by restraints and edges (see, for example, Figure 6 and [36,37]). However, research studies on the structural performance of delaminated LG sections are still limited (i.e., [38–41]). In addition, literature investigations are focused on the stress response of simple LG members under laboratory conditions, which is with artificially imposed delaminations.

**Figure 6.** Examples of severe delamination in LG member, in the region of point mechanical restraints. Reproduced with permission from (**a**) [36] and (**b**) [37].

In this regard, Equation (2) is expected to provide only approximate estimates of the actual dynamic performance of in-service LG members, thus suggesting the need of more refined methods of analysis.

#### **3. Experimental Study on In-Service LG Beams**

A series of field experiments was carried out on existing LG member in order to further explore the actual dynamic behaviour of glass structures in operational conditions.

#### *3.1. Specimens and Test Methods*

The experimental study was carried out in May 2019, on a selection of 32 LG beams composed of two, *h*<sup>1</sup> = *h*<sup>3</sup> = 10 mm thick fully tempered glass layers and a middle PVB foil (*h*<sup>2</sup> = 0.76 mm its thickness). The width of LG beams was fixed (*b* = 0.14 m), while the variations were represented by the total span *L*. For the majority of them, *L* was in the range of 2.4 m and 2.7 m. The minimum span—even with identical nominal section properties and restraints—was in the order of 1.45 m. Such a marked variability in the span was required—at the design stage—to accommodate some geometrical irregularities of the primary structure hosting the specimens. The examined LG beams are, in fact, currently part of an in-service glass walkway (in the form of handrails, see Figure 7), being constructed in the early 2000 in the context of a Roman age Basilica monument in Aquileia, Italy (see also [15,16] for further details on the pedestrian system).

**Figure 7.** Example of service loads for the tested LG beams (photos by C. Bedon, courtesy of So.Co.Ba.).

For all of the LG members, the end restraints were realized in the form of stainless steel point-fixings according to Figure 8a, with holes having 42 mm nominal diameter and positioned at a distance *d* = *b*/2 = 70 mm from the edges (see Figure 8b).

**Figure 8.** Experimental LG specimens: (**a**) detail of the typical restraint (photo by C. Bedon, courtesy of So.Co.Ba.) and (**b**) schematic representation of the test setup test setup.

At the time of the on-site experiments, the selected LG beams were subjected to induced vibrations and their acceleration in time was monitored via a single tri-axial sensor [42], glued on the top surface of the mid-span section (Figure 8b). All of the tests were carried out with a mean temperature of 23 ◦C and a 62% relative humidity.

For each specimen, multiple measurements were collected (minimum three test repetitions), and then post-processed to predict the corresponding (mean) frequency.

In this regard, the advantage and potential of Operational Modal Analysis techniques is represented by the possibility to derive even relevant mechanical parameters for in-service structures that cannot be subjected to destructive (or laboratory) experiments. Otherwise, for the LG structure object of analysis, the serviceability of the religious monument did not allow for performing massive experimental

measurements, and required the use of a minimum number of instruments, thus involving some further uncertainties for the interpretation of test results.

A key influencing parameter, for example, was represented by a certain degradation of the PVB interlayers (due to repeated non-controlled ambient conditions and time), with visible delaminated regions, close to the restraints and along the edges of each beam (detail views are proposed in Figure 9 for some of the tested specimens).

Another major issue in the experimental study and vibration serviceability assessment consisted in the actual life-time of the selected LG specimens, thus in additional difficulties for the reliable estimate of PVB mechanical properties. Most of them were characterized by a mean service life of ≈15 years at the time of the research study. However, some of them have been replaced during the years, without any track of maintenance/replacement interventions.

Finally, a further uncertainty was represented by the actual stiffness contribution of the steel point-fixings in use (Figure 8), thus its effects on the overall dynamic response of the LG members.

**Figure 9.** Example of delaminations (i.e., bubbles and shadows) for a selection of tested LG beams (photos by C. Bedon, courtesy of So.Co.Ba.).

#### *3.2. Derivation of Experimental Fundamental Frequencies*

The analysis of the experimental results was based on post-processing of the collected acceleration-time data, see Figure 10. Given the availability of a single control point only for each test specimen, special care was spent for the fundamental frequency of LG beams, disregarding higher experimental modes, or vibration shapes and damping related issues.

**Figure 10.** Examples of test records: (**a**) time-acceleration data and (**b**) Power Spectral Density (PSD) function for selected specimens.

All of the LG specimens proved to offer a beam-like behaviour, but strictly related to the effect of mechanical supports and PVB layers. In Figure 11, for example, the experimental frequencies are proposed for the tested beams, as a function of *L*0. Analytical estimates from Equation (2) are also collected, as obtained for LG members having nominal experimental dimensions, but classical boundaries (S-S or C-C). Disregarding the uncertain PVB shear stiffness of the test specimens, moreover, three configurations are analytically taken into account in Equation (2) and Figure 11, that is the limit "*layered*" and "*monolithic*" theoretical configurations, and the PVBA properties of Figure 3.

It is worth of interest that most of the test predictions are comprised within the lower limit of S-S beams with "rigid" shear connection ("*monolithic*" curve of Figure 11a) and the upper limit of C-C beams with "weak" mechanical bonding between the glass layers ("*layered*" curve of Figure 11b).

**Figure 11.** Experimental and analytical frequency estimates for the examined LG specimens, grouped by *L*0, with (S-S) or (**b**) clamped (C-C) configurations.

#### **4. Analysis of Relevant Influencing Parameters**

#### *4.1. Sti*ff*ness Contribution of Point-Fixings*

The actual stiffness contribution of the joints in use was first assessed, with the support of FE numerical models (ABAQUS). In Figure 12, the reference numerical model is shown, being representative of the nominal geometry for a connection detail. A small portion of glass (0.14 × 0.14 m) was also taken into account to reproduce the actual joint region and interactions. In doing so, a set of three-dimensional (3D) solid brick elements was used to describe the steel restraint (C3D8R type from ABAQUS element library). The mesh size and pattern was chosen to ensure a refined description of the examined system, with 3300 solid elements for the steel joint and the portion of LG plate.

The mechanical interaction between the steel connector and the LG portion was accounted in the form of a *penalty* & *normal behaviour* surface-to-surface contact algorithm (ABAQUS library), so as to allow for possible relative sliding between the steel and glass components in the region of glass hole, but also the possible separation under tensile loads. Rigid nodal restraints were applied at the bottom face of the steel connector. At the same time, the LG plate was restrained in its thickness and width (end section), so as to avoid possible rotations, being a part of a full beam in bending about the minor axis.

A preliminary static nonlinear analysis was carried out to estimate the expected stiffness contributions under the assumption that the so assembled FE model can be representative of the actual end region for one of the tested specimens. Based on Figure 12, the bending performance of the small scale FE model was, in fact, explored by imposing a linearly increasing bending moment *Mz* at the top face of the steel point-fixing (via a reference RP node, and hence distributed on the full surface of steel with a "*coupling*" constraint).

Hence, the longitudinal (i.e., *x* direction) and vertical (i.e., *y* direction) displacements of four selected control points were monitored, so that the corresponding elastic stiffnesses (both rotational and translational/axial) could be properly calculated. Two of these control points, see Figure 12, were set at the edges of the hole, to monitor the bending rotation of the portion of LG plate, with respect to the steel connector.

**Figure 12.** Local numerical model of the reference joint (ABAQUS). (**a**) Assembly details (with hidden mesh pattern) and (**b**) bending deformation.

Globally, the translational stiffness *Ks* was found to be relatively high, hence justifying the assumption of an ideal rigid joint for comparative purposes (i.e., *Ks* = ∞). Given the reference system of Figure 12, otherwise the rotational stiffness was estimated in *Kr* ≈ 150 kN m/rad.

The so-calculated stiffness was first used form some preliminary FE calculations, carried out in ABAQUS with the support of Equation (3), to explore the sensitivity to *Kr* of the experimentally predicted frequencies. In Figure 13a, comparative estimates are shown for selected LG specimens. Frequency estimates are also proposed for the limit configurations of C-C and S-S beams. In general, it is possible to notice that "ideal" analytical calculations tend to result in marked percentage scatter, with respect to on-site experiments (Figure 13b).

For the S-S calculations, the FE predictions resulted in a mean Δ = −29%, as compared to the test data (with Δmin = −14%, Δmax = −40%, Dev.St = ±7). Given that the S-S boundary fully disregards the rotational role of the point-fixings in use, such an assumption is generally expected to severely underestimate the actual bending stiffness of a given composite system, thus the corresponding frequency. On the other side, see Figure 13b, the C-C assumption gave evidence of a marked overestimation of test results, with a mean Δ = +59% (Δmin = +35%, Δmax = +92%, Dev.St = ±11).

For the calculations that were carried out with partially rigid restraints, finally the mean percentage scatter was found in the order of Δ = +50% (with Δmin = +27%, Δmax = +82%, Dev.St = ±10), thus suggesting an improved agreement with, but still recommending more detailed analyses to assess, the actual dynamic performance of the tested LG specimens. The poor correlation of Figure 13 could be, in fact, justified by the actual shear contribution of the PVB layers, which, besides the availability of shear/loss moduli of literature (i.e., Figure 3), can be strongly affected by severe operational conditions, including high temperature and humidity variations, and medium/long-term degradation phenomena. The presence of visible delaminations that were visually detected for most of the tested LG beams (i.e., Figure 9) is another key parameter that affects the composite stiffness of the specimens, but it can be hardly quantified with detail.

**Figure 13.** Numerical (with PVBA interlayer properties) and experimental frequencies for the tested LG beams: (**a**) comparative frequency values and (**b**) corresponding percentage scatter Δ, as a function of *L*0.

#### *4.2. Derivation of Practical Fitting Curves for LG Members with Flexible Restraints*

According to literature, the definition of closed-form formulations accounting for the actual mechanical restraint in viscoelastic LG beams in free vibrations can involve complex mathematical problems (i.e., [43]). Otherwise, simple analytical fitting curves can offer robust support for reliable frequency estimations.

Assuming that a given glass member is restrained via mechanical supports having partial rotational stiffness *Kr* 0 (in kNm/rad) and *Ks* = ∞, the restraints act as flexible clamps that must be properly taken into account for design. Moreover, following Section 4.1, it is convenient to express these rotational and axial stiffness contributions as:

$$R\_r = \frac{K\_r L\_0}{EI} \tag{9}$$

and

$$R\_s = \frac{K\_s L\_0^3}{EI} \tag{10}$$

where *Rr* = 0 corresponds to the limit condition of a S-S beam and *Rr* = ∞ denotes the C-C configuration. As far as a monolithic or LG beam is taken into account in Equations (9) and (10), its bending stiffness *EI* can be expressed as a function of the total *b* × *h* or equivalent *b* × *hef* section (Equation (3)).

Thus, from a practical point of view, the vibration frequency of a given glass beam with non-ideal restraints can be conveniently calculated as (with *Rr* ≥ 0 and *Rs* = ∞):

$$
\overline{f}\_n = k\_f f\_n \tag{11}
$$

where *f <sup>n</sup>* is given by Equation (2) for a S-S beam and *kf* represents a magnification factor, depending on *Rr*.

In this paper, a series of parametric FE numerical simulations was carried out in ABAQUS for glass beams with different geometrical features (*b* = 0.1–0.5 m, *h* = 0.005–0.04 m, *L* = *L*<sup>0</sup> = 1–4 m) and variable rotational stiffness *Kr* (with *Ks* = ∞). In doing so, beam (B31) type elements were used, with equivalent springs being able to reproduce the desired rotational stiffness *Kr* (Figure 14).

**Figure 14.** Parametric numerical analysis of monolithic glass beams (ABAQUS): (**a**) axonometric view of the expected vibration shape (*n* = 1) and (**b**) qualitative shape variation, as a function of *Rr*.

Analytical fitting curves are thus proposed in this paper, in support of practical calculations (i.e., based on Equation (11)) on general LG beams with flexible restraints, see Figure 15 and Table 2.

As shown in Figure 15, in particular, it was observed that the use of flexible restraints and their combination with different geometrical/mechanical LG properties can result in even marked variations of the expected vibration frequencies. This is especially the case of the fundamental mode (*n* = 1), where relatively stiff restraints can amplify up to ≈2.25 times the S-S estimations. Certainly, the vibration shape is also expected to modify with the stiffness variation of restraints (i.e., Figure 14b). However, the current investigation was specifically focused on frequency estimates, and the accurate analysis of shape sensitivity would require more detailed experimental methods, as compared to the available test predictions.

**Figure 15.** Analytical fitting curves (see also Table 2) for the frequency estimation of glass beams with flexible restraints (*Kr* > 0 and *Ks* = ∞).


**Table 2.** Key input parameters for the magnification factor *kf* (Equation (11)), given that *Kr* > 0, *Ks* = ∞, and *Rr* in Equation (9)).

Given that *hef* in Equation (3) modifies with *G*2(ω), another relevant outcome of Figure 15 is that the role of real restraints can be properly taken into account in the iterative analytical procedure recalled in Section 2. Marked variations of *hef* (and thus frequency estimates), as in Figure 16, can in fact be obtained within the "*layered*" and "*monolithic*" limit conditions, as far as *Rt* modifies.

Figure 16a, in particular, shows the evolution of *hef* (Equation (3) for selected LG members bonded with PVBM foils. The fundamental mode (*n* = 1) is the most sensitive to restraint parameters. As far as stiff interlayers are used for a given LG beam, a mostly "*monolithic*" equivalent section *hef* can be also expected, with minimum sensitivity to restraints, as in Figure 16b. On the other side, the thickness and frequency sensitivity to the restraints in use progressively increases in presence of weak interlayers (as it is in presence of severe operational conditions). Moreover, from Figure 16, it is also possible to perceive that frequency estimates based on the simplified assumption of "*layered*" or "*monolithic*" thicknesses with ideal supports (*R*<sup>R</sup> = 0 or *R*<sup>R</sup> = ∞) can result in mostly rough (and even unconservative) calculations, with respect to the real performance of a given LG system.

**Figure 16.** Variation of the dynamic thickness *hef* (Equation (3)) of LG beams, as a function of (**a**) the order of the vibration mode *n* or (**b**) by changing the interlayer properties.

#### *4.3. E*ff*ect of Delaminations*

The vibration analysis of composite beams with delaminations is another complex issue, which can be generally solved with the use of advanced FE models (i.e., [44–47]), or with the support of coupled numerical-experimental studies on artificially damaged specimens (see, for example [48–50]). In both the cases, besides different composite laminates are taken into account, a common outcome of existing literature projects is represented by the high sensitivity of vibration frequencies to delaminations. Otherwise, dedicated studies for delaminated LG structures are still missing.

#### 4.3.1. Analytical Description of the Problem

Given a composite beam with delaminations, the mathematical problem can be solved in accordance with [51], i.e., in the form of an effective longitudinal modulus (*Eef*) representative of the actual mechanical interaction between the constituent layers. For a laminated section. according to Figure 17a, such a modulus is given by:

$$E\_{cf} = \frac{8}{h^3} \sum\_{j=1}^{m/2} (E\_x)\_j \left( z\_j^3 - z\_{j-1}^3 \right) \tag{12}$$

with *Ex* the MoE of the *j*-th layer, *m* is the total number of layers composing the beam section (with *b* × *h* total dimensions); *zj* is the distance between the outer face of the *j*-th layer and the neutral plane of the section. When *m* is an odd number—as in the case of the examined double LG sections, see Figures 2 and 17b—Equation (12) is still valid, given that the external *j*-th layer is subdivided into two symmetric parts.

**Figure 17.** LG beams with delaminations: (**a**) reference analytical model for simply supported, delaminated composite beams (reproduced from [51]) and (**b**) detail example for a double LG section.

The effect of delaminations is, in fact, accounted for by reducing *Eef*, as a function of the number and size of delaminated regions. Based on [51], the resulting longitudinal MoE of a partially delaminated composite section is, in fact, given by:

$$E\_d = \left(E\_{zd} - E\_{ef}\right)\frac{A\_d}{A\_t} + E\_{ef} \tag{13}$$

with:

$$E\_{zd} = \frac{\sum\_{j=1}^{s} E\_{cf} z\_j}{z} \tag{14}$$

Moreover:

*Ezd*—is the longitudinal MoE of a totally delaminated section (along one or more interfaces), representative of the so-called imperfect effective MoE;

*S*—is the number of sub-layers detected by the delamination;

*zj*—represents the thickness of the *j*-th sub-layers; and,

*Ad*, *At*—are respectively the delaminated and total interfacial area between the bonded layers.

Key assumptions of Equations (12)–(14) are that:


Based on point (d), it is expected that the analytical method herein recalled could overestimate the bending stiffness decrease of the delaminated LG beams, thus resulting in conservative frequency predictions. Another issue can be related to the location of delaminated regions (disregarded by the surface parameter *Ad*).

For general LG members, which were included the tested specimens, delamination phenomena typically occur along the edges and close to restraints (Figure 9), where the LG layers are not protected and/or properly sealed. In some other cases, the spotted delaminated regions can be recognized over the surface of a given LG member. However, in both the cases, natural delaminated regions (and their quantitative effects for structural calculations) can be hard to visually detect and properly quantify.

4.3.2. Reliability of Frequency Calculations for Delaminated LG Specimens

The application of Equations (12)–(14) to the examined LG beams, where even marked frequency reductions can be expected, as compared to the theoretical values of LG beams with uniform PVB bonding.

Generally, the MoE of partially delaminated beams (*Ed*) is linearly dependent on the delaminated area, when compared to the total bonding surface (see Equation (13)). As far as the effective MoE is taken into account for the examined LG beams (i.e., Figure 17b), it is thus possible to expect a MoE decrease in the order of 30%, for even a limited portion of delaminated interlayer. Such a MoE variation turns out in a marked variation of the bending stiffness for the composite LG beam, and at the same time reflects on a different shear stiffness from the interlayer in use (i.e., Figure 3), thus representing a relevant influencing parameter for calculations.

A parametric study was hence carried out on selected configurations in order to assess the reliability of Equations (12)–(14) for simple analytical calculations on delaminated LG beams. In accordance with Figure 18, more in detail, schematic delaminations were defined, to reproduce—even in a simplified way—some of the on-site qualitative observations from the experimental tests.

In doing so, symmetry was taken into account for all of the possible configurations, thus resulting in:


**Figure 18.** Selected delamination schemes for the parametric study on LG beams. Dashed area are representative of delaminations.

The analytical estimates of fundamental frequencies were carried out with the dynamic thickness *hef* and iterative approach of Section 2, with the additional set of iterations due to a decreasing MoE (with increasing the delamination surface *Ad*). The same iterative approach includes the effect of mechanical restraints, in accordance with Equation (11).

The support of full 3D solid models from ABAQUS was also taken into account. In the latter case, the choice of full 3D models was suggested by the need of including possible delaminations in several regions of interest (see Figure 19). The glass and interlayer foils were rigidly connected via "*tie*" constraints, in order to ensure a rigid mechanical connection for the involved surfaces. In the presence of delaminated regions, moreover, additional contact interactions were used between the involved layers. In Figure 19a, an example is proposed for the D4 scheme. The mechanical restraints were then reproduced in the form of equivalent axial and rotational springs (with *Kr* > 0 and *Ks* = ∞). These springs were connected to reference RP nodes, and then restrained to each LG beam with the use of additional "*coupling*" constraints, thus reproducing the actual effect of the steel point-fixings in use (see Figure 19b)). Frequency analyses were thus carried out on a wide set of configurations of LG beams. In Figure 19c, an example of typical deformed shape is shown (*n* = 1).

Through the FE parametric calculations, variations were made in terms of span (*L*) and delaminated region (*Ad*).

Figure 20 reports some of the collected results, for a LG beam with *L* = 2.65 m and *d* = 0.07 m, 10 mm and 0.76 mm the thicknesses of glass and PVB layers, respectively. The interlayer properties were defined as for the PVBA material law reported in Figure 3.

From Figure 20a, in particular, it is possible to notice that delaminations close to the restraints (i.e., schemes D1 to D3 of Figure 18) can be mostly disregarded. Frequency variations with respect to the vibration frequency of the undamaged (fully bonded) LG members are, in fact, proposed as a function of *Ad*/*At*, as obtained from Equations (12)–(14), giving evidence of an expected frequency decrease in the order of −0.2%. This finding is also in line with literature efforts (for different typologies of composite laminates, see, for example [48–50]), where it was proven that (artificially imposed) delaminated regions close to restraints have minimum effects.

**Figure 19.** Numerical modelling of LG beams with delaminations: (**a**) delaminated regions, (**b**) end detail, and (**c**) typical deformed shape (ABAQUS).

Otherwise, as far as the damage location moves along the LG beam span (i.e., scheme D4), a mostly linear frequency decrease was observed from the collected parametric results, see Figure 20b. In this sense, the parametric investigation summarized herein proved that:


**Figure 20.** Frequency variation of delaminated LG beams with flexible mechanical restraints (with PVBA interlayer, *Kr* = 150 kNm/rad and *Ks* = ∞). (**a**) D1-to-D3 or (**b**) D4 scheme results, for selected LG beams.


#### 4.3.3. Final Remarks on Practical Analytical Calculations for Design

In conclusion, some final analytical calculations were carried out, in order to assess—even in the lack of more detailed experimental methods—the actual bonding effect of the PVB foils in use for the tested LG specimens, with respect to practical literature recommendations. Given that the fundamental frequency of the selected LG beams is strictly related to several parameters, the degradation of shear stiffness for the PVB foils in use—as a direct effect of long-term performances—represents another relevant parameter, as also discussed in the previous sections.

Literature contributions addressing the actual dynamic shear modulus of PVB foils (including both the storage and loss moduli), however, are available for dynamic estimates at different vibration frequencies (i.e., Figure 3), or different temperature scenarios, which can be hardly adapted to specific case-studies. Long-term effects of design loads (up to 50 years) are mostly referred to the relaxation of the storage modulus of PVB (and other common interlayers) under permanent loads only (see for example [23] and others).

In this regard, Figure 21 collects some analytical estimates for the tested LG beams. Following the research outcomes that are partly summarized in Figure 20, and the real experimental scenarios of Figure 9, three levels of delamination were taken into account, i.e., corresponding to (i) un-delaminated LG beams, or to a MoE degradation up to (ii) 5% and (iii) 15% the nominal value (based also on Equation (13)). In addition, two tentative storage moduli were considered for the PVB foils in use, namely represented by:


As shown in Figure 21 and Table 3, the best correlation between analytical frequencies and most of the experimental results was obtained in the presence of a relatively weak PVB layer (*G*<sup>2</sup> = 0.07 MPa).


**Table 3.** Percentage scatter Δ (mean value) of analytical frequency estimates, compared to experimental data.

Besides that, the collected comparisons allow for further emphasizing the key role of multiple input parameters on the reliable assessment of in-service glass structures, and, in particular, the role of restraints (compared to C-C or S-S conditions) and the presence of possible delaminations, which should be properly taken into account.

**Figure 21.** Comparison of experimental and analytical frequencies for the tested LG beams, as a function of deamination severity and interlayer stiffness (selection).

#### **5. Conclusions**

Laminated glass (LG) elements are largely used in buildings and civil engineering infrastructures, in the form of simple members (i.e., columns, beams and plates), but also combined together to realize stand-alone assemblies, where glass can interact with other traditional constructional materials. Besides the need of safe design methods for glass structures under ordinary loads, special care is increasingly spent by researchers in the dynamic response of LG systems under impact, moving loads due to pedestrians (in the case of walkways and roofs) or seismic events. As a part of buildings or complex systems, the dynamic parameters of LG elements should be properly taken into account, for the design of independent members or complex systems.

In this regard, several analytical methods available in the literature, for practical estimates of the fundamental vibration frequencies of simple LG members with viscoelastic interlayers. Otherwise, in most of the cases, reliable dynamic predictions can be obtained for the early-design stage. The intrinsic limit of literature methods is, in fact, represented by the assumption of ideal theoretical boundaries that often do not capture the mechanical effect of typical restraints in in use for LG systems. At the same time, they do not include the potential effects due to material degradations that can often manifest in in-service LG structures, like delaminations of the bonding interlayers.

In this paper, major issues for the vibration analysis of in-service (double) LG members were discussed and explored, with the support of analytical calculations, on-site experimental tests, and parametric Finite Element (FE) numerical analyses. As shown—as compared to laboratory research studies or early-stage design calculations—major uncertainties can be represented by the lack of accurate input material properties, especially for characterizing the shear stiffness of the bonding interlayers. The latter is strictly affected by the vibration frequency, but also by operational and

ambient conditions to properly assess. Such an issue can in fact reflect on unsafe equivalent thickness assumptions for analytical and numerical calculations, thus dynamic estimates. Further relevant influencing parameters that can magnify the material uncertainties and approximate assumptions are then represented by the actual boundary condition of real LG structures, which can be markedly different with respect to the ideal restraints of practical use for design. Accordingly, it was shown that even minor flexibility contributions of joints can turn out in marked stiffness for the examined LG members, thus requiring separate calculations. The restraints themselves can, in fact, modify the vibration frequency of a given LG members, thus further affecting the PVB stiffness and require iterative calculations. Finally, another relevant issue for in-service LG members derives from the degradation of common interlayers in use, especially delaminations that can be hardly quantified, but having marked effects on dynamic estimates.

Based on on-site experimental tests on in-service LG beams, these influencing parameters were separately explored in the paper. Analytical methods were then proposed for reliable estimates on general LG beams. Their advantage is that both restraint features and delaminations can be taken into account when compared to classical theories for slender composite beams of literature proposals for LG members in ideal conditions.

**Author Contributions:** The paper results from original research investigations carried out by the author.

**Funding:** No external financial funding was received for the study.

**Acknowledgments:** The So.Co.Ba Foundation "Società per la conservazione della Basilica" (Bergamin) is gratefully acknowledged for facilitating the on-site experimental measurements.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


© 2019 by the author. 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 (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
