**Seismic Response of Steel Moment Frames (SMFs) Considering Simultaneous Excitations of Vertical and Horizontal Components, Including Fling-Step Ground Motions**

#### **Shahrokh Shahbazi 1, Armin Karami 2, Jong Wan Hu 3,4,\* and Iman Mansouri <sup>5</sup>**


Received: 16 March 2019; Accepted: 14 May 2019; Published: 20 May 2019

**Abstract:** Near-field (NF) earthquakes have drawn considerable attention from earthquake and structural engineers. In the field of earthquake engineering, numerous studies have identified the devastating nature of such earthquakes, and examined the characteristics related to the response of engineering structures to these types of earthquakes. Herein, special steel moment frames (SMFs) of three-, five-, and eight-story buildings have been examined via a nonlinear time history analysis in OpenSees software. The behavioral seismic differences of these frames have been evaluated in two states: (1) under the simultaneous excitation of the horizontal and vertical constituents of near-field earthquakes that have Fling-steps in their records; and (2) under simultaneous excitation of the horizontal and vertical constituents of far-field (FF) earthquakes. In addition, during modeling, the effects of panel zones have been considered. Considering that the simultaneous effects of the horizontal and vertical constituents of near-field earthquakes were subjected to a fling-step resulting in an increased inter-story drift ratio, the horizontal displacement of stories, an axial force of columns, created the moment in columns, base shearing of the structure, and velocity and acceleration of the stories.

**Keywords:** near-field earthquake; fling-step; far-field; simultaneous excitation; special moment frame (SMF)

#### **1. Introduction**

According to the elastic rebound theory, the sudden release of gradually accumulated stress and displacement in the earth's crust leads to the formation of a large permanent ground displacement over a few seconds and creates a fling-step [1] (Figure 1). The result would be a long-period pulse which is different from the better-known directivity pulse produced by the constructive interference of propagating seismic waves. Ventura et al. [2], and also Grimaz and Malisan [3], investigated the influences of fling-step on buildings through a parametric study; their study proves that the motion with the fling effect creates a much larger response than those without the fling effect (and thus potentially greater damage). Moreover, several studies have shown a considerable difference between the seismic demands of near- and far-field records. This is due to the high frequency content of near-field records and pulse-like behavior of these types of ground motions [4–8].

**Figure 1.** Ground motion depicting fling-step (i.e., static offset) from the YPT-N/S station in the 1999 Kocaeli, Turkey earthquake. Dp: Displacement amplitude; Tp: Period or duration; T1: arrival time.

The main purpose of this paper is to develop a study of reference [8]. Kalkan and Kunnath [8] considered only the horizontal component of ground motion, while in this paper the simultaneous excitations of the horizontal and vertical constituents are considered. Besides, the panel zones are modeled in this paper, which were ignored in reference [8]. Among the seismic demands, only inter-story drift has been selected in the study of Kalkan and Kunnath, however in the present paper, in addition to inter-story drift, story displacement, columns' axial force, columns' moment, story shear, velocity and acceleration of story through an extensive post-processing are investigated.

Numerous studies have sought to clarify the response of buildings to diverse kinds of abbreviated pulses and create a design method for pulse-type records. For instance, references [9–11] evaluated the response of buildings and predicted the near-field response spectrum. References [12–17] examined the effects of pulse-like motions on building response. Moreover, the impacts on strength reduction factors were assessed in references [9,15,18,19]. Recently, Shahbazi et al. [14] studied the effects of soil classification on the seismic behavior of steel moment frames (SMFs), considering soil structure interactions and near-field ground motions.

An exclusive sine wave model for fling-step pulse waveforms was recently introduced by Burks and Baker [20]. However, a single sine wave cannot fully describe fling-step pulses, as they may consist of more than one frequency. Burks and Baker [20] presented a number of equations to predict the parameters of their model without the need for producing fling-step pulses.

Since these techniques take a lot of time, Farid Ghahari et al. [21] extracted the pulse component of a ground motion by adopting the moving average method. Soil structure interactions (SSIs) can also affect the response of buildings. The seismic evaluation of SMFs subjected to near-field ground motions with a forward directivity, including SSI impacts, were investigated by Shahbazi et al. [22], and the effect of fling-steps were ignored.

Moment frame systems are simply defined as assemblages of rigidly connected columns and beams [23–26], with its lateral strength mainly being achieved by the creation of bending moments and shear forces in the members and joints [27–30]. Owing to their pulse-like properties, near-field earthquakes result in a fast transmission of energy from the ground to the structure. Therefore, the elements in charge of increasing structural ductility fail to optimally perform. Thus, this study selected special steel moment frames, which have high ductility expectations. Moreover, in order to consider more realistic conditions in the modeling process, not only were the horizontal and vertical effects of near- and far-fault records simultaneously applied to the frames (which is the case in real earthquakes), but the effects of panel zone modeling were also incorporated.

#### **2. The E**ff**ects and Characteristics of Near-Field Earthquakes**

Past studies show that near-field earthquakes have the critical energy of pulses. However, while these earthquakes may have a smaller magnitude than Richter, they have a high potential for damage. In general, near-field earthquakes have specific characteristics, which are briefly mentioned as follows:


These characteristics arise from the many phenomena that occur near the seismic zone, which will be discussed later.

#### *Fling-Step E*ff*ect*

One of the recorded characteristics of the near-field records in last two decades, such as the 1999 earthquakes in Turkey and Taiwan, is the fling-step. This property occurred in the range of a few seconds of slip in the direction of the fault slip and is independent of the dynamic displacement originating from the directivity pulse of the ground's rupture. Accordingly, this phenomenon will affect components that are in the line of strike (parallel to the slip fault in strike-slip earthquakes, and in the direction of the slope in the dip-slip earthquakes).

For dip-slip earthquakes, fling-step and directivity effects occur in one direction and fling-step ground simultaneously occurs with maximum dynamic displacement, which must be considered as loads perpendicular to each other. The impact of the loads at one time is associated with a high potential for damage [32]. Briefly, fling-step ground is associated with a long-range velocity pulse and a uniform step in the time history of displacement. For instance, the Yarimka station experienced an approximately 2-meter permanent displacement during the Kjaili earthquake [33] (see Figure 2).

**Figure 2.** NS direction Yarimca station recordings during the Kocaeli earthquake [33].

The below figure is schematic of a strike-slip fault, the pulse is caused by the forward orientation and subsequent deformation which occurs in the perpendicular parameter and parallel to the major parameter of the fault, respectively. However, these two impacts are considered in the perpendicular parameters to the fault for dip-slip faults (see Figure 3), and time intervals where these factors are

demonstrated with each other and separately are represented in Figure 4. Figure 3 displays the trend of pulse movement along with remnant displacement for strike-slip and dip-slip faults. The time interval plots caused by strike and reverse faults can be observed in Figure 4.

**Figure 3.** Conceptual plot representing the trends of the fling-step and directivity pulse for strike-slip and dip-slip faulting [32].

**Figure 4.** Conceptual plot showing the time intervals of strike-slip and dip-slip faulting in which the fling-step and directivity pulse are displayed together and discretely [32].

#### **3. Numerical Models**

#### *3.1. Selected Numerical Models*

In the current research, steel moment skeletons of 3-, 5- and 8-floor with a particular ductility were applied for investigating the seismic behavior of three models founded on soil type 2 (375 m/s ≤ Vs (shear-wave velocity) ≤ 750 m/s). Seismic design of these three models was conducted using ETABS software and Iranian national construction codes [34]. The height of stories is 3.2 m and each of them was founded with three linear and lateral spans by the same width of 5 m [35]. In addition, corresponding columns and beams were chosen on both sides with stories undergoing similar efficiencies. Based on European standard cross sections, various kind of profiles were assumed for beams and columns, yielding cross sections which were employed in this research for beams and box-shapes were suggested for columns. As a result, the schematic model of the reference skeletons is illustrated in Figure 5. The profiles of the beams and columns for three-, five- and eight-story skeletons are presented in Table 1. Furthermore, leaning columns were applied to modeling second order impacts (Figure 5).

**Figure 5.** Topology and grouping details of stories [26].


The following symbol is applied for the columns place in the analysis result section (in Section 4 of paper):

Cij\* is the code for position of results in the columns;

i = story number;

j = number of columns located in the left-hand of buildings.

There are different assumptions being considered in this research. Dead (fixed) load on a building is due to the weight of constant elements including beams, floor segments, columns and walls. The fixed load is used as a monotonic distributed load parallel to each beam, column, or wall and as a surface load for the deck/segment. The magnitude of the beam and column loads equals the surface area times the material weight density. The magnitude of the wall load equals wall's thickness times the wall height times the material density. The live load is calculated based on the standard loading code. The live and live loads of 20 kN/m<sup>2</sup> and 65 kN/m2, respectively [36], were applied to all stories. But different loads were applied for roofs, at 54 kN/m2 and 15 kN/m2, respectively [26,37]. All models were fulfilled based on the rigid bed which was used for all stories. Elastic elements were considered for all beam and columns in OpenSees software which has been employed for modeling these structures. Bilin non-multi linear materials were used in order to describe the behavioral properties of elements [26]. The Krawinkler panel zone model was selected in this research (see Figure 6). The parameters of the Bilin material were obtained based on the study by Lignos and Krawinkler [38].

Panel zone distortion occurs mainly in shear because of contrary moments in the columns and beams. Panel zone modeling is explicitly performed through the method described by [39]. In this approach, modelling of the panel zone is performed in a rectangular pattern with eight very stiff elastic beam-column components. Shear deformation in the panel zone is shown by one rotational source (Figure 6).

**Figure 6.** Conceptual display of an ordinary panel zone [40].

The Bilin material is defined as the corrected deterioration model of Ibarra–Medina–Krawinkler with a bilinear disturbed response. Figure 7 presents the components of the Bilin material. Lignos and Krawinkler presented the relationships between the variables [38]. The 3-, 5- and 8-story buildings have fundamental periods equal to 0.48, 0.91, and 0.78 s, respectively.

**Figure 7.** The corrected deterioration model of Ibarra–Medina–Krawinkler [38].

For simulating the nonlinear behavior of the investigated structures, their modeling was performed by using elastic beam-column components in association with rotational sources. According to the corrected IMK (Ibarra–Medina–Krawinkler) deterioration approach, the sources are supposed to represent a bilinear hysteretic reaction. Rotational springs are employed to represent the nonlinear behavior of moment frames using the concentrated plasticity concept.

A rotational spring, located at the center of the reduced beam sections (RBS) and connected to the panel zone through an elastic beam-column member, is used to model the plastic hinge. Each frame member is modeled as an elastic element with two rotational springs at either end to connect it to other elastic elements. Therefore, in order to simulate the stiffness of the actual frame member, the stiffness of different components of the structure should be adjusted [41].

#### *3.2. Records Selection*

Choosing the appropriate record is one of the most important parts of time history analysis. In this research, 12 earthquake records have been selected for nonlinear time history analysis, meanwhile, near-field earthquakes have the effects of fling-step earthquakes. The characteristics of these seismographs have been presented in Table 2.

#### **4. Evaluating the Seismic Response of Structures**

The 5% damped response spectrum should be produced and combined for each couple by applying the square root of the sum of the square (SRSS) approach, so that a specific spectrum is built for each couple. Ground movements should be measured such that the mean value of their SRSS spectra does not decrease less than 1.4 times the standard design spectra during 0.2*T*–1.5*T* seconds, where *T* indicates [34].

The present research analyzed 36 nonlinear time histories according to 12 chosen ground motions. The maximum inter-story drift ratio was used as the first criterion for assessing the seismic demand on the structures. Figure 8 shows the maximum drift ratio stemming from the nonlinear time history analysis under the simultaneous effects of the horizontal as well as vertical components of ground motions for buildings with a special moment framework.



a Faulting Mechanism = TH: Thrust; REV: Reverse; SS: Strike-slip; OB: Oblique.

**Figure 8.** The maximum inter-story drift ratio of stories in studied frames under the simultaneous effects of the vertical and horizontal constituents of far- and near-field ground motions with fling-step records.

The graphs suggest that, in a three-story building, the maximum inter-story drift ratio was 0.03 under the simultaneous effects of the horizontal and vertical constituents of near-field earthquakes with fling-step and 0.01 for far-field earthquakes. These maximum ratios are associated with Big Bear (#(#EQF(V+H)1) and Chi-Chi-52 records (#EQN(V+H)2), and emerged on the second floor of the structure in both earthquakes. In addition, the maximum inter-story drift ratio was three times larger in the near-field earthquake compared to that in the far-field earthquake.

This maximum ratio was 0.01 subjected to the influence of far-field records and 0.06 subjected to the effect of near-field records in a five-story building, which were respectively associated with Chi-Chi-52 #EQN(V+H)2) and Kern County (#EQF(V+H)2). The inter-story drift ratio was maximized on the second floor of the structure subjected to the influence of Chi-Chi-52 (#EQN(V+H)2) and on the third floor in Kern County (#EQF(V+H)2). This maximum ratio was six times larger in the near-field earthquake compared to in the far-field earthquake. Although this trend is similar in both buildings, the inter-story drift ratio was maximized on the sixth floor of the five-story building in the near-field earthquake of Big Bear (#EQF(V+H)1) and on the fifth floor of the eight-story building in Chi-Chi-52 (#EQN(V+H)2). Table 3 compares the average maximum inter-story drift ratio in all three structures for far-field and near-field earthquakes. In this table, "Near/Far" denotes the ratio of near- and far-field.


**Table 3.** Comparison of the average of seismic response parameters of studied frames under to the simultaneous effects of the horizontal and vertical constituents of far- and near-field earthquakes with fling-step records.

In the three-story building, the maximum horizontal displacement of the roof under the simultaneous influence of the horizontal and vertical constituents was 0.28 m for the near-field earthquake and 0.10 m for the far-field earthquake. This displacement was caused by the Big Bear(#EQF(V+H)1) and the Chi-Chi-52 (#EQN(V+H)2) records. The displacement imposed on the structure by the near-field earthquake was 2.8 times greater than that caused by the far-field earthquake.

According to the Kern County (#EQF(V+H)2) and Chi-Chi-52(#EQN(V+H)2) records, in the five-story building, the maximum horizontal displacement of the roof was 0.18 m under the simultaneous effects of the horizontal and vertical constituents of the near-field earthquake and 0.14 m for the far-field earthquake with the fling-step earthquake. The horizontal displacement of the roof in the five-story building subject to the near-fault ground motions was therefore 5.86 times greater than that of the far-field earthquake.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-52 (#EQN(V+H)2) records, the maximum horizontal displacement of the eight-story building was 1.15 m when subjected to the influence of both components of the near earthquake and 0.21 m when subjected to the influence of the far-fault ground motions. The horizontal displacement of the roof in the near-field earthquake was therefore 5.48 times higher than that of the far-field earthquake. Table 3 compares the average maximum horizontal displacement subjected to influences of far-field and near-field earthquakes.

The maximum axial force developed in the C12 column on the first floor of the three-story building was 582.10 kN when subjected to the influences of both the horizontal and vertical constituents of the near-fault record with fling-step and 289.50 kN when subjected to the influences of the far-field earthquake (Table 1). The Big Bear (#EQF(V+H)1) and Chi-Chi-29 (#EQN(V+H)3) records respectively imposed these forces on the columns. The peak axial force created in the column subjected to the influence of both components of the near-field earthquake was therefore 2.01 times higher than that of the far-field earthquake.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-52 (#EQN(V+H)2) records, the peak axial force developed by the near-fault record in the C13 column of the first floor of the eight-story building was 3109.9 kN and that by the far-field earthquake was 1598.50 kN. The maximum force developed by the near-field earthquake in the column was therefore 1.95 times higher than that developed by the far-field earthquake. Table 3 compares the average maximum axial force in the column subjected to the influence of far-field and near-field ground motions.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-52 (#EQN(V+H)2) records, the maximum moment developed in the C11 column of the first floor of the three-story building subjected to the simultaneous effects of the horizontal and vertical constituents of the near-fault records with a fling-step earthquake was 613.60 kN·m, and that of the far-field earthquake was 265.60 kN·m. The maximum moment developed by the near-field earthquake was therefore 2.31 times higher than that of the far-field earthquake.

According to the Kern County (#EQF(V+H)2) and Chi-Chi-52 (#EQN(V+H)2) records, the maximum moment developed in the C13 column subjected to the effects of the horizontal and vertical constituents of the near-field earthquake was 1578.00 kN·m and that of the far-field earthquake was 522.50 kN·m. The moment created in the column by the near-field earthquake was therefore 3.02 times higher than that created by the far-field earthquake.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-52 (#EQN(V+H)2) records, the maximum moment developed in the C13 column of the first floor of the eight-story building by the near-field earthquake was 2316.90 kN·m and that developed by the far-field earthquake was 796.80 kN·m. The maximum axial force developed by the near-field earthquake was therefore 2.91 times higher than that developed by the far-field earthquake. Table 3 compares the average maximum moment of the column subjected to the influence of far-field and near-field records.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-52 (#EQN(V+H)2) records, the maximum story shear developed in the three-story building subjected to the influences of the near-field and far-field earthquakes was 3940.80 kN and 3519.30 kN, respectively. The maximum story shear developed in the three-story building by the near-field earthquake was therefore 12% higher than that developed by the far-field earthquake.

According to the Big Bear (#EQF(V+H)1 (and Chi-Chi-52 (#EQN(V+H)2) records, the maximum story shear developed in the five-story building was 11,727.70 kN when subjected to the near-field earthquake and 9929.40 kN when subjected to the far-field earthquake. The maximum story shear developed by the near-field earthquake was therefore 1.18 times higher than that developed by the far-field earthquake.

According to the Big Bear (#EQF(V+H)1) and Chi-Chi-29 (#EQF(V+H)3) records, the maximum story shear created in the eight-story building by the near-field earthquake was 42,066.80 kN, and by the far-field earthquake 27,993.00 kN. The maximum story shear created in the eight-story building subjected to the influence of the near-field earthquake was 1.5 times higher than that of the far-field earthquake. Table 4 compares the average maximum story shear subjected to the influence of far-field and near-field earthquakes.


**Table 4.** Comparison of the mean story shear created in the stories of studied frames subjected to the simultaneous effects of the horizontal and vertical constituents of both far- and near-fault records with fling-step ground motions.

According to research and the Big Bear (#EQF(V+H)1) and Chi-Chi-29 (#EQN(V+H)3) records, the maximum acceleration developed on the floors of the three-story building under the simultaneous effects of the horizontal and vertical constituents of the near-field earthquake was 28.38 m/s2 and 10.81 m/s<sup>2</sup> for the far-field earthquake. The maximum acceleration created by the near-field earthquake was therefore 2.63 times higher than that created by the far-field earthquake.

In the five-story building, the maximum acceleration developed by the near-field earthquake of Big Bear (#EQF(V+H)1) was 19.49 m/s2 and that by the far-field earthquake of Chi-Chi-29 (#EQN(V+H)3) was 12.74 m/s2. The maximum acceleration was developed on the fourth floor in the near-field earthquake and on the roof of the building in the far-field earthquake. The maximum acceleration developed by the near-field earthquake was therefore 1.53 higher than that developed by the far-field earthquake.

The maximum acceleration developed on the roof of the eight-story building by the near-field earthquake of Big Bear (#EQF(V+H)1) was 23.34 m/s2, and that developed by the far-field earthquake of Chi-Chi-52 (#EQN(V+H)2) was 12.02 m/s2. The maximum acceleration developed by the near-field earthquake was therefore 1.94 times higher than that developed by the far-field earthquake. Table 3 compares the average maximum acceleration subjected to the influence of far-field and near-field earthquakes.

The maximum velocity developed on the floors of the three-story building was 1.57 m/s under the simultaneous effects of the horizontal and vertical constituents of the near-field earthquake and 0.91 m/s by the far-field earthquake, which were imposed on the roof of the structure by the records of Big Bear (#EQN(V+H)1) and Chi-Chi-84 (#EQN(V+H)1). The maximum velocity caused by the near-field earthquake was therefore 1.73 times higher than that caused by the far-field earthquake.

The maximum velocity developed on the roof of the five-story building by the near-field earthquake of Kern County (#EQF(V+H)2) was 2.50 m/s and by the far-field earthquake of Chi-Chi-52 (#EQN(V+H)2) was 1.18 m/s. The maximum velocity caused by the near-field earthquake was therefore 2.12 times higher than that caused by the far-field earthquake.

The maximum velocity developed on the roof of the eight-story building by the near-field earthquake of Big Bear (#EQF(V+H)1) was 3.56 m/s and that developed by the far-field earthquake of Chi-Chi-52 (#EQN(V+H)2) was 1.40 m/s. The maximum velocity caused by the near-field earthquake was therefore 2.54 times higher than that caused by the far-field earthquake. Table 3 compares the average maximum horizontal displacement subjected to influence of far-field and near-field earthquakes.

Moreover, the fling-step versus near-field no fling-step motions is investigated. The results show that the average of maximum inter-story drift, story displacement, axial force of columns, end moment of columns, velocity and acceleration of the story under fling-step ground motions are less than those with near-field no fling-step motions (Table 5). However, Table 6 reveals that the story shear is higher in the case of fling-step records.


**Table 5.** Comparison of the average of seismic response parameters of the studied frames under the simultaneous effects of the horizontal and vertical constituents of near-field earthquakes with fling-step and non-fling-step records.

**Table 6.** Comparison of the mean story shear created in the stories of studied frames subjected to the simultaneous effects of the horizontal and vertical constituents of both near-fault records with fling-step and non-fling-step ground motions.


By comparing the below table results, it can be said that the response of structural excitation is subjected to the influence of the horizontal component of the near-field earthquake, given the permanent displacement of the earthquake is greater than the response of the structure subjected to the influence of the horizontal component of the far-field. This conclusion was achieved from evaluating structural excitation under the simultaneous effects of the horizontal and vertical components of both far- and near-field earthquakes (see Tables 7 and 8).


**Table 7.**

Comparison

 of the maximum means of the inter-story drift ratio, horizontal

displacement,

 acceleration

 and velocity in studied frames subjected to just the

*Appl. Sci.* **2019** , *9*, 2079


*Appl. Sci.* **2019** , *9*, 2079

**Table 8.**

Comparison

 of the maximum mean of the inter-story drift ratio, horizontal

displacement,

 acceleration

 and velocity in studied frames subjected to

#### **5. Conclusions**

The fling-step property due to a stationary ground displacement is generally performed through a combination of a unilateral velocity pulse and a uniform stage in the displacement time interval. The separate stage in the displacement time interval happens along the fault slip (i.e., in the direction of strike and dip for strike- and dip-slip occurrences, respectively).

According to existing studies, near-field earthquakes are characterized differently from far-field earthquakes in terms of the frequency content and amplification range. Thus, it is necessary to study the effects of such parameters on structures. With respect to near-fault fields, forward directivity and fling-step pulses may appear in the case of each component with a parallel or perpendicular position to the fault line, either separately or together in one record component. This study evaluates the seismic behaviors of three-floor, five-floor and eight-floor special steel bending frames that undergo near-field earthquakes with fling-step effects and also far-field earthquakes. This study also modeled the panel zone in order to render a more realistic behavior of the structure. Non-linear time-history analysis of the structures demonstrated that 'floor displacement', 'inter-story drift ratio', 'floor acceleration', 'axial and anchor force on the columns' and 'floor velocity' subjected to influence of the horizontal component of near-field earthquakes with a fling-step effect are greater than the corresponding values of these parameters in far-field earthquakes.

**Author Contributions:** S.S. and I.M. figured out the given idea. S.S. and A.K. presented the theory and conducted the computations. I.M. and J.W.H. investigated and confirmed the analytical approaches. J.W.H. and I.M. persuaded S.S. and A.K. to study both orientation and fling-step impacts and controlled the results of this research. All researchers considered the findings and provided the ultimate work.

**Funding:** Basic Science Research Program of the National Research Foundation of Korea (NRF) provided facilities for this research, for which the authors are grateful. Financial support of the Ministry of Science, ICT & Future Planning (2017R1A2B2010120) is highly appreciated.

**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* **Modeling Palletized Products: The Case of Semi-Filled Bottles under Top-Load Conditions**

**Ana Pavlovic 1, Cristiano Fragassa 1,\*, Luca Vegliò 1, Felipe Vannucchi de Camargo 1,2 and Giangiacomo Minak <sup>1</sup>**


Received: 21 August 2019; Accepted: 7 October 2019; Published: 2 January 2020

**Abstract:** An investigation on numerical methods able to simplify the mechanical behavior of PET bottles, partially filled with liquid, under compression loadings is presented here. Compressive stress conditions on bottles are very common during their transportation and can be accompanied by large deformations and instabilities that can compromise the integrity of the pack, with the risk of significant damages. The present paper proposes two approaches, both based on finite elements, and with elastic-plastic material models properly defined with the scope of investigating the complex phenomena that take place during these loading conditions. Although not perfect in terms of accuracy, these numerical methods have been proven to be capable to predicti the transport-related integrity risks, showing results that agree with experimental data, especially during the initial phases of load compression.

**Keywords:** transport; palletized goods; damage; bottle; buckling; Polyethylene terephthalate (PET)

#### **1. Introduction**

The irreversible damage of goods during transportation, caused mainly by sloshing, vibration, and shock mechanical solicitations [1], is a common scenario nowadays in the logistics sector, resulting in significant economic losses to both the retailer and manufacturer [2]. The environmental aspect is particularly affected by the greenhouse gas (GHG) emissions to produce and ship goods that, being damaged, will not satisfy the customer and will end up in a landfill. On the other hand, the economic losses can be ultimately managed in some way. For instance, an increase in the price of goods, as estimated with proper methods, as in [3], can be added to ensure the economic sustainability of the industrial investment. On the contrary, the environmental cost inferred by deteriorated goods, as claimed in [4], is often definitive and articulated on multifaceted aspects, such as materials no longer usable, energy and water used for product processing, and fuel and carbon dioxide used during the transport phases.

The loss in functionality and usability of goods during transportation can be considered a non-negligible fraction of the whole number of products moved around the world. For instance, a study commissioned by Walmart [5] estimated the phenomenon in 1,176,260 pallets damaged each year, equivalent to 58,813 truckloads of product going into landfill that, at the financial level, equates to about \$2400M and 18% earnings. Many leading companies are trying to identify and stem the causes of the deterioration of goods due to transport [6]. The evaluation of the economic repercussions of these phenomena through in-depth analyses [7] can lead to the identification of solutions to be adopted by the company to minimize losses.

Several recent studies [8,9] have been directed to optimize transport conditions, developing solutions for what is commonly referred to as 'secondary packaging', which does not take into consideration packaging materials that are not in direct contact with the product (e.g., blisters, sachets, boxes). In particular, in the case of bottles, the secondary packaging refers to the action of arranging the bottles together into bundles covered with stretch wrap films, positioned side by side and overlapped to form the pallets (Figure 1a). The pallets, which, in the case of bottles, can also contain tens of bundles arranged on several levels, are protected by the deposition of several layers of stretch film. This plastic and thin film is deposited with special implants, combining both the wrapping and pre-tensioning actions of the film itself (Figure 1b).

**Figure 1.** A multilevel system: (**a**) a single bottle; (**b**) bottles grouped in bundles; (**c**) pallet and palletizing plant.

In this way, a stabilization effect is created with respect to external dynamic stresses, which is added to the general protection offered by the presence of the plastic layer outside. Through a suitable choice of the type of film and the wrapping parameters it is possible to guarantee suitable levels of protection with respect to the most common stress conditions [10,11] that are generated during the transport phases. This is possible as long as the bottles retain their structural integrity. The external packaging can do very little if the stressors on the individual bottles are such that can damage them. This effect is not so rare, considering the attempt to use progressively lighter and slimmer bottles close to the limit conditions, in the face of always large and heavy pallets. The result is that bottles are subjected to various types of loads during transport and with various possibilities of damage, with one of the most dangerous being that which leads to buckling.

Plastic bottles have been widely used since 1989 as containers for water, milk, and other beverages. The body is usually made of PET (polyethylene terephthalate) from pre-printed plastics (Figure 2a) and then blown, while the bottle cap is made of PE (polyethylene). The form is not only conceived for an aesthetic purpose, but it also aims to guarantee adequate stiffness during transport (being stacked on top of each other) and specific features during use (Figure 2b).

**Figure 2.** (**a**) Preform in PET (polyethylene terephthalate), from which the bottle is obtained through a thermoforming process; (**b**) example of the final geometric complexity of a bottle.

The understanding of what happens during the crushing of the bottle involves a number of complications, and, to detail this complexity, it must be considered that plastic bottles are usually:


The bottles are then filled, not completely, with the liquid to be contained and transported (natural water in the present case). By its turn, the partial filling represents one of the most significant reasons able to add complexity to the problem [12]. The air cushion inside the bottle creates a gas–liquid–structure interaction in response to any external mechanical (or thermal) stresses. This is the main reason why assembling a numerical model becomes challenging.

Several studies are available where experimental tests were conducted with methods and results similar to those presented in this study. For instance, in [13,14], the bottle was characterized by a monoaxial compression test, extrapolating the trend of the force as a function of the lowering. In [15], a system was specifically developed for the evaluation of the contact pressures between the bottom of the container and the fixed base, using pressure transducers in order to be experimentally correlated with a drop test. In [16], a probe for monitoring the internal pressure of the bottle was considered, stating the decrease in the volume of air during the deformation phase. In [17], a system of circumferential bands driven by a traction system was described to determine the radial resistance of the component, being able to characterize the body also orthotopically.

It must be also considered that there are widespread standard tests, for example, top-load tests (according to ASTM D 2659-95 [18], as in [17,19,20]), that lead this kind of test as a common practice, also in the case of industrial applications. The difficulty is not related to the implementation of these tests or in collecting pertinent experimental data, but in the possibility of generating a complete model (theoretical or numerical) able to interpret the information.

There are several studies that attempt to tackle this theoretical or numerical modeling through two different approaches. The first type of method, closely linked to the bottle-producing companies, consists of modeling the body according to the geometry, type, and quantity of material necessary to guarantee a predetermined resistance. In numerous studies (e.g., [21,22]) the simulations were totally focused on the minimization of the thickness to reduce the cost or, as analyzed in [23], to introduce new types of bottles to the market. Once implemented, the characteristics of the material and the discretization of the profiles with a small variable thickness represent a computational challenge, which can be addressed with the models presented in [24–26].

The second method consists of evaluating the interaction between the liquid, gas, and wall, described in [27]. This type of approach, based on both the Eulerian and Lagrangian models, is used to define the volume of liquid and air present inside the bottle, losing precision in the definition of the container itself. An example of this approach can be found in [28].

The various papers presented in [19,23] on static and dynamic simulations describe the numerical methods used for the simulation of the bottles. The effect of the interaction between water and air analyzed in [23,29] was considered only from a phenomenological point of view, that is to say, by representing it by a pressure acting on the inner surface of the container, as in [16,30,31]. In some cases, this phenomenon was modeled by means of a membrane consisting of elements that can be infinitely extended and which simulate a contact layer between the contents of the bottle and the wall. Unfortunately, such an approach, as used in [17,19], is not described in detail, making this type of simplification useless. Thus, it is possible to say that, at the moment, there are no studies that consider all the different phenomena of interest altogether. Furthermore, the complexity that has emerged would produce formulations and models that are difficult to apply in an industrial context, whereas the present interest relies on having rapid information on the risk of damage compared to the variation

in shape, thickness, material, and stressors. All this seems possible only through the realization of specific experimental campaigns, one for each of the cases of interest, which represents an expensive and not always feasible approach.

#### **2. Materials and Methods**

#### *2.1. Aim and Scope*

The present study aimed to analyze the behavior of a thin plastic bottle when subjected to vertical compression loads by integrating a numerical approach and an experimental one. For this purpose, the study presented a needful simplification of the phenomena acting on the bottle [32], characterized by a geometrically complex profile and made of PET, with the orthotropic characteristics given in Table 1. The same table also reports the main isotropic properties of PE, used for the bottle cap.


**Table 1.** Material characteristics of PET and polyethylene.

The bottle was considered to be partially filled with liquid, closed by the cap, constrained in the base by a friction contact, and subjected to a compression load along the longitudinal axis. Every study took into consideration:


**Figure 3.** Modeling the bottle: (**a**) main elements; (**b**) example of characteristic thickness [33].

With this scope, two approaches toward the simplification are proposed here.

#### *2.2. Dual Approach*

The first approach involved the use of an exactly faithful external geometry, as obtained from the 3D model supplied by the bottle manufacturer. In the case this geometry is not available, it could be obtained by reverse engineering, while the tendency of the thickness could be derived by analogy with respect to similar cases, available in the literature (Figure 3b). In terms of numerical discretization, the bottle was approximated by two-dimensional shell elements. The material properties were found in the vast specialized literature, taking into account the general characteristics of PET and PE. In particular, in the present case a bilinear isotropic hardening model was used for modeling PET (in accordance with [19,23,34]). The thin walls of the bottle were subjected to compression stress, to which a pressure component linked to the progressive compression of the air gap was added. The trend over time of the pressure during the test could be estimated through the equation of real gases (or, given the low pressures involved, through the perfect gas equation). However, according to the literature, which foresees the explosion of the bottle in the presence of internal overpressures around 1 MPa in the case of 1.5 L bottles, a linear increase was used up to this value (Figure 4) [23]. The liquid was considered incompressible and capable both of transferring the pressure exerted by the gas to the contact walls and of adding the effect of its own pressure due to the water column (Stevino's law).

**Figure 4.** Bilinear isotropic hardening material model (in the case of 75 MPa for the yield strength and 1 MPa for the tangent modulus).

The second approach, of lower complexity, involved the use of a simplified bottle geometry, obtainable from the available databases, and chosen only by analogy with the bottles of real interest. This geometric model can be considered a single solid block, discretized by 3D solid elements (including air and liquid). The properties of the materials were assigned in such a way as to respect the macroscopic behavior of the bottle, as emerged from the experimental tests. To make this approach possible, different models and materials had to be considered and compared with the experimental evidence to verify when the top load tests were correctly reconstructed. Although this simplification may seem extreme in terms of simplification, compared to a plain investigation of the physical phenomena, sometimes it represents the only possibility. In fact, it made it possible to model the behavior of the pallet, of which the bottle represents only one of the elements. Once the model of the bottle was brought back with a black-box approach, to which numerical, literature, and experimental elements were integrated, it was possible to study the interaction between these black-boxes, grouped in bundles, which in turn interact in the composition of the pallet, all mediated through the action of polymeric coatings.

#### **3. Experiments**

#### *3.1. Experimental Evidence*

During the experimental tests, five 1.5 L water bottles (h = 310 mm, φ = 87 mm) were considered as specimens. Compressive tests were realized on an Instron 8032 servo-hydraulic universal testing machine with a 2517-005 load cell with ±250 kN dynamic and ±500 kN static range at a 0.25 mm/s ramp accounting for different displacements. The output data given by the testing machine were force, time, and displacement. Tests were run along the bottle axis (vertical direction), with a behaviour as shown in Figure 5. Force–displacement diagrams were also measured and displayed in Figure 6 where results from five tests are reported. These trends show an initial proportional zone (a), the bottleneck completely dented (b) and, lastly, the buckling phenomenon (c).

It can be noted that an exceptional event occurred at 8–10 mm of displacement, equal to around 2.5–3.0% of axial deformation and 0.5 kN, able to split the graph in two different regions, one with linear and the other with nonlinear behaviors. It corresponds to the indentation of the bottle, with the cap collapsing inside the neck (Figure 7a). Once the elastic response zone was passed, various phenomena became visible, increasing the load. In particular, over 2 kN, buckling happened in the specimens due to the geometry of the bottle: the tapered shape in the midways transversal plane prompted this behavior, whereas axial compression was left aside. At the end of the test, in the absence of a load, a residual deformation was also present (Figure 7b).

**Figure 5.** Predominance of buckling on single-bottle compression tests: initial part of the test with a proportional response (**a**), as soon as the bottleneck is completely dented (**b**), and near to the 3% strain region when buckling takes place (**c**).

**Figure 6.** Force–displacement diagrams of specimens: (**A**) whole trend and (**B**) details showing the proportional response zone (**a**), the bottleneck completely dented (**b**) and final buckling (**c**).

**Figure 7.** Residual deformations with tests ending before (**a**) and after buckling (**b**).

#### *3.2. An Interpretation*

Within the scope to anticipate the simplification proposed later, the single bottle can be treated as a column. In this case, from the experimental measures it was possible to build an equivalent stress–strain chart to figure out an equivalent elastic modulus. The force was transformed into stress by dividing it by the circular transversal area of the bottle, with a nominal diameter of 87 mm. On the other hand, the strain was calculated based on the initial length of the specimen, i.e., the height of the bottle, equal to 310 mm.

This stress–strain diagram is reported in Figure 8a, where it is possible to verify the same unexpected variation in the trend. Before this change, a clear elastic tendency is shown until at least 1% of deformation, in which the slope can be represented using the Young's modulus (similarly to the case of the elastic response of materials). Then, the elastic modulus could be also calculated as the slope of a highly correlated linear fitting curve from 0% to 1% (Figure 8b). In other terms, if it was possible to consider the system of a bottle as an isotropic material, then the Hooke's formula should be taken into account, yielding a value equal to:

$$E = \text{--} \pounds 08 \pm 0.47 \,\text{GPa.} \tag{1}$$

**Figure 8.** (**a**) Experimental stress–strain curve; (**b**) calculation of elastic modulus on the linear zone.

Analogously, treating the single bottle as a column, the theoretical threshold force for buckling can be calculated. This theoretical analysis, given in Equations (2)–(4), assumes the moment of inertia *I* (mm4) of a solid circular cylinder (because after the bottleneck is dented the fluid should be considered as incompressible), the elastic modulus *E* (MPa), the column length *L* (mm) taken as 310 mm, and finally the buckling factor *n* which is equal to 1 in this condition where both ends are pivoted, and not fixed.

$$F = \frac{n\pi^2 EI}{L^2} \,\prime\,\tag{2}$$

where:

$$I = \frac{\pi d^4}{64},\tag{3}$$

so,

$$F\_{buck,theory} \cong \text{750 N.}\tag{4}$$

At the same time, it is evident the sharp decrease in the strain deformation diagram (shown by the transition line in Figure 8a) can be attributed to the variation in the resistant section during compression, coinciding with the collapse of the bottleneck (Figure 5b), before reaching the column of liquid. During the first compression phase, the force acts on the section corresponding to the cap, as the bottle deforms longitudinally, the resistant section widens until it reaches the maximum resistant section, which coincides with the meeting of the fluid column, equal to 5942 mm2.

#### **4. Numerical Models**

#### *4.1. Shell (2D) Model*

Two-dimensional shell elements were used in this first finite elements (FE) study to evaluate the behavior of PET bottles under compression in the case of top-load tests. In particular, both a non-linear static and a buckling analysis were done by ANSYS Workbench v.18.1.

The geometry was simplified as displayed in (Figure 9a): all secondary details like fittings, chamfers, and notches were removed, since they marginally affect the system; the cap was merged to exclude contact problems; the thickness was simply defined as uniform and equal to 0.3 mm.

Four-node shell elements, 13,258 in total, were used to create the FE grid. This mesh, optimized for buckling analysis, can be observed in the area of cap and neck in Figure 9b.

In the base of the bottle all the translations and rotations were constrained (Figure 9b). The effect of water and air inside the bottle was replaced by an internal pressure, increasing over time roughly in accordance with the equation of state of perfect gases. The contact between bottle and plate was defined as a contact with friction, with a coefficient of friction equal to 0.2 representing the friction between PET and steel. Material properties were defined in line with Table 1 with a bilinear isotropic hardening model with a plastic zone starting from a yield value of 75 MPa (Figure 4).

**Figure 9.** (**a**) Bottle geometry before and after simplification; (**b**) detail of geometry and mesh by shell elements; (**c**) application of forces and boundary conditions (case of non-linear static analysis).

With the scope to perform the quasi-static analysis, the bottle was loaded by a rigid plate moved by a displacement of 10 mm in the vertical (Y) direction (Figure 9c). This value was chosen by considering the maximum displacement in the elastic region as measured by top-load tests (Figure 6). In some ways, it represents the experimental evidence of the maximum load that this specific bottle can resist without moving its status to the plasticity.

The simulation of maximum load is widely used in the PET bottle industry and the vertical load test determines the weak points inside the bottles. Thus, the numerical analysis limited to this value can be extremely useful, for example, in detecting problems in the blow molding line and to check if the component failures are repeated in the same points.

Also, in this case, a non-linear structural simulation was performed with the maximum load for bottle verification in order to obtain a reliable model of an empty PET bottle. The simulation was divided into two steps: load in terms of displacement up to 10 mm, and pressure increase up to 1 MPa. Non-linearity due to geometry was taken into account through large displacements.

The stress versus strain status following a 10 mm displacement is shown in Figure 10 in terms of (equivalent) von-Mises stress (Figure 10a) and total deformation (Figure 10b). In particular, the maximum stress was 164 MPa. This stress state exceeds the yield point, showing the bottle in plastic conditions.

**Figure 10.** Quasi-static analysis in top-load conditions for a 10 mm vertical displacement: diagram of stress (**a**) and strain (**b**).

The buckling analysis was performed using the same conditions of simulation. It made it possible to detect the geometrical configuration of instability together with the related buckling load factor (BLF). This factor acted as a load multiplier with the scope to evaluate the buckling load for each condition of instability. When the BLF is greater than one, the structure is safe apart from when the external load is increased. It also means that when the buckling factor is less than one, the structure could collapse under a given load.

This consideration is exactly in line with what is highlighted by the results shown in Figure 11. With BLF < 1 in respect to several configurations, the thin bottle is not able to resist to a 10 mm displacement without buckling. At the same time, as reported by the simulation, the deformations seem very limited, suggesting the buckling is locally restricted.

**Figure 11.** Buckling analysis in top-load conditions for a 10 mm vertical displacement, limited to the initial four modes of instability and for a load multiplier (*k*) of: (**a**) *k* = 0.19; (**b**) *k* = 0.16; (**c**) *k* = 0.12; (**d**) *k* = 0.09.

#### *4.2. Solid (3D) Model*

Three-dimensional solid elements were used to evaluate an extreme hypothesis of simplification consisting of a bottle discretized by a homogeneous object, made of a single material. The presence of liquid, air, and the plastic polymers of small thickness were not taken into account separately, but the bottle was analyzed globally, as a single solid body. This body was considered as made by an equivalent material, whose characteristics were set in order to guarantee a mechanical response in line with the experimental tests. From this perspective, the geometry should be also simplified. Starting from the 2D profile, a 3D revolution solid was generated and discretized by solid elements (Figure 12a).

**Figure 12.** Finite Element modeling by solid elements: (**a**) 3D bottle model; (**b**) tetrahedral mesh; (**c**) boundary conditions.

A critical requirement for efficient and accurate explicit analysis is a high-quality mesh, to be also balanced with the availability of computational resources. Explicit analysis benefitted from the integration of the design phase in the ANSYS Workbench environment, which includes powerful automatic mesh generation. In general, there are two options, one is to create high-quality hex (brick) elements using multi-zone meshing, while the second one automatically uses the swept method to create hex elements. Some geometries, especially those created for manufacturing and imported from CAD data, are too complex to be swept to produce a full hex mesh. Tetrahedral elements can accurately represent small portions of a part's geometry that cannot be swept. On the contrary, tetrahedral elements can be not accurate enough in investigating large deformations and, in addition, mesh size can affect stress results. At least, different mesh density must be tuned to confirm the repeatability of this model. For this reason, it would be advisable to optimize the validity of the mesh, both in terms of element types and dimensions. At the same time, the ultimate goal of the research was to try to represent the behavior of a bottle inside a pallet in the most elementary way possible. For this reason, it was decided to limit the discretization to an automatic meshing by ANSYS tools.

In particular, quadrilateral dominant elements were used and the final mesh for explicit analysis consisted of 14,331 nodes and 13,275 elements (Figure 12a).

With respect to the first case, similarities and differences existed in the way loads and constraints were applied. The rigid plate on the top (Figure 12b) was set with a 2 mm/s speed (instead of a displacement) while bottle nodes on the bottom were fixed avoiding translations and rotations.

In an extreme simplification of this kind, which considers the entire system as homogeneous and discretized by finite elements of the same type, would be possible when a model of 'equivalent material' can offer the same responses as the real system (Figure 13).

**Figure 13.** Experimental stress–strain curve discretization by linear trends, namely bilinear isotropic hardening.

Ansys Workbench 18.1 also makes it possible to insert the material behavior by additional complex methods, as in the case of the multilinear isotropic hardening or Chaboche kinematic hardening models (Figure 14). These hardening models can be generated by fitting the experimental data and are able to provide rather different responses to the simulation. In Table 2 the accuracy of predictions is evaluated by comparing real and simulated strains in the case of different hardening models.

**Table 2.** Comparing measured and simulated residual strains (ε) for different hardening models.


**Figure 14.** Discretization by advanced hardening models: (**a**) multilinear isotropic; (**b**) Chaboche kinematic.

In particular, it is shown that simulations implementing the Chaboche kinematic hardening model are able to offer a better approximation of the system (with estimated errors in strains around 7%), as is also evident in Figure 14b, where an asymmetrical effect is correctly detected despite other models.

In Figure 15 the effect of deformation of the bottle during the top-load test with 90 mm forced displacement is shown as a simulation that adopted Chaboche kinematic hardening in comparison with the reality.

**Figure 15.** Effect of a top-load test with 90 mm forced displacement: simulation (**a**) of deformations in the case of the use of a Chaboche kinematic hardening model compared with the reality (**b**).

#### **5. Discussion**

Disregarding the uncertainties, it is fair to assume that the experimental values of 550 N for the threshold buckling force and 3.1 MPa of elastic modulus can be used for simplified numerical simulations. Furthermore, it is also evident that both the shell model and the solid elements present interesting results, such as to make them useful for a simplified modeling of the physical system. This general aspect assumes a certain importance in consideration of the many variables involved, which make an exact solution practically unreachable. Regarding this overall complexity, it can be noted, for instance, that 188 different materials are listed in the category of unreinforced polyethylene terephthalate (PET), together with their slightly different properties to be considered.

At the same time, it should also be highlighted that the experimental and simulated values are actually closer to each other for the following reasons:


#### **6. Conclusions**

The mechanical response and the stability conditions of a single bottle under compression loads were considered in this paper, as the first fundamental step of a multilevel approach aiming at the complete discretization of a pallet behavior (Figure 16).

**Figure 16.** Discretization by a multilevel approach: (**a**) single bottle; (**b**) bottles in bundles; (**c**) pallet.

In terms of methodology, experimental measures were used to investigate two different models of discretization, each one able to take into account some essential aspects of the overall problem, simplifying the others. Their assumptions were validated through a comparison between the numerical predictions and the experimental evidence. In particular, a top load test was used on five partially filled bottles, which revealed a mechanical behavior articulated in different phases.

As general results, it is possible to say that both models acceptably predicted the initial phases of bottle deformation, characterized by an elastic response of the system. They also provided useful information regarding the effects of instability within an acceptable margin of error, making it possible to use these data for improving the bottle design.

These simplified models were also very efficient in terms of computation time, representing a first concrete step towards the development of a simplified model able to predict the static and dynamic behavior of an entire pallet during the transport phases, which represents the real challenge for the packaging industry.

**Author Contributions:** Conceptualization, C.F. and A.P.; methodology, A.P. and C.F.; software, L.V. and A.P.; validation, A.P. and C.F.; formal analysis, G.M.; investigation, F.V.d.C., L.V. and C.F.; resources, C.F.; data curation, A.P. and F.V.d.C.; writing—original draft preparation, C.F. and L.V.; writing—review and editing, C.F. and F.V.d.C.; visualization, A.P.; supervision, G.M.; project administration, G.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the European Union and Region Emilia Romagna inside the POR-FESR 2014–2020 action.

**Acknowledgments:** This research was implemented as part of the '*TechnoLab 4.0*' project, managed by AETNA Group S.p.a. (Italy), thanks to the financial support from the European Union and Region Emilia Romagna inside the POR-FESR 2014–2020 action.

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

#### **References**


© 2020 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/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18