NUMERICAL MODELING FOR THE SEISMIC RESPONSE OF CONCRETE TILT-UP BUILDINGS by XAVIER TELLIER B.Sc., Ecole Polytechnique (France), 2010 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Civil Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) January 2013 © Xavier Tellier, 2013 ii ABSTRACT Concrete tilt-up building is a prevalent construction technique used for industrial and commercial applications in North America. This construction technique offers many significant advantages over conventional cast-in-place construction. This includes the reduction in construction time and the amount of formworks. Despite the large array of buildings that has been constructed using such technique, the nonlinear behaviour of the concrete tilt-up buildings is still not well understood. The nonlinear behaviour of the concrete tilt-up building has been studied in this thesis. The nonlinear response of the concrete tilt-up building is largely affected by the nonlinear behaviour of the connectors between the panels and the slab, and between the panels. Past researches have been conducted to experimentally examine the nonlinear behaviour of the tilt-up panel connectors. The experimental results were used in this thesis to develop an empirical numerical model capable of reproducing the force-deformation response of the tilt-up connectors under combined axial and shear deformation. The numerical model takes the shear strength and stiffness degradation into account after axial cycles of inelastic deformation. A finite-element software was developed specifically to study the nonlinear static and dynamic behaviour of concrete tilt-up buildings. A typical tilt-up building designed in the study of Olund (2009) was modeled. Incremental dynamic analysis was performed using the developed finite element software to assess the seismic performance of the prototype tilt-up building. The results of the incremental dynamic analysis provided valuable information to understand the nonlinear behaviour of the concrete tilt-up buildings under seismic load. Detailed parametric studies were carried out to examine the nonlinear behaviour of tilt-up buildings. Parameters such as connector configurations; variation of the roof stiffness and strength; and coefficient of friction between the panels and slab were studied. iii TABLE OF CONTENTS ABSTRACT .......................................................................................... ii TABLE OF CONTENTS.................................................................... iii LIST OF TABLES.............................................................................. vii LIST OF FIGURES........................................................................... viii 1 – INTRODUCTION.......................................................................... 1 1.1 General description of tilt-up construction method ........................................................... 1 1.2 Work needed ...................................................................................................................... 2 1.3 Research goals and scope................................................................................................... 4 1.4 Thesis overview ................................................................................................................. 4 2 - LITERATURE REVIEW .............................................................. 6 2.1 Review of the current construction techniques and design approach ................................ 6 2.1.1 Roof diaphragm .......................................................................................................... 6 2.1.2 Wall panels.................................................................................................................. 7 2.1.3 Connectors ................................................................................................................ 12 2.2 Review of previous experimental tests conducted on connectors ................................... 19 2.2.1 Experimental testing by Lemieux, Sexsmith and Weiler (1998).............................. 20 2.2.2 Experimental testing by Devine (2009) .................................................................... 23 2.3 Previous research on roof diaphragm .............................................................................. 26 2.4 Review of previous analytical models ............................................................................. 28 2.4.1 Olund (2009)............................................................................................................. 28 iv 2.4.2 Devine (2009) ........................................................................................................... 31 2.4.3 Need for a new model ............................................................................................... 32 3 - NUMERICAL MODEL FOR THE EM5 CONNECTOR ....... 33 3.1 Summary of experimental tests used in this study........................................................... 33 3.2 General mechanical behaviour of the EM5 connector..................................................... 35 3.2.1 Monotonic and cyclic envelope ................................................................................ 35 3.2.2 Cyclic behaviour and damage................................................................................... 37 3.2.3 Coupling.................................................................................................................... 38 3.2.4 Symmetry.................................................................................................................. 38 3.3 Proposed analytical model ............................................................................................... 39 3.3.1 Mathematical rules for the uncoupled model............................................................ 39 3.3.2 Coefficients of the uncoupled model ........................................................................ 48 3.3.3 Model of coupling..................................................................................................... 48 3.3.4 Model of failure ........................................................................................................ 51 3.3.5 Model of paired connectors ...................................................................................... 53 3.4 Validation of the model ................................................................................................... 56 3.4.1 Shear behaviour ........................................................................................................ 56 3.4.2 Axial behaviour......................................................................................................... 58 3.4.3 Coupled behaviour .................................................................................................... 60 3.4.4 Comparison with model of Devine (2009) ............................................................... 62 3.4.5 Sensitivity ................................................................................................................. 63 v 4 - NUMERICAL MODELING OF THE SEISMIC BEHAVIOUR OF TILT-UP BUILDINGS................................................................ 64 4.1 Description of the finite element model........................................................................... 64 4.1.1 Modeling assumptions .............................................................................................. 65 4.1.2 Modeling of the roof diaphragm............................................................................... 66 4.2 Mathematical Model ........................................................................................................ 66 4.2.1 Static model .............................................................................................................. 66 4.2.2 Dynamic model......................................................................................................... 83 4.3 Model calibration ............................................................................................................. 87 4.3.1 Typical building geometry........................................................................................ 87 4.3.2 Roof diaphragm model calibration ........................................................................... 89 4.4 Model Validation ............................................................................................................. 92 4.4.1 Hand calculation ....................................................................................................... 92 4.4.2 Comparison with model of Devine (2009) ............................................................... 99 4.4.3 Comparison with another numerical model............................................................ 101 5 - NONLINEAR RESPONSE OF TILT-UP BUILDINGS ........ 106 5.1 Seismic performance evaluation of the prototype building ........................................... 107 5.1.1 Seismic hazard ........................................................................................................ 107 5.1.2 Structural response.................................................................................................. 108 5.2 Parametric Study............................................................................................................ 112 5.2.1 Connection configuration - Pushover analysis ....................................................... 112 5.2.2 Connection configuration - Nonlinear dynamic analysis........................................ 119 5.2.3 Coefficient of friction ............................................................................................. 121 vi 5.2.4 Roof diaphragm stiffness and strength ................................................................... 122 6 - SUMMARY AND CONCLUSIONS ......................................... 125 6.1 Summary........................................................................................................................ 125 6.2 Discussion of results ...................................................................................................... 126 6.3 Conclusion ..................................................................................................................... 128 6.4 Recommendation for future research............................................................................. 128 REFERENCES ................................................................................. 130 APPENDIX A. DESIGN NOTES ON TYPICAL BUILDING.... 132 APPENDIX B. TIME-HISTORY ANALYSIS RESULTS .......... 135 APPENDIX C. MODEL CODE...................................................... 141 vii LIST OF TABLES Table 2.1 Connections shear strength (Modified from Lemieux et al., 1998)........................... 22 Table 2.2 Steel deck diaphragm configurations (Essa et al., 2003)........................................... 27 Table 2.3 Shear strength and stiffness in monotonic tests (Essa et al., 2003) ........................... 27 Table 2.4 Shear strength and stiffness in cyclic tests (Essa et al., 2003)................................... 28 Table 3.1 Summary of the experimental tests used to develop the EM5 connection model ..... 34 Table 3.2 Coefficients used for the uncoupled model ............................................................... 48 Table 4.1 Roof model parameters .............................................................................................. 91 Table 4.2 Peak rocking drift reduction .................................................................................... 103 Table 5.1 Summary of the selected ground motions................................................................ 107 viii LIST OF FIGURES Figure 2.1 Formwork layout before casting of concrete tilt-up panels (Devine, 2009)............... 8 Figure 2.2 EM1 joist seat (left - CAC Concrete handbook 2006 , right - Devine, 2009).......... 13 Figure 2.3 EM2 shear plate (left - CAC Concrete handbook 2006; right - Devine, 2009)........ 14 Figure 2.4 Use of EM2 shear plate to connect the diaphragm angle to the panels (CAC Concrete handbook 2006 ) ......................................................................................................... 14 Figure 2.5 EM3 shear plate (CAC Concrete handbook 2006)................................................... 15 Figure 2.6 EM4 shear plate (CAC Concrete handbook 2006)................................................... 15 Figure 2.7 EM5 connection (left - CAC Concrete handbook 2006 , right - Devine, 2009) ...... 16 Figure 2.8 Slab to Panel Connection (Lemieux, Sexmith and Weiler, 1998)............................ 17 Figure 2.9 Free body diagram of tilt-up panels (Olund, 2009).................................................. 18 Figure 2.10 Shear test setup (Lemieux et al., 1998) .................................................................. 20 Figure 2.11 Load-displacement plots of shear tests (Lemieux et al., 1998) .............................. 21 Figure 2.12 Rendering of testing setup (Devine, 2009)............................................................ 24 Figure 2.13 Steel deck roof diaphragm testing setup (Essa et al., 2003)................................... 26 Figure 2.14 IDA drift results from 8 ground motion records (Olund 2009):............................. 30 Figure 3.1 Shear behaviour of the EM5 connector .................................................................... 35 Figure 3.2 Axial behaviour of the EM5 connector .................................................................... 36 Figure 3.3 Trilinear model for cyclic behaviour ........................................................................ 37 Figure 3.4 Unsymmetric connection behaviour......................................................................... 38 Figure 3.5 Force-deformation response regions of the EM5 connector model ......................... 40 Figure 3.6 Optimum coefficient rDisp for each cycle ............................................................... 43 ix Figure 3.7 Model for the increase of the coefficient rForce when maximum displacement is increased .................................................................................................................................... 46 Figure 3.8 Optimum backbone parameters F1, F2 VS Maximum vertical displacement for the coupled tests............................................................................................................................... 49 Figure 3.9 Optimum degradation rate parameter Wref Vs Maximum vertical displacement ... 50 Figure 3.10 Horizontally dissipated energy at failure (left) and shear displacement at failure (right) for all the coupled tests as a function of the maximum vertical displacement ............... 52 Figure 3.11 Backbones of the paired connector model.............................................................. 55 Figure 3.12 Comparison of the shear model with monotonic shear test results ....................... 56 Figure 3.13 Comparison of uncoupled shear model with experimental curves......................... 57 Figure 3.14 Comparison of axial model with monotonic test results ........................................ 58 Figure 3.15 Comparison of model with pure axial cyclic tests.................................................. 58 Figure 3.16 Comparison of axial model with seven uplift tests ................................................ 59 Figure 3.17 Comparison of model with coupled tests where the connection is held-up ........... 60 Figure 3.18 Comparison of model with coupled tests where the connection is brought down to zero axial displacement.............................................................................................................. 61 Figure 3.19 Comparison of the proposed EM5 model with the model of Devine (2009) ......... 62 Figure 4.1 Simplified roof model............................................................................................... 66 Figure 4.2 Free body diagram of a two panel system................................................................ 67 Figure 4.3 Conversion of external forces at the contact corners................................................ 68 Figure 4.4 Local degrees of freedom ......................................................................................... 70 Figure 4.5 The two rocking cases and their coordinates............................................................ 72 Figure 4.6 Inputs and outputs of the contact algorithm ............................................................. 74 Figure 4.7 Inputs and outputs of the uplift algorithm ................................................................ 76 Figure 4.8 Functioning of the uplift algorithm on a 2 panel example ....................................... 77 x Figure 4.9 Inputs and outputs of the rocking algorithm............................................................. 79 Figure 4.10 Functioning of the rocking algorithm..................................................................... 79 Figure 4.11 Functioning of the contact algorithm ..................................................................... 81 Figure 4.12 Modification of external force for equilibrium check ............................................ 82 Figure 4.13 Degrees of freedom defined at the middle of the panels ........................................ 83 Figure 4.14 Use of contact algorithm in dynamics .................................................................... 86 Figure 4.15 General layout of the typical building used for the analysis (Olund, 2009)........... 88 Figure 4.16 OpenSees 2D roof model........................................................................................ 90 Figure 4.17 Simulation of the test using Opensees.................................................................... 91 Figure 4.18 Displacement of configuration 1 ............................................................................ 93 Figure 4.19 Displacement of configuration 2 ............................................................................ 93 Figure 4.20 Two panel model for hand calculation ................................................................... 96 Figure 4.21 Comparison of pushover behaviour of the model with F.Devine's model (2009) for five configurations ................................................................................................................... 100 Figure 4.22 Perform 3D panel model ...................................................................................... 102 Figure 4.23 Comparison of the proposed model with the Perform model using time history analysis with one ground motion record scaled at four amplitudes......................................... 104 Figure 5.1 Design spectrum and spectra of the selected ground motions with 5% damping .. 108 Figure 5.2 IDA results for the 10 selected ground motion records.......................................... 109 Figure 5.3 Pushover loading .................................................................................................... 112 Figure 5.4 Connection failure sequence................................................................................... 114 Figure 5.5 Pushover curves for configurations using two panel-slab connectors per panel.... 115 Figure 5.6 Pushover curves for configurations using three panel-slab connectors per panel.. 116 Figure 5.7 Pushover curves for configurations using four panel-slab connectors per panel ... 116 xi Figure 5.8 Demand parameters for various connector configurations..................................... 119 Figure 5.9 Sliding and rocking demand for various values of the coefficient of friction µ .... 121 Figure 5.10 Demand parameters for various roof diaphragm characteristics .......................... 123 1 1 – INTRODUCTION 1.1 General description of tilt-up construction method Tilt-up construction is a method widely used in North America to build warehouses, schools, and offices. Inspired by timber construction, where the façade of a building was often assembled on the ground before being tilted in place, the first tilt-up buildings appeared in the 1940's. The process is the following: a concrete slab is cast first and is then used as a formwork to cast concrete panels. Once the panels are set, they are lifted and put at their final place using a crane. Compared to regular concrete construction methods, tilt-up requires a very small amount of formwork. In addition, all wall elements can be cast at once. This makes tilt-up buildings generally faster to build and more economical than cast-in-place buildings. Tilt-up construction was primarily used for large, single story box-type warehouse structures with few doors and small openings. Advances in the method allow nowadays for the construction of multiple story buildings with more varied shapes and large openings. In this way, a broad range of structures can now be conceived using tilt-up construction. Moreover, tilt-up structures can be more functional and architecturally pleasing than in the past. In tilt-up buildings, panels are used to resist both the vertical load from the roof and floor and the lateral loading from wind and earthquakes. When it comes to using tilt-up panels as lateral force resisting systems, panels can be divided into two categories: panels with large openings, which behave like moment frames, and solid rectangular panels. The latter have a very particular behaviour compared to regular concrete shear walls. Indeed, solid panels can rock 2 and slide on their foundation, and their behaviour is greatly affected by the type and number of connectors used to link them to the slab and to each other. The specificity of concrete construction has only been taken into account recently in the Canadian building code. Their design was first addressed in 1994 in the Design of Concrete Structure Standards CAN/CSA A23.3-94. Currently, tilt-up buildings in Canada are designed following the recommendations of Weiler (2006) in the Cement Associations of Canada Concrete Design Handbook - Third Edition. 1.2 Work needed The seismic behaviour of solid tilt-up panels systems is currently not well understood. Displacement of panels is mainly due to two motions: sliding and rocking. The deformation of the panel is relatively small as compared to the sliding and rocking motions. The modeling of the sliding and rocking motions is highly nonlinear. Due to the presence of frictional and contact forces, a minimum force is required to trigger panel rocking or sliding. These motions are also responsible for a significant amount of energy dissipation during cyclic loading. The current connection design procedure aims at preventing both sliding and overturning failure. The capacity calculation for those two failure modes are made independently. The necessary number of panel-slab connections is calculated first to resist the base shear of the roof and the panels, without consideration of degradation due to axial connection damage during rocking. Overturning capacity is then calculated by using the panel weight and the strength of the panel-panel connectors, without taking into account the overturning capacity of the panel-slab connections. This procedure is made using static analysis. The ductility of the rocking motion is much higher than the one of the sliding motion. Indeed, on the one hand, the sliding capacity is limited by the panel-slab connection whose failure in shear occurs between 2 mm and 16 mm depending on the connector embedded in the panel. On the other hand, the overturning capacity is limited by the vertical displacement capacity of the panel-slab connection and by 3 the shear capacity of the panel-panel connectors. The vertical displacement capacity of the panel-slab connection is about 100 mm, and the shear capacity of the panel-panel connector is above 33 mm. The ability for the system to deform inelastically is usually taken into account by using a higher force reduction factor (Rd) for the overturning demand (typically 1.5) than for the sliding demand (for which Rd=1.0 is usually used). However, this choice has not been justified by detailed dynamic analysis. The main obstacle facing dynamic analysis of tilt-up buildings using solid panels is the model of the connections. The nonlinear response of the connections plays a key role in the seismic behaviour of the building, and until very recently, limited experimental data was available to understand their mechanical behaviour. Experimental tests have been conducted at the University of British Columbia by Lemieux et al. (1998) to investigate the ductility and the strength of the most common tilt-up connectors. Although these tests provided very valuable information about the cyclic response of the connectors, experimental results were still missing in order to understand the cyclic behaviour of connectors used in tilt-up buildings with solid panels. In particular, no test was performed in the axial direction of the EM5 connector, this direction being the one resisting overturning in the panel-slab connection. A new test series on the connectors has been performed by Devine (2009) at the University of British Columbia. This series includes uplift tests on the panel-slab connection, as well as coupled tests, where the shear behaviour of panel-slab connectors previously damaged by uplift was studied. It showed that uplift of the panel-slab connection can reduce the shear strength significantly. In the work of Devine (2009), the experimental data was used to get the pushover response of many connection configurations. However, dynamic effects as well as cyclic behaviour were not investigated. Olund (2009) studied some potential dissipating energy mechanisms using incremental dynamic analysis. Because of the lack of knowledge about the connector’s behaviour, this study focused on buildings where the panel-slab connectors were ignored. As a result, in order to better understand the seismic behaviour of tilt-up building using rigid panels, a nonlinear dynamic tilt-up model incorporating a realistic connection model is needed. 4 1.3 Research goals and scope The research objectives of this study are to: • Use the available experimental data to develop a numerical model of the EM5 connector. This model has to take the coupling between vertical and horizontal motions into account. • Incorporate this connection model in a building model to study the seismic response of tilt-up structures with rigid panels. • Use the developed finite element models to: - Assess the current design procedure. - Investigate alternative connection configurations. - Conduct parameters study to confirm the sensitivity of important building parameters (such as roof stiffness, coefficient of friction between the panels and the foundation, coupling of vertical and horizontal connection behaviours) on the structural design. 1.4 Thesis overview The remaining chapters of this thesis will address the objectives listed above. Chapter 2 serves as literature review of the current design procedure and the state-of-the-art analytical and experimental work conducted today. Chapter 3 presents a numerical model for the EM5 connector. The hysteresis of the numerical model was calibrated using the experimental data conducted at the University of British Columbia. Coupling between shear and axial behaviour is included. Chapter 4 describes a numerical tilt-up building model that incorporates the connector model presented in Chapter 3. A typical building designed in Olund (2009) is described and implemented. The static and dynamic behaviour of the model were verified. Chapter 5 summarizes the parameter studies of the prototype tilt-up model. Both nonlinear 5 static and dynamic analyses were conducted. Chapter 6 presents the conclusions and recommendations resulting from this study. 6 2 - LITERATURE REVIEW This chapter first summarizes the construction techniques and design requirements of tilt-up construction in Canada. Then, the main research studies performed on tilt-up construction using rigid panels and steel roof diaphragms are presented. Previous research focused primarily on connector response via experimental testing and on building response. 2.1 Review of the current construction techniques and design approach Design approach and construction techniques for concrete tilt-up buildings were summarized in Olund (2009) and Devine (2009). Design guidelines for wall panels and connectors are provided by Weiler (2006) in the Concrete Design Handbook. 2.1.1 Roof diaphragm Design approach and construction techniques for the roof diaphragm are summarized in Olund (2009) and Devine (2009). Typically, two kinds of diaphragms are used in tilt-up buildings in North America. In the United States, the most common is panelized plywood roof diaphragms supported by glulam beams. In Canada, only corrugated cold-formed steel decks are used. The steel decking is supported by steel beams and/or open web steel joists. A steel angle is added along the perimeter of the deck to insures it is supported along the edges. 7 The design of steel diaphragms in Canada follows a document called “Design of Steel Deck Diaphragms – 3rd Edition” (Canadian Sheet Steel Building Institute, 2006). Two methods are described in this document: the Tri-Services method developed by S.B. Barnes and Associates (1973), and the Steel Deck Institute (SDI) method (Luttrell, 2004). The Tri-Services method used to be popular in Canada. However, it can be only applied to steel decks having button punched and seam welded side-lap connections. Recent testing by Tremblay et al., (2002) at the Ecole Polytechnique in Montreal showed that button punched and welded deck connections have extremely limited ductility compared to nail and screw fastened decks. This resulted in a change in design habits, nails and screws being much more common now. The SDI method provides design guidelines as well as equations to compute the stiffness, the strength, and the period of the diaphragm. For seismic and wind design, the main requirement is to provide sufficient capacity in order to avoid fastener failure and shear buckling of the deck. The capacity is adjusted by choosing the adequate deck thickness and fastener type and spacing. Lateral force demand on the diaphragm is computed using a beam model: the maximum shear occurs at the ends of the diaphragm, and the maximum moment occurs at the middle of the roof diaphragm (resulting in compressive and tension forces). Since the maximum shear and tension always occur close to the perimeter of the roof diaphragm, the perimeter of the diaphragm is usually constructed using thicker metal and the interior roof diaphragm is usually constructed using thinner metal. For seismic design, the design forces are typically reduced by a reduction factor which depends on the fastener type used. If screws are used for side laps, and the deck is attached to beams and joists with pin connections, a reduction factor of 1.95 (RdRo = 1.5 x 1.3) can be used. 2.1.2 Wall panels Tilt-up panels are aimed at providing both vertical and lateral load resistance in a building. Design of tilt-up panels for vertical and out-of-plane loading is governed by flexure, since tilt- up panels are often very slender in the out-of-plane direction. Regarding the in-plane direction, 8 the behaviour is usually governed by sliding and rocking mechanisms, which are resisted by friction, panel weight, and the connectors. This section summarizes the construction practices and the current design practices for out-of-plane and in-plane loading of tilt-up panels. i) Construction technique Panel wall construction is described in Devine (2009). During tilt-up building construction, a slab on grade is cast first. Once this slab is set, a chemical that prevents concrete bonding is spread on the top: it will prevent the future panels to bond with the slab. Then formwork and rebars are laid as shown on Figure 2.1. Figure 2.1 Formwork layout before casting of concrete tilt-up panels (Devine, 2009) Concrete tilt-up panels are then cast. Usually, the interior face of the panel is facing up, in order to allow a more convenient placing of the connectors. Once panels are set, they are lifted and put at their final place using a crane. Temporary bracings are needed to hold them in place before the roof diaphragm is installed. 9 ii) Design tilt-up panels for vertical and out-of-plane loading Since tilt-up panels are usually very slender in the out-of-plane direction, the design procedure typically pays special attention to the out-of-plane loading due to both wind and earthquake forces, as well as the eccentricity of the vertical load and the related P-∆ effect. Design guidelines are presented by Weiler (2006) in the Concrete Design Handbook. Eq. 2.1 shows the equation used to check the flexural capacity of tilt-up panels. br MM ≥ (2.1) Where: • Mr is the factored out-of-plane moment capacity; • bbf MM δ= is the factored moment including P-∆ effect; • btf b KP /1 1 − =δ is the moment magnification factor to account for the P-∆ effect; • 0 2 )( 28 ∆+++= tfwftf f b PP eP lw M is the factored moment not including P-∆ effect; • Pwf is the factored panel weight above mid height; • Ptf is the factored axial load at the top of the panel; • e is the axial load eccentricity at the top of the panel; • ∆0 is the initial deflection at panel mid height; • Kb is the bending stiffness of the wall panel; • l is the clear vertical span of the panel; • wf is the factored lateral load. 2010 NBCC requires that a building should be designed to resist seven combinations of different loading types, including dead, live, snow, wind and earthquake loads. For out-of- plane flexure, the following two combinations have to be taken into account: 10 1.25 Dead Load + 1.4 Wind load + (0.5 Live Load or 0.5 Snow Load) (2.2) 1.0 Dead Load + 1.0 Earthquake Load + (0.5 Live Load or 0.5 Snow Load) (2.3) Wind forces are calculated in accordance with 2010 NBCC Clause 4.1.7. In seismic areas, Eq. 2.3 governs. The earthquake forces on panels are calculated using Eq. 2.4 from 2010 NBCC Clause 4.1.8.17 for Elements of structure: ppEaap WSISFV )2.0(3.0= (2.4) Where: • Sa(0.2) is the 5% damped spectral response acceleration at a period of 0.2 seconds. Value based on published climatic data within 2010 NBCC. • Fa is the acceleration-based site coefficient that is adjusted based on the soil conditions and Sa(0.2). • p xrp p R AAC S = where 0.47.0 ≤≤ pS • n x x h hA 21+= is a height factor, where hx, hn are the height above the base level. For out-of-plane forces, hx is taken as the center of mass of the panel. • EI is the importance factor. For normal importance structures, EI equals 1.0. • pW is the panel weight Category 1 of 2010 NBCC Table 4.1.8.17 applies for the design of tilt-up panels for out-of- plane bending: • Cp is the component risk factor. Usually taken as 1.0. • Ar is the dynamic amplification factor. For short period buildings with flexible walls, this is equal to 1.0. If the natural frequency of the component is close to the fundamental period of the building, this factor could increase to as much as 2.5 • Rp is the response factor associated with the ductility of the component. Its value is normally 2.5 for design of reinforced tilt-up wall panels out-of-plane. 11 The slenderness of panels is also limited to 50 and 65 for panels having a single and double reinforcing mat, respectively. iii) Design of tilt-up panels for in-plane loading In the in-plane direction, panels act as the lateral-force resisting system for the building. Therefore, they take lateral loads from roof, floors and out-of-plane walls. As for in-plane loading, 2010 NBCC Clause 4.1.7 is used to calculate wind load, and clause 4.1.8 is used for earthquake load. According to this last clause, the design base shear of the building should be computed using Eq. 2.5: od evava RR WIMTSF V )( , = (2.5) In this equation: • V is the design base shear in kN. • Fa,v is a site coefficient. It depends on the soil condition and the building period. • Sa(T) is the spectral acceleration in unit of g, given by the uniform hazard spectrum at the building’s fundamental period, T. • Mv is a coefficient that takes the shear force due to higher modes into account. In tilt-up buildings, higher modes are usually not considered, and this coefficient is chosen as 1. • Ie is the importance factor. Since tilt-up construction can be used to build either low occupancy warehouses or schools, this coefficients can vary between 0.8 and 1.3. • W is the weight of the building in kN. • Rd is the force reduction factor based on ductility. • Ro is the force reduction factor based on overstrength. This equation is valid as long as the building is either medium size and regular, or small size, or located in an area of low seismicity. Tilt-up buildings usually fall in the "Regular structure" category and medium sized, so this equation can be used most of the time. Otherwise, dynamic 12 analysis would be required. If Rd is higher than 1.5, the base shear does not have to be higher than 2/3* Sa(0.2). Eq. 2.6 is obtained from 2010 NBCC and can be used to assess the fundamental period of the building. 4 305.0 na hT = (2.6) Where hn is the height of the shear wall in meters. This equation does not take into account the flexibility of the diaphragm and is thus giving a period much lower than the one given by modal analysis. This is a conservative estimate as a smaller period increases the design shear forces on the tilt-up building. 2010 NBCC allows engineers to use a higher period (up to a period of 1.5 times the period calculated using Eq. 2.6) calculated by model analysis. The weight of the building, W, includes the weight of the in-plane panels, the roof, half of the out-of-plane panels, and 25% of the design snow load. The bottom of the out-of-plane panels is restrained from sliding in the out-of-plane direction. This restraint is assumed to take half of the forces applied on those panels, that is why only half of out-of-plane panel weight is considered. Panels without connectors have limited lateral capacity. This capacity is provided by contact forces between the panels and the foundation slab. Connection technique and design are described in the next section. 2.1.3 Connectors Contrary to cast-in-place concrete buildings, tilt-up construction requires connectors to link all the structural elements to each other. According to Lemieux et al. (1998), these connectors are used for the transfer of: • Gravity loads from beams and joists to the tilt-up panels. • Shear forces from roof joists and steel deck to the tilt-up panels. 13 • Shear forces from intermediate floors to the tilt-up panel. • Shear forces between tilt-up panels. • Shear forces from tilt-up panels to the floor slab or foundation. • Forces from mechanical equipment or architectural items to the tilt-up panels. Weiler (2006) provides guidelines for the use of five standard connectors in the Concrete Design Handbook. These five connectors are now the most common in Canada. These connectors and their use are described below. EM1 - Joist Seat The EM1 is a joist seat connector made of a steel L89 x 89 x 6 x 300 mm angle which can be embedded in concrete with two 15M Gr 400 rebar anchors. This connection is used as a seat for the openweb steel joists. Figure 2.2 shows a drawing of this connector and a picture before concrete casting. Figure 2.2 EM1 joist seat (left - CAC Concrete handbook 2006 , right - Devine, 2009) 14 EM2 - Shear Plate The EM2 connector is a shear plate connector made of a steel PL150 x 9.5 x 200 plate with two 16mm studs welded. Figure 2.3 shows a drawing and a picture of the EM2 connection. Figure 2.3 EM2 shear plate (left - CAC Concrete handbook 2006; right - Devine, 2009) This connection is usually welded with an edge angle to be connected to the steel deck. Figure 2.4 shows the details of the EM2 connection welded with the steel angle. Figure 2.4 Use of EM2 shear plate to connect the diaphragm angle to the panels (CAC Concrete handbook 2006 ) 15 EM3 - Shear Plate The EM3 is a shear plate connector which is a stronger than the EM2 connection. It is made of a PL 200 x 9.5 x 200 plate with four 16mm studs. Figure 2.5 shows the details of the EM3 connector. Figure 2.5 EM3 shear plate (CAC Concrete handbook 2006) It can be used for the same functions as the EM2 connection. Typical use includes attaching an angle seat for an Open Web Steel Joists (OWSJ) or providing an anchorage in a panel in order to attach this panel to the slab. EM4 - Shear Plate The EM4 is a shear plate connector which is even stronger than the EM3 connector. It is made of a PL 225 x 9.5 x 460 mm plate with eight 16mm studs. Figure 2.6 shows the details of the EM4 connector. The EM4 connector is typically used to attach girders of the steel decking system. Figure 2.6 EM4 shear plate (CAC Concrete handbook 2006) 16 EM5 - Edge angle The EM5 is an edge angle connection, which is made of an L 38 x 38 x 6 x 200 mm steel angle welded to a bent 20M rebar anchor, as shown in Figure 2.7. Figure 2.7 EM5 connection (left - CAC Concrete handbook 2006 , right - Devine, 2009) This connection is used to connect edges of concrete elements. Its long rebar provides sufficiently deep and wide anchorage to avoid cone failure at the angle. It is most commonly used in the panel-panel connection and panel-slab connection. Figure 2.8 shows how the EM5 connection can be used in a slab along with an EM2 or an EM3 connection in a panel in order to make a panel-slab connection. The EM5 and EM3 connectors are welded together. A filler bar is added in between to make the weld stronger. 17 Figure 2.8 Slab to Panel Connection (Lemieux, Sexmith and Weiler, 1998) Connection layout design The choice of the number and the configuration of connectors is based on static analysis. For seismic design, roof shear and panel shear are computed using Eq. (2.5) where half of the mass of the out-of-plane panels is included in the roof shear. The value of RdRo is chosen as 1.0. Using a free body diagram shown in Figure 2.9, the forces in the connectors can be calculated. 18 Figure 2.9 Free body diagram of tilt-up panels (Olund, 2009) The base shear capacity is calculated by adding the frictional force at the base of the panels and the shear capacity of each of the panel-slab connectors. The number of connectors is chosen based on the horizontal capacity needed to resist the base shear generated by the earthquake force. It should be noted that during an earthquake shaking the shear capacity of the bottom connectors might be damaged due to vertical motion during panel rocking. Such behaviour is observed in tests from Devine (2009). However, this phenomenon has not been taken into account in the current design procedure. The overturning demand of the panel is calculated using Eq. 2.7. od wallsplabeinroof of RR hVhV M 2 − + = (2.7) Where: • h is the panel height. • Vroof is the seismic force due to the roof, the out-of-plane panels, and the design snow load. Since there are usually one row of panel on each side of the building, only half of this force is applied in the free body diagram. 19 • Vin-plane walls is the seismic force due to the in plane panels of the row. • Rd and Ro are chosen as 1.5 and 1.3 for the EM5-EM5 connection. The overturning capacity is computed using Eq. 2.8. ( ) bVNbVbNWWM downHoldippppipproofpanelr −+ −++= 1 222 / (2.8) Where: • Nipp is the number of in-plane panels in the building, half of them being in the considered row. • Vp/p is the sum of the shear capacity of all the panel-panel connections between two panels. • VHold-down is the shear force applied by the out-of-plane panels to the in-plane panels through the edge connectors. • b is the width of panels. This overturning capacity is adjusted by adding sufficient panel-panel connectors. It should be noted that the overturning capacity of the panel-slab connectors is not taken into account. 2.2 Review of previous experimental tests conducted on connectors Connectors play an important role in the seismic response of solid wall tilt-up panel buildings. This section summarizes the state-of-the-art experimental tests that have been conducted on tilt-up connectors. 20 2.2.1 Experimental testing by Lemieux, Sexsmith and Weiler (1998) Lemieux, Sexsmith and Weiler (1998) investigated the strength and ductility of the common tilt-up connectors described in Section 2.1.3. Those connectors were subjected to tension and shear tests, under both monotonic and cyclic loading history. The results were then used to recommend the shear and tension strength of those connectors in the Concrete Design Handbook. The tilt-up model developed in this thesis requires only the uplift behaviour of the EM5 connector, and the shear behaviour of the EM2~5 connectors. The former was not investigated in this experimental program, so only the latter is presented here. Figure 2.10. shows the experimental setup for the shear strength tests. Connections were embedded in a 140 mm thick concrete panel with specified strength of 30 MPa. This value was chosen because it is the minimum thickness recommended by the Concrete Design Handbook. Those specimens were cast in a construction site, and then bolted to the strong floor. Figure 2.10 Shear test setup (Lemieux et al., 1998) For each connection type, one monotonic shear test and two cyclic tests were conducted. For the cyclic tests, three cycles were performed at an increasing displacement target. Depending on the tests, between six and eight displacement increments were applied. Figure 2.11 shows 21 the load-displacement plots of the all the shear tests. Test numbers are summarized in Table 2.1. Figure 2.11 Load-displacement plots of shear tests (Lemieux et al., 1998) 22 Based on the results presented in Figure 2.11, the cyclic envelope appears to follow the monotonic plot in the elastic range, but the monotonic envelop appears to overestimate the cyclic envelop in the large deformation range. The EM5 connection appears to be the most ductile connection. It can sustain a monotonic displacement of 15 mm, whereas monotonic tests on EM2, EM3 and EM4 failed between 2 and 5 mm. For cyclic loading, the maximum displacement of the EM2, EM3 and EM4 connectors is between 2 and 4 mm, whereas the EM5 fails at 7 mm. Strength results (in units of kN and kips in brackets) and failure modes of the connectors are summarized in Table 2.1. Table 2.1 Connections shear strength (Modified from Lemieux et al., 1998) The failure of the embedded connectors often implies cone failure of the concrete. So the thickness of a panel is an important parameter. 23 The average strength of the EM5 connector is higher than the one of the EM2 and EM3 connectors. This means that when those connectors are used in series (e.g. in a panel-slab connection, as shown in Figure 2.8), the lower force capacity and ductility from the EM2 and EM3 connector must be used in design. Tension tests were performed on the EM5 connector, but not in the axial direction. 2.2.2 Experimental testing by Devine (2009) Further experimental tests on the EM5 connection were carried out by Devine (2009) at the University of British Columbia. There were two main objectives for this new test series. The first was to determine the uplift behaviour of the EM5 connector, which is crucial to model the panel rocking behaviour. The second was to study the coupling between the axial damage and the shear response. Figure 2.12 shows the experimental setup used by Devine (2009). Such setup allowed both uplift and shear loading. The testing setup was made in a way that reproduces some of the constraints that connections would have in a tilt-up building: it allowed controlled out-of-plane movement, and rotation was restrained. The EM5 connectors were embedded in a piece of concrete whose specified strength was 25MPa. Testing was conducted more than 100 days after casting to reproduce field conditions. 24 Figure 2.12 Rendering of testing setup (Devine, 2009) Four types of tests were performed: • Pure shear tests with both monotonic and cyclic displacement load history. • Pure axial tests with both monotonic and cyclic displacement load history. • Coupled tests, held up: three axial cycles were applied first, the connection was then held up, and finally displacement cycles were applied in the shear direction with increasing amplitude until failure. • Coupled tests, brought down: same as the previous tests, except that the connection was brought back to zero axial displacement before applying the shear load cycles. This protocol for coupled tests made clearer how vertical displacement affects the shear response of the connector than a test where the two directions are tested simultaneously. Tests showed that the axial cyclic ductility of the EM5 connection can be easily improved by changing the type of filler bar, which is an element used to strengthen the weld that attaches the EM5 connector to another element. When the typical round filler bar was replaced by a thicker square bar, the cyclic displacement capacity was increased from 75 mm to more than 100 mm. 25 Pure shear tests were performed in order to validate results from the previous testing program. A different cyclic loading protocol was applied, using higher displacement increments. The results are consistent with previous testing. In all three pure shear tests, failure occurred due to the fracture of the EM5 rebar. Eight coupled tests were performed. Four levels of vertical displacement (25 mm, 50 mm, 75 mm and 100 mm) were tried. For each level, two tests were conducted, one where the connection is held up, one where the connection is brought back to zero vertical displacement. The shear strength and stiffness were found to be greatly affected by the previous axial damage. In particular, the strength during cyclic loading dropped from an average of 220 kN for the undamaged connection to 140 kN for connections having undergone an uplift displacement of 100 mm. 26 2.3 Previous research on roof diaphragm The roof diaphragm plays an important role in the seismic behaviour of tilt-up buildings. The base shear due to the roof diaphragm is greatly affected by its strength and stiffness. Many experimental studies have focused on the static response of the roof diaphragm in the past. A recent testing series at the Ecole Polytechnique in Montreal by Essa, Tremblay and Rogers (2003) investigated the shear stiffness and strength of steel deck diaphragms for both static and monotonic loadings. The aim of this testing series was to validate the design and modeling recommendations from the Steel Deck Institute. Figure 2.13 shows the experimental setup used in this study. Load was applied parallel to the deck direction. The tested dimensions of the tested specimen were 6.10 m x 3.66 m. Figure 2.13 Steel deck roof diaphragm testing setup (Essa et al., 2003) Nine common deck configurations were tested. This includes two deck thicknesses, two deck profiles, six kinds of deck-to-frame fasteners and three kinds of sided lap fasteners. For each 27 configuration, two specimens were tested: one with a monotonic loading and one with a cyclic loading. Table 2.2 shows the characteristics of the specimen tested. Table 2.2 Steel deck diaphragm configurations (Essa et al., 2003) Strength and stiffness of the monotonic tests and cyclic tests are shown on Table 2.3 and Table 2.4, respectively. Table 2.3 Shear strength and stiffness in monotonic tests (Essa et al., 2003) 28 Table 2.4 Shear strength and stiffness in cyclic tests (Essa et al., 2003) The results show that the diaphragm with welds and no washers as deck-to-frame connections have the lowest ductility. On the other hand, diaphragms using B-deck profiles with screwed side lap fasteners and nails or welds with washer have the highest strength during cyclic loading. 2.4 Review of previous analytical models 2.4.1 Olund (2009) Olund (2009) proposed a concrete tilt-up model using the finite element software PERFORM 3D Version 4.0.3 (Computers and Structures, 2007). An analysis was performed following the main steps of the ATC-63 methodology (ATC, 2008). This methodology is a systematic way to assess the performance of a seismic force-resisting system. It is based on the Incremental Dynamics Analysis method proposed by Vamvatsikos and Cornell (2002). The first part of this study was to design archetypical tilt-up buildings. A study amongst tilt-up designers in North America was performed in order to determine the most common characteristics of the roof diaphragm (material, kind of fastener, deck profile, deck thickness, ratio of the roof mass to the mass of the out-of-plane panels, girder and joist span, dead load), 29 wall panels (dimensions , number, arrangement), and connections (type, arrangement). Results were employed to design two typical buildings: one using solid wall panels and one using panels with large openings. Panels were modeled using concrete and steel nonlinear fibre. Roof diaphragm and its connections to the panels were assumed to be linear elastic. Contact and friction at the panel- slab and at the panel-panel interfaces were modeled using nonlinear materials (a compression only elastic spring material for contact and an elastic perfectly plastic material for friction). The monotonic behaviour of the panel-slab and panel-panel connectors was modeled using the results from the first test conducted by Devine (2009). Coupling between horizontal and vertical motions was not modeled. Since the connection model was monotonic, the seismic behaviour of the building could not be investigated. As a result, most of the work on solid panels was performed without panel-slab connectors, a situation that would occur when all the connectors have failed. Three models were derived to study possible energy dissipating mechanisms: • A panel sliding model, where panel rocking was restrained and the panel-slab connections of the in-plane panels were removed. • A panel rocking model, where panel sliding was restrained and all the connectors of the in-plane panels were removed. • Frame action, for panels with wide openings. A set of 8 ground motion records was applied to the building. Records were scaled such that the spectral acceleration at the fundamental period ranges from zero to 3g. Figure 2.14 shows the median maximum drift of panel top and the drift at the middle of the roof. The displacements are presented for both the sliding model and the rocking model. The sliding drift was found to be smaller than the rocking drift. However, it also leads to a high residual displacement that would be nearly impossible to repair. 30 Figure 2.14 IDA drift results from 8 ground motion records (Olund 2009): The model was also used to propose a value of the RdRo factor for the rocking and the frame mechanism. The 2005 NBCC method assesses the RdRo value based only on the connectors behaviour without taking the building mechanisms into account. This study showed that an RdRo value of 2.1 for the rocking mechanism resulted in a collapse probability lower than 10% for a hazard of 2% in 50 years. It was concluded that a rocking mechanism would be more practical because of the residual displacement of the sliding mechanism. For the frame mechanism, the value of Ro =1.3 and Rd=1.5 that are commonly used in design were found to be unconservative because they resulted in high shear demand in the roof diaphragm. The natural period of the sliding model was 0.58s, which is within 10% of the period predicted by ASCE41-06 when half of the mass of the out-of-plane panels is added to the mass of the roof. This period is much higher than the one predicted by 2010 NBCC in Eq (2.6), which is 0.26s. 31 2.4.2 Devine (2009) In a study by Devine (2009), strength, ductility, and failure mode of tilt-up panel systems using various connection configurations were studied. The results of the experimental tests described in 2.2.2 were used to develop a backbone model of tilt-up connectors. This model was then incorporated in a panel row model. A monotonic pushover study was performed using this model. The model was two dimensional with the assumption that the panels are rigid, and the effect of the steel angle and of the connections to out-of-plane panels were neglected. Studied configurations include the combination of 1 to 4 panels, 2 to 4 panel-slab connectors, and 2 to 4 panel-panel connectors. The strength of the different systems was compared. It was found that adding panel-panel connectors increases significantly the strength of a system. However, this strength was found to be very brittle. A total of four failure mechanisms were predicted using this static monotonic pushover model, each of them having resulting in a specific behaviour: • Overturning without uplift: this failure mode occurs when there are not sufficient panel- panel connectors to uplift any panel before they fail. • Overturning with panel uplift: this failure mode occurs when there are enough panel- panel connectors top uplift one or two panel during the rocking motion. Contrary to the previous failure mechanism, all the panel-panel connectors do not fail simultaneously. As a result, they increasing the system resistance in a larger drift range. • Sliding: this failure mode typically occurs when there is a lot of panel-panel connectors and few between panels and slab. This failure mode occurs in 13 configurations out of 25: this is the most common one. • Extremely brittle failure: this failure mode occurs when there are so many panel-panel connectors that panel-slab connectors fail before the panel-panel ones. This results in failure at a very low panel drift. 32 A study was then led on the performance of ideal connectors. It was suggested that two characteristics would be desirable: increasing the uplift displacement capacity of the connectors, and having a minimum coupling between the shear and the axial behaviour of the connectors. 2.4.3 Need for a new model In the past, Devine (2009) and Olund (2009) have each developed a numerical model to simulate the seismic behaviour of tilt-up buildings. Devine’s model was developed using simple empirical equations to estimate the force- deformation backbone of the EM5 connector. This model was used to simulate the peak force- deformation response of the tilt-up panels with EM5 connections under monotonic loading condition. The model, however, was unable to simulate the force-deformation response of the EM5 connection under cyclic response. Olund (2009) used the Perform3D software to developed a full three dimensional nonlinear model for a tilt-up building, however the EM5 connectors were not implemented. To overcome deficiencies of the past models, a numerical model of the EM5 connector was developed in Chapter 3 using the experimental tests conducted by Devine (2009) and Lemieux et al. (1998). This model was then incorporated in a building model using the software MATLAB, as explained in Chapter 4. 33 3 - NUMERICAL MODEL FOR THE EM5 CONNECTOR The seismic behaviour of concrete tilt-up buildings is highly influenced by the cyclic behaviour of the connectors. Amongst all the connectors used in the tilt-up industry to link concrete elements to each other, the EM5 connector is the only connection recognized by the Concrete Design Handbook to have a significant ductility. In this chapter, a numerical model of the EM5 connector is developed based on the results obtained from the experimental tests from Devine (2009) and Lemieux et al. (1998). 3.1 Summary of experimental tests used in this study All the available tests on the EM5 connectors in the shear and uplift direction were used to calibrate the nonlinear EM5 model. They are summarized in Table 3.1. 34 Table 3.1 Summary of the experimental tests used to develop the EM5 connection model Specimen name Loading protocol Source Axial only tests 2 Monotonic Devine (2009) 3 Cyclic Devine (2009) 8 Cyclic Devine (2009) 11 Monotonic Devine (2009) Shear only tests S1 Monotonic Lemieux et al. (1998) S2 Cyclic Lemieux et al. (1998) S3 Cyclic Lemieux et al. (1998) 5 Monotonic Devine (2009) 6 Cyclic Devine (2009) 12 Monotonic Devine (2009) Coupled tests 15 3 vertical cycles at 25mm, brought back to 0, then cyclic shear Devine (2009) 14 3 vertical cycles at 50mm, brought back to 0, then cyclic shear Devine (2009) 10 3 vertical cycles at 100mm, brought back to 0, then cyclic shear Devine (2009) 16 3 vertical cycles at 25mm, held up, then cyclic shear Devine (2009) 13 3 vertical cycles at 50mm, held up, then cyclic shear Devine (2009) 7 3 vertical cycles at 75mm, held up, then cyclic shear Devine (2009) 9 3 vertical cycles at 100mm, held up, then cyclic shear Devine (2009) Four tests were available to calibrate the model coefficients for the axial behaviour, and six for the uncoupled shear behaviour. No test was available from the study of Lemieux et al. (1998) regarding the axial behaviour of the EM5 connector. Seven tests were available to calibrate the coupled model. In the experimental study of Devine (2009), two connection details were tested: one using a typical rounded filler bar, and one using a square filler bar. Those two connection details were treated in the same way in the model. In fact, the force-deformation plots were very similar, the only significant difference being the axial displacement cyclic capacity. 35 3.2 General mechanical behaviour of the EM5 connector In tilt-up buildings, EM5 connections are mostly utilized in the shear and axial directions. The responses in these two directions are quite different and hence the response in these two directions will be calibrated separately. It should be noted that, based on the experimental data from Devine (2009), the shear response is affected by the axial response. However, the axial behaviour is assumed to be unaffected by the shear response. 3.2.1 Monotonic and cyclic envelope Figure 3.1 and 3.2 show the force-deformation response of the EM5 connection under shear and axial load, respectively. Figure 3.1 Shear behaviour of the EM5 connector 36 Figure 3.2 Axial behaviour of the EM5 connector As shown in these figures: • The force-deformation envelope of cyclic tests is very close to the monotonic tests. For shear tests under cyclic load, the force envelope tends to cap beyond the yielding displacement, while the monotonic force-deformation plot tends to increase gradually with increasing deformations. On the other hand, for axial tests, the cyclic envelope is almost identical to the monotonic envelope. A difference can only be observed between 40 and 90 mm. This difference is due to spalling of the concrete cover of the connector. This phenomenon occurs at very random displacements depending on the axial test, so it was ignored in the modelling. • For both monotonic and cyclic cases, the force-deformation envelope can be approximated well using a bilinear model. 37 3.2.2 Cyclic behaviour and damage The force-deformation response of the EM5 connection under cyclic loading can be approximated using a tri-linear line as shown in Figure 3.3. The three parts of this trilinear model can be described as unloading, transition, and reloading. Axial and shear behaviour both follow similar patterns, hence they can be modeled in a very similar way. The main difference comes from the fact that negative displacement in the axial behaviour is impossible, hence restricted. Figure 3.3 Trilinear model for cyclic behaviour The strength and the stiffness of the tri-linear lines tend to degrade significantly after repeated cycles. When multiple cycles are performed at the same displacement, the reloading part starts at a higher displacement and a lower force. This shows the EM5 connection is not capable of retaining the strength after it suffers the inelastic damage. Nevertheless, when the connection is pushed to a new maximum displacement, the connection seems to be able to resist higher force similar to the case of the monotonic loading condition. This shows that modelling of the damage is a key part of the numerical model. 38 3.2.3 Coupling Based on the experimental tests conducted by Devine (2009), when the EM5 connection is damaged under axial loading, then the stiffness, strength, and degradation rate in the shear direction are affected. For example, the cyclic shear strength, which is about 230kN in uncoupled tests, drops to 140kN in tests that have experienced an axial uplift of 100 mm beforehand. 3.2.4 Symmetry Most of the tests experience nearly symmetrical response, where the maximum negative strength is within 10% of the maximum positive strength. However, in some special cases such tests 15 and 16 from Devine (2009) are shown on Figure 3.4, the behaviour can be unsymmetrical. (a) , (b) Test 15 and Test 16 from Devine (2009) Figure 3.4 Unsymmetric connection behaviour 39 In Test 15 (Figure 3.4 (a)) the maximum positive force is 25% stronger than the maximum negative force (+239kN vs. -191kN). In Test 16 (Figure 3.4 (b)), the maximum negative force is 50% higher than the maximum positive force (-249kN vs. +166kN). However, since the design and loading history for specimen 15 and 16 is relatively symmetrical and without additional data to back up the difference, it is assumed that the force-deformation response of the EM5 connector will always be symmetrical. 3.3 Proposed analytical model A numerical model was developed to model the cyclic response of the EM5 connector. Both the coupled and uncoupled response were modeled. 3.3.1 Mathematical rules for the uncoupled model i) General description The model takes a displacement as the input, and returns the corresponding resisting force. The resisting force will account for the strength and stiffness degradation from the preceding cycles. To get the restoring force for any displacement, the force-displacement domain is separated in three regions, namely Region 1, Region 2 and Region 3. The force-deformation response in each region is calculated based on the stiffness rule defined in each region. Figure 3.5 shows the detail as where regions 1, 2 and 3 are defined. It can be noted that the regions presented in Figure 3.5 represent the case where the displacement is increasing (positive increase). Symmetrical rules apply if the displacement decreases (i.e. goes towards the negative direction). 40 Figure 3.5 Force-deformation response regions of the EM5 connector model The regions are defined as follows: • Region 1 includes the following three sub areas: 1) all the area below the force level FB; 2) the area where displacement is higher than DA and force is lower than FA; 3) The triangular area below the point A, which is defined using the elastic stiffness. • Region 2 is the area where the displacement is lower than DA and the force is higher than FB, except the triangular area described in Region 1. Because of the model response rules, the force in this region is always lower than FA. • Region 3 is the area where the displacement is higher than DA and the force is higher than FA. In general, if the force-deformation response is in region 1, the response follows the elastic stiffness kel. If the force-deformation response is in region 2, the response goes toward point A. If the force-deformation response is in region 3, the response follows the reloading stiffness kr. 41 The force should neither exceed the backbone curve, and it should also not be outside the coloured areas of Figure 3.5. In order to deal with decreasing displacements, a second point A is defined on the negative side. The name of A in the positive and negative direction is defined as A+ and A-, respectively. The position of the points A+ and A- depends on the displacement history and the dissipated energy. Eq. 3.1 to 3.4 show the rules to calculate the position of point A+ and A-: DA+ = rDisp * Dmax+ (3. 1) DA- = rDisp * Dmax- (3. 2) FA+ = rForce+ * FBB+ (3. 3) FA- = rForce- * FBB- (3. 4) Where: • rDisp, rForce+ and rForce- are dimensionless coefficients which are calculated using Eq. (3.7), (3.8) and (3.9), respectively. The values of these two parameters depend on the dissipated energy and the past displacement history. • Dmax+ and Dmax- are the maximum positive and negative displacement experienced by the connector, respectively. • FBB+ and FBB- are the values of the backbone at Dmax+ and Dmax-, respectively. The same rules apply to the modeling of the connector in the axial direction, except DA- is set to zero and FA- is set to at a constant value FAn. Using the rules described above, the force-displacement point follows the backbone during the first cycle. Indeed, displacement and forces are always higher than the ones of the point A during the first cycle. As a result, the stiffness is given by the reloading slope, which is the same as the elastic slope in the first cycle. When the yielding displacement is reached, the force becomes limited by the backbone force. 42 As shown on Figure 3.5, the backbone is defined by its yielding point (force F1, displacement D2) and a second point on the yielding slope (Force F2, displacement D2 , D2 being arbitrary). ii ) Cyclic damage parameters Two parameters have been found to be relevant to model damage: the energy dissipated by the connection and the maximum displacement. The dissipated energy is stored on an internal variable W. At every displacement step, this energy is calculated by the area under the force-deformation response curve, which can be simplified using a trapezoidal rule as shown in Eq. 3.5. W(t+1) = W(t)+0.5*(F(t+1)+F(t))*(D(t+1)-D(t)) (3. 5) where F(t), D(t), and W(t) are the force, displacement, and dissipated energy at time step t, respectively. At t=0, W(0) is set to zero. The value of the dissipated energy affects the reloading slope (in region 3) as well as the position of the points A+ and A-.The value of Dmax+ and Dmax- affects only the position of the points A+ and A-, respectively. A constant reference energy, Wref, is defined to normalize W every time it appears in a formula. This parameter enables the change of the rate of all the phenomena due to degradation at the same time. It will be used to model the coupled behaviour of the connectors, as explained in section 3.3.3. iii) Reloading slope When multiple cycles are performed, the reloading slope in Region 3 decreases. This phenomenon is computed using Eq. 3.6 to get the reloading slope kr: Wref* W1 21 2 α + − += krkrkrkr (3. 6) In this equation, kr1 and kr2 are constant coefficients which represent the reloading slope in the initial and final cycles, respectively. In short, as intended in the equation, at the earlier phase of 43 the loading, the value of W is small, kr reduces to kr1. As W increases, kr reduces to kr2. The values of kr1 and kr2 can be easily measured from the force-deformation curve of the experimental data. The coefficient α represents the degradation rate of the reloading slope. iv) Unloading slope The unloading slope is modeled using the elastic stiffness kel, which is constant. v) Position of the point A * Coefficient rDisp The coefficient rDisp represents the ratio of the displacement of the point A to the maximum displacement. To simplify the model, the same coefficient rDisp is used for both the positive and negative direction. Figure 3.6 (a) and 3.6 (b) show how the coefficients rDisp was fitted as a function of W/Wref for the axial and shear tests, respectively. These plots are obtained from Test 5 and Test 3 from Devine (2009), respectively. (a) Axial test #3 ; (b) Shear test #5 Figure 3.6 Optimum coefficient rDisp for each cycle As shown in this figure, the value of rDisp increases almost linearly as a function of W/Wref. Linear regression was used to compute the relationship of rDisp as a function of W/Wref , as shown on Eq 3.7: 44 rDisp = W/Wref*rDisp1+rDisp2 (3. 7) where rDisp1 and rDisp2 are the coefficients obtained from the linear regression. Their value is shown in Table 3.2. * Coefficient rForce The force level of the point A was found to depend on more than one parameter, and a linear regression model (as for the coefficient rDisp) gives a poor approximation. The force position of the point A gets lower when multiple cycles are performed at the same maximum displacement, but it gets higher when a new maximum displacement is achieved. As a result, the parameter rForce needs to account for both the dissipated energy W and the variation of the maximum displacement. Equations 3.8 and 3.9 shows the equation developed to calculate rForce+ and rForce-, respectively: rForce+ = rForceRef*(1-rForceM+) (3. 8) rForce- = rForceRef*(1-rForceM-) (3. 9) where rForceRef is a constant coefficient which represents the ratio of the force A with respect to the value on the backbone curve. The value of rForceM+ and rForceM- is calculated using Eq 3.10 and 3.11. The value of rForceM - and rForceM+ depends on the maximum displacement, Dmax, and the energy dissipated by the connector, W. + + ++ +++ ∂ ∂ + ∂ ∂ +=+ max max max )()()(),,1( D D rForceMW W rForceM trForceMDWtrForceM δδ (3. 10) − − −− −−− ∂ ∂ + ∂ ∂ +=+ max max max )()()(),,1( D D rForceMW W rForceM trForceMDWtrForceM δδ (3. 11) Where: • rForceM+(t) and rForceM-(t) are the values of rForceM+ and rForceM- at time step t, respectively. The initial value of those two coefficients is zero. 45 • W rForceM ∂ ∂ + )( is the derivative of rForceM+ with respect to W. • Wδ is the increment of W between time t+1 and time t. • + + ∂ ∂ max )( D rForceM is the derivative of rForceM+ with respect to Dmax+. • +maxDδ is the increment of Dmax+ between time t+1 and time t. The force level of the point A was found to decrease quickly in the first cycle, and more slowly in the later cycles. As a result Eq. 3.13 is used to describe the variation of rForceM+ and rForceM - with W: rForceM+ = 1 - exp(-W/βWref)+f(Dmax+) (3. 12) rForceM - = 1 - exp(-W/βWref)+f(Dmax-) (3. 13) Where: • β is a dimensionless coefficient; • f(Dmax+) is a term unrelated to W. Taking derivative of Eq 3.12 and 3.13 with respect to W: Wref WrefW W rForceM β β )/exp()( − = ∂ ∂ + (3. 14) Wref WrefW W rForceM β β )/exp()( − = ∂ ∂ + (3. 15) The variation of rForceM+ and rForceM- with an increase of the maximum displacement is set using the following rule: if the maximum displacement is increased, the force increases as well; the ratio of the force increase to the force on the backbone (at that displacement) is used to increase rForce+. This is explained on the example sketched in Figure 3.7, where γ is a fixed coefficient. 46 Figure 3.7 Model for the increase of the coefficient rForce when maximum displacement is increased During the increase of the maximum displacement, the force increases by a value dF. The point A, which was at the position Aold before the increase, goes up to the position Anew. For example, if the value of dF is 20% of the backbone force, then A goes up by γ*20%. The mathematical formula that appears on Figure 3.7 is written for rForce+ and rForce- in Eq 3.16 and 3.17: ++ + + + = ∂ ∂ BBF dF rForce D D rForceM γ δ max max )( (3. 16) −− − − − = ∂ ∂ BBF dF rForce D D rForceM γ δ max max )( (3. 17) Taking the derivative of Eq 3.8 and 3.9 with respect to Dmax gives: + + + + + + ∂ ∂ −= ∂ ∂ max max max max )( * )( D D rForceMfrForceRD D rForce δεδ (3. 18) − − − − − − ∂ ∂ −= ∂ ∂ max max max max )( * )( D D rForceMfrForceRD D rForce δεδ (3. 19) 47 Eq 3.20 and 3.21 are derived by substituting Eq 3.18 and 3.19 into Eq 3.16 and 3.17: frForceR rForce F dF frForceR D D rForce D D rForceM BB ε γ ε δ δ + + + + + + + + −= ∂ ∂ −= ∂ ∂ maxmax max max )( )( (3. 20) frForceR rForce F dF frForceR D D rForce D D rForceM BB ε γ ε δ δ − − − − − − − − −= ∂ ∂ −= ∂ ∂ maxmax max max )( )( (3. 21) The variable Dmax is then introduced using the reloading slope: dF=kr*Dmax (3. 22) This leads to these final equations: + + + + + + −= ∂ ∂ maxmax max * **)( DfrForceRF rForcekD D rForceM BB r δ ε γδ (3. 23) − − − − − − −= ∂ ∂ maxmax max * **)( DfrForceRF rForcekD D rForceM BB r δ ε γδ (3. 24) In this process, the value of rForceM+ and rForceM- is forced to remain positive. As a result, the point A is never higher than the backbone force multiplied by rForceRef. One could argue that this model is quite complicated. However, a basic linear model where only the dissipated energy is taken into account was also tried. It resulted in a very poor approximation of the force level of the point A. * Coefficient FAn The force level of the point A- was assessed for all the cyclic axial tests. The average value was used to define the coefficient FAn. 48 3.3.2 Coefficients of the uncoupled model Model coefficients for the uncoupled shear and axial behaviour were calibrated using the tests described in section 3.1. They are summarized in Table 3.2. Table 3.2 Coefficients used for the uncoupled model Shear behaviour Axial behaviour Backbone parameters F1 (kN) 207 34 F2 (kN) 254 80 D1 (m) 0.0013 0.002 D2 (m) 0.015 0.06 Normalization coefficient for dissipated energy Wref (kN.m) 21.1 16.7 Reloading slope parameters kr1 (kN/m) 175 17 kr2 (kN/m) 13.2 3.4 α 0.05 0.025 Displacement of point A rDisp1 0.00064 0.012 rDisp2 0.34 0.77 Force of point A rForceRef 0.5 0.4 β 1 1 γ 3 4 FAn (kN) - -38 Upper force limit of region 1 FB (kN) -14 0 3.3.3 Model of coupling As shown in the experimental results presented in Devine (2009), the shear force-deformation response is affected by the axial deformation. Hence, a coupling model is developed. Note that this coupling model does not describe how axial response is affected by the shear response as no experimental test is available in this aspect. 49 The uncoupled shear behaviour is modified in the following way to take coupling into account: • The backbone curve is lowered to model the strength and stiffness reduction. This is done by modifying the value of the envelope parameters F1 and F2. • The reference energy Wref is lowered to take into account the fact that shear damage occurs more quickly after the connection has undergone some axial damage. • The height of the point B is adjusted to fit the average of all the shear tests instead of fitting only the uncoupled tests. i) Modification of the backbone For each coupled test from Devine (2009), the values of F1 and F2 that fit best the envelope of the cycles were determined. Those optimum parameters F1 and F2 have been plotted versus the maximum vertical displacement and the vertical dissipated energy. Using either of those vertical damage parameters, a linear approximation of the coefficients F1 and F2 was found to work well. In order to keep the model simple, the maximum vertical displacement was used as the axial damage parameter. Figure 3.8 shows the values of the optimum parameters F1 and F2 plotted versus the maximum axial displacement. 0 50 100 150 200 250 0 25 50 75 100 Maximum axial displacement (mm) O pt im u m F1 (kN ) 0 50 100 150 200 250 300 0 20 40 60 80 100 Maximum axial displacement (mm) O pt im u m F2 (kN ) Brought back to zero Held up Linear (Brought back to zero) Linear (Held up) Figure 3.8 Optimum backbone parameters F1, F2 VS Maximum vertical displacement for the coupled tests 50 The difference between tests where the connector is kept at the maximum vertical displacement and tests where the connector is brought back to zero vertical displacement was found to be very small. This is a good validation for the choice of the maximum vertical displacement as the only axial damage parameter. The following equations for F1 and F2 were derived using linear regression: F1 = -1.08 kN/mm * max(Dv) + 205kN F2 = -1.15 kN/mm * max(Dv) + 258kN (3. 25) Where max(Dv) is the maximum axial displacement of the connector. ii) Modification of the shear damage rate For each test from Devine (2009), the degradation rate Wref that gave the best fit was determined. Here are the optimum values of Wref plotted versus the maximum vertical displacement: Figure 3.9 Optimum degradation rate parameter Wref Vs Maximum vertical displacement These points can be well approximated by a linear model. A least square fitting gives the following relation: Wref = -0.12 * max(Dv) +21.1 kN.m (3. 26) 51 iii) Adjustment of the unloading end The ideal position of the point B has been estimated for all the shear tests. In tests by Devine (2009), the height of B does not change in a noticeable way when the damage due to uplift is increased. Therefore, this parameter is chosen independently of damage. The average heights on the negative and positive sides are close to each other and therefore the average height is chosen: FB = 33.1kN This average force is higher than the one derived in the uncoupled model. However, this is unlikely to be due to axial damage. Indeed, FB is nearly the same for all of F.Devine's tests, and the value of FB in those tests is much higher than in K.Lemieux's tests. As a result, this force increase is due to the fact that more of F.Devine's tests are taken into account in this new average. 3.3.4 Model of failure All of the specimens described in Table 3.1 were tested until failure. For most of them, failure happens when the displacement is increased to a value higher than the maximum displacement of the previous cycle (tests 14, 16, 10, 13, 6 and 7). However, failure sometimes occurs in the middle of a cycle (test 15). Such failure can be simulated in the numerical model in two ways: • Based on the amount of dissipated energy; • When a certain displacement is reached. Those modelling approaches were studied for both the shear behaviour and the axial behaviour. i) Coupled shear failure 52 In the coupled tests, the shear failure might be affected by the axial force-displacement history. Potential trends for the two failure parameters described above (dissipated energy in shear and maximum shear displacement) were studied using the following axial response variables: • The maximum vertical displacement, • The current vertical displacement, • The energy dissipated vertically. Using the maximum vertical displacement as the axial response variable gave the best results. No trend could be observed when the current vertical displacement or the energy dissipated vertically was used. Hence, the maximum vertical displacement was chosen as the variable describing the axial response in the coupled shear failure model. Figure 3.10 shows how the dissipated energy at failure and the maximum horizontal displacement at failure vary with the maximum axial displacement for all the shear tests (coupled and uncoupled): 0 5 10 15 20 25 30 0 20 40 60 80 100 Max Dv (mm) W (kN . m ) y = 0.0202x + 15.436 0 5 10 15 20 25 0 20 40 60 80 100 120 Max Dv (mm) D h (m m ) Figure 3.10 Horizontally dissipated energy at failure (left) and shear displacement at failure (right) for all the coupled tests as a function of the maximum vertical displacement Points are very sparse when the dissipated energy in shear is used as the variable describing failure. However, a trend can be observed if the maximum shear displacement is used instead. This parameter seems much more reliable to predict the horizontal failure of the connector. Therefore, shear displacement is chosen as the failure parameter. 53 The maximum displacement does not depend significantly on the vertical displacement since the trend line is nearly horizontal. As a result, shear failure was modeled to occur at a shear displacement that does not depend on the vertical displacement history. This displacement was chosen as 16.5mm, which is the average value of the maximum displacement at failure for all the shear tests considered in this model. ii) Axial failure The study of axial failure is simpler since the axial behaviour is uncoupled. Only four axial tests were conducted until failure (namely tests 1, 2, 3 and 11 from F.Devine (2009)). Similarly to the shear behaviour, using the maximum axial displacement as the failure parameter resulted in a smaller dispersion than when the dissipated energy is used. As a result, the failure displacement was chosen as the average failure displacement of those four tests, which is 98mm. iii) Failure implementation If the connection displacement is higher than the failure displacement at any time step, then the force returned for any displacement is zero. However, the post-failure stiffness is bounded in order to avoid having to deal with a very high negative stiffness in the building numerical model. The value chosen to limit the post-failure stiffness is twice the elastic stiffness. 3.3.5 Model of paired connectors In buildings, connectors are always assembled two by two. Typically, EM3-EM5 connections are used to link panels to the slab, and EM5-EM5 connections are used to link panels to each other. Tests have been performed only on EM3 and EM5 connectors alone, so the behaviour of paired connectors has to be assumed. 54 It is assumed that, in an EM5-EM5 connection, each connector has the same force-deformation response as observed in the experimental tests. Most likely, the behaviour will be a little bit different because of the different restraints between an EM5-EM5 connection and the experimental setup. However, this difference is likely to be small, and no experimental data is available in this aspect. As a result, the EM5-EM5 force-deformation response was obtained using the EM5 model with doubled backbone displacement parameters (D1 and D2). A detailed model of the EM5-EM3 connection would be complicated since a model of the EM3 connection would have to be developed first. However, testing results from Lemieux et al. (1998) show that the EM3 connector has a very small displacement capacity and that its shear strength is lower than the EM5's shear strength - the EM3 shear strength being however higher than the EM5 axial strength. As a result, the EM3 connector will much likely fail before the EM5, thus triggering failure at a much lower displacement than the one predicted by the EM5 model. Contrary to the EM3 connector, the EM4 connector is significantly stronger than the EM5 connector. Hence, it can be assumed that the EM4 connector will always remain elastic when paired in series with an EM5 connector. As result, the EM5-EM4 connection model can be obtained from the EM5 model by increasing the backbone parameters D1 and D2. The EM5-EM4 is very rarely used in practice to connect the panels to the slab. However, this model enables the study of the seismic response of tilt-up buildings where the whole shear ductility capacity of the EM5 connector is used in the panel-slab connector. It can be noted that using EM4 connectors to link the panels to the slab would be expensive. However, in the two following situations, the ductility of the EM5 connection could be better used without having to resort to EM4 connectors: • Using connectors with six studs. According to the tests from Lemieux et al. (1998), a connector having a strength intermediate between the EM3 and the EM4 connectors would be sufficient to keep the EM5 connector as the ductile fuse. 55 • Having thick enough panels. In the tests from Lemieux et al., panels are 140mm thick. One of the tested EM3 connectors was found to fail in a cone failure mechanism. Having a thicker panel would thus increase the EM3 connector strength. Figure 3.11 shows the backbone of the actual EM5-EM3, EM5-EM4, and EM5-EM5 connections. The axial behaviour of the EM5-EM5 model is not plotted since it will not be used in the tilt-up model. Figure 3.11 Backbones of the paired connector model 56 3.4 Validation of the model The analytical model was validated for each of the test in Table 3.1. The displacement history of each test was extracted and smoothed. This displacement was used as an input in the model developed and calibrated in the previous sections. The force given by the model was then compared with the experimental force. 3.4.1 Shear behaviour Figure 3.12 shows how the envelope of the shear model compares to results from monotonic shear tests from Devine (2009) and Lemieux et al. (1998): -50 0 50 100 150 200 250 300 0 10 20 30 40 Displacement (mm) F o rc e ( k N ) Devine - Test 5 Devine - Test 12 Lemieux - Test S1 Numerical model Figure 3.12 Comparison of the shear model with monotonic shear test results Figure 3.13 shows the comparison of the force given by the model with the cyclic shear response of three uncoupled tests. 57 Top left: Test 6 from Devine (2009) ; Top right: Tests S2 from Lemieux et al. (1998) Bottom left: Test S3 from Lemieux et al. (1998) Figure 3.13 Comparison of uncoupled shear model with experimental curves The model response follows closely enough the experimental curve. The comparison to the monotonic test is satisfying, and the strength degradation at each cycle is well approximated. In test S3 (Figure 3.13 (c)), the degradation occurs too fast on the positive side. However, it occurs also a bit too slowly in the other tests. 58 3.4.2 Axial behaviour Results from Test 2 from Devine (2009) are plotted on Figure 3.14 along with the axial force- deformation envelope of the model: Figure 3.14 Comparison of axial model with monotonic test results Figure 3.15 and 3.16 show the comparison of the axial model response with the axial cyclic tests from Devine (2009). Tests on Figure 3.16 were coupled, so the axial loading was not applied until failure. The numerical curves are close enough to the experimental ones. In particular, the strength at each cycle and the initial stiffness are well approximated. left: test 3 - right: test 8 Figure 3.15 Comparison of model with pure axial cyclic tests 59 Top left: test 16 - top right: test 15 2nd line left: test 13 - 2nd line right: test 14 3rd line left: test 7 - 3rd line right: test 9 Bottom left: test 10 Figure 3.16 Comparison of axial model with seven uplift tests 60 3.4.3 Coupled behaviour Figure 3.17 and 3.18 show how the force output of the model compares to the held-up and brought down test results from Devine (2009), respectively. The x-axis shows the displacement history of the connector, in meters. Top left: test 16 (Uplift 25mm) - Top right: test 4 (Uplift 50mm) Bottom left: test 7 (Uplift 75mm) - Bottom right: test 15 (Uplift 100mm) Figure 3.17 Comparison of model with coupled tests where the connection is held-up 61 Top left: test 6 (Uplift 0mm) - Top right: test 15 (Uplift 25mm) Bottom left: test 14 (Uplift 50mm) - Bottom right: test 10 (Uplift 100mm) Figure 3.18 Comparison of model with coupled tests where the connection is brought down to zero axial displacement 62 For each test, the general behaviour is well approximated and the force from the model is satisfyingly close to the experimental force. In some tests, the strength is not approximated well on one side. For some tests, one side is too strong and the other is not strong enough. This is due to the fact that a symmetric model was chosen. 3.4.4 Comparison with model of Devine (2009) A simple model of the connector was defined in Devin (2009). The model is only able to simulate the backbone response of the shear behaviour using a trilinear curve. The axial response was modeled using a bilinear curve. Thus, the model was limited to only monotonic loading. Figure 3.19 shows how the comparison of the force-deformation response between the proposed model and the model presented in Devine (2009). (a) Uncoupled shear model; (b) Axial model Figure 3.19 Comparison of the proposed EM5 model with the model of Devine (2009) 63 The two models give very close responses. The main difference is the shear strength: the proposed model is 45kN stronger than the model of Devine (2009). This is due to the fact that Devine's trilinear model was based on the envelope of cyclic tests, whereas the envelope of the proposed model is fitted using both the envelope of cyclic test and the monotonic tests. 3.4.5 Sensitivity Any change of the coefficient α must be at least 0.005 in order to be really noticeable on the experimental curves. Regarding the coefficient β, changes whose amplitude is lower than 0.05 do not change results noticeably. The effect of the coefficient γ is the most subtle. Changes below 0.5 do not significantly affect the response. Any change of γ affects most shear specimens that have undergone a high axial damage. Those tests get damaged too quickly if the value of γ is too low. 64 4 - NUMERICAL MODELING OF THE SEISMIC BEHAVIOUR OF TILT-UP BUILDINGS Using the numerical EM5 connection model developed in Chapter 3, a detailed finite element program was developed to study the nonlinear response of tilt-up buildings. The nonlinear dynamic behaviour of this model was validated statically and dynamically against hand calculation, previous tilt-up building models, and a second numerical model. 4.1 Description of the finite element model The model is aimed at understanding the seismic behaviour of a box-type tilt-up building using solid panels. In order to use the coupled connection model presented in Chapter 3, the software MATLAB (MATLAB R2010a, MathWorks, 2010) was used. 65 4.1.1 Modeling assumptions A simplified 2D tilt-up buildings model was developed using the finite element method. Here are the main assumptions that were made: • All the tilt-up panels have the same dimensions. • The in-plan panels are assumed to be rigid and adequately braced in the out of plane direction. • The roof diaphragm is assumed to behave elastically and is hence modeled using a calibrated elastic spring. • The in-plane panels are not connected directly to the out-of-plane panels. • All the in-plane panels are linked using a steel angle at the top. As a result, all the panels have the same horizontal displacements. The steel angle also served as a rigid element to transfer the lateral force evenly to all the in-plan panels. It is assumed that this steel angle and its connections to the panels do not fail. • Since all the panels have the same horizontal displacement, the gap between the panels remains the same. This gap, which is typically about 1 or 2cm, is too small to allow a different rocking for each panel. Therefore, all the panels are assumed to have the same rocking angle. • All the panels can move vertically. Uplift of panels is made possible by the panel-panel connectors. • As a result the model has a total of 3+n degrees of freedom, where n represents the number of panels. Those degrees of freedom are called sliding, rocking, panel uplift, and roof translation. • Both contact and Coulomb friction between the slab and the panels are modeled. • Standard EM connectors are used to connect the panels to the slab and to connect the panels to each other. 66 4.1.2 Modeling of the roof diaphragm The seismic behaviour of the roof diaphragm is complex. Development of a detailed nonlinear model of the roof diaphragm is beyond the scope of this thesis. In this thesis, the roof diaphragm is modeled using an elastic spring and a lump mass as shown on Figure 4.1. The calibration of this model is explained in part 4.3.2. Figure 4.1 Simplified roof model 4.2 Mathematical Model 4.2.1 Static model Sections 4.2.1.1 to 4.2.1.6 describe successively the degrees of freedom of the system, the derivation of the stiffness matrix, the description of the contact forces, the equilibrium equation, the algorithm developed to determine displacement and contact forces, and finally the way equilibrium is insured. 67 4.2.1.1 Definition of degrees of freedom Figure 4.2 shows the free body diagram of a two panel system having two panel-slab connections and two panel-panel connections. The mechanical equations of the software are explained using this example, but similar approaches can be derived for higher number of panels. In static analysis, the roof is not considered since it can be readily replaced by a static force. Figure 4.2 Free body diagram of a two panel system In this figure: • wp is the panel weight • S1 and S2 are the applied lateral loads on panel 1 and 2. These loads can represent wind or earthquake forces. • Yc1 and Yc2 are the contact forces applied by the slab on panels 1 and 2. • Fr1 and Fr2 are the frictional forces at the panel-slab interface. 68 When the panels are rocking, the lower corner of each panel is referred to as the contact corner. Since the panels are assumed to be rigid, the external forces on each panel can be simplified to a resultant force and moment at the contact corner, as shown on Figure 4.3. The panels displacement is described by the motion of the contact corner. Figure 4.3 Conversion of external forces at the contact corners On this figure: • Sext1 and Sext2 are the horizontal external forces on panel 1 and 2, respectively; • Mext1 and Mext2 are the external moments on panel 1 and 2; • Yext1 and Yext2 are the vertical external forces on panel 1 and 2, respectively; • s is the horizontal displacement of the contact corners, called sliding • θ is the rotation angle of the panels, called rocking angle • y1 and y2 are the vertical displacements of the contact corner of the panel 1 and 2, they will be often referred to as uplift of panel 1 and 2. Since this system has four degrees of freedom, only four external forces are needed to solve the associated equilibrium matrix problem. The parameters Mext and Sext are chosen to describe the sliding and the rocking motion. They are defined using Eq 4.1 and 4.2: 69 Sext = Sext1 + Sext2 (4. 1) Mext = Mext1 + Mext2 (4. 2) Eq. 4.3 shows the displacement and force vectors: = 2 1 ext ext ext ext ext Y Y M S X and = 2 1 y y s x θ (4. 3) Where: • Xext is the external force vector. It is important to not that this vector includes both the applied loads and the contact forces; • x is the displacement vector; 4.2.1.2 Stiffness matrix The global stiffness matrix for the free body diaphragm shown in Figure 4.3 is computed by multiplying three matrices, as shown in Eq. 4.4. K = B k A (4. 4) Where: • A is the transformation matrix that converts the panel displacement (global) into connection displacements (local). Linear geometry transformation is used; • k is the stiffness matrix that represents the connector’s degrees of freedom. It converts the local displacement into local forces ; • B is the matrix that converts the local forces from each connection into global forces. The local degrees of freedom are shown on Figure 4.4: 70 Figure 4.4 Local degrees of freedom In this figure, the displacements xLi,, yLi, are the displacements of the ith connector. Because of the linear geometry used in the computation, the panel-panel connectors do not experience any axial elongation. Local displacements are computed from the global displacement vector using matrix A, as shown in Eq 4.5. = 2 1 8 7 6 6 5 5 4 3 2 2 1 1 y y s A y y y x y x y y y x y x L L L L L L L L L L L L θ (4. 5) 71 In the case of panels rocking to the right (i.e. θ≤0), matrices k, A and B for the two panel are computed using Eq 4.6 to 4.8. − − − − = pppp pppp ax sh ax sh pppp pppp ax sh ax sh kk kk k k k k kk kk k k k k k (4. 6) − − − − − − = 100 100 1020 0001 1010 0001 0100 0100 0120 0001 0110 0001 d d d d d d A (4. 7) tAB = (4. 8) Where: • ksh and kax are the shear and axial stiffness of the bottom connections, respectively; • kpp is the shear stiffness of the panel-panel connection; 72 • d is the panel width; • d1 and d2 are the distance from the contact corner to bottom connectors 1 and 2, respectively. Determination of rocking case It is important to note that the determination of the global stiffness matrix depends on which way the panels are rocking. Figure 4.5 shows the two possible rocking cases. Figure 4.5 The two rocking cases and their coordinates At every time step, the rocking angle is assumed to be the same as in the previous time step. If the computed rocking angle is compatible with this assumption, then the angle is confirmed. If the angle is opposite to the assumed angle, then the angle is set to zero. If the panel was not rocking at the previous time step, the rocking case is determined by the direction of the applied forces. Nonlinear connectors The stiffness of a nonlinear material depends on whether the material is loaded or unloaded. In order to determine whether each connector is being loaded or unloaded, the following procedure is used: 73 1- A trial stiffness matrix is calculated using an average of the loading slope and the unloading slope for each connector. 2- A trial displacement is calculated using this matrix. 3-.This trial displacement is used to derive the displacement direction of each connector. This information is used to calculate the stiffness matrix with correct slopes. 4- The displacement is calculated using the adequate stiffness matrix. This method may not work in some other problems. However, in this model, it was found to give the correct stiffness in all the cases encountered. 4.2.1.3 Contact forces The external force vector in Eq 4.3 contains both the applied loads and the contact forces: CFXX ext += (4. 7) Where: • X is the applied load vector (e.g. weight, wind, earthquake) • CF is the contact force vector. The contact force vector can be decomposed as follows: = 2 1 Yc Yc Mc Fr CF (4. 8) Where: • Fr is the Coulomb frictional force; • Mc is the contact moment. The value of this moment is zero when the panels are rocking. • Yc1 and Yc2 are the vertical reaction forces of the slab on panel 1 and 2, respectively. The contact force vector often contains some unknown parameters. This is why Eq 4.7 is shown to decompose the external force vector. 74 4.2.1.4 Equilibrium equation At every time step the equilibrium equation must be satisfied: extconnectors Xf = (4. 9) Where connectorsf is the internal force vector due to the connectors. Internal forces are approximated at every time step using Eq 4.10: )(*)()1( 1 nnnconnectorsconnectors xxKnfnf −+=+ + (4. 10) Where: • )(nf connectors is the internal force vector at time step n. This force is zero at the first time step: 0)0( ==nf connectors ; • Kn is the stiffness matrix at time step n; • xn is the displacement vector at time step n. Introducing the equilibrium equation and Eq 4.7 into Eq 4.10, Eq 4.11 is derived. )()(* 111 nnnnnnn CFXCFXxxK +−+=− +++ (4. 11) 4.2.1.5 Contact algorithm Eq 4.11 contains two unknown variables: the displacement and the contact forces at time step n+1. An algorithm has been written to solve this problem. Its outputs and inputs are summarized on Figure 4.6: Figure 4.6 Inputs and outputs of the contact algorithm 75 The main task of this algorithm is to answer the following questions: • Do the panels slide? • Do the panels rock? • Are any of the panels uplifted? In order to answer those questions, a procedure of successive assumptions and checks has been developed and implemented in the contact algorithm. Eq 4.11 is rewritten as Eq 4.12 to show the unknown variables. )( 2 1 2 1 1 1 1 1 1 1 1 1 1 nn n n n n nn n n n n n CFX Yc Yc Mc Fr Xx y y s K +− += − + + + + + + + + + θ (4.12) At time step n+1, sn+1, θ n+1, y1 n+1, y2 n+1, Fr n+1, Mc n+1, Yc1 n+1, and Yc2 n+1 are unknown. To solve Eq 4.12 for those 8 unknown variables, 4 assumptions have to be made. Assumptions about y1 n+1, y2 n+1, Yc1 n+1, and Yc2 n+1 are made by the uplift algorithm; assumptions about rocking and sliding are made by the rocking and the contact algorithm, respectively. The functioning of each of those three algorithm is now described. Uplift algorithm A first algorithm has been developed in order to determine whether some panels are uplifted or not. This algorithm requires two assumptions: one about the sliding and one about the rocking. Each of this assumption provides a numerical value for either the force or the displacement. For example, if the assumptions "No sliding" and "Rocking" are made, then the sliding displacement is known (sn+1 = sn), and the contact moment is known (Mc n+1=0). The inputs and outputs of this algorithm are described on Figure 4.7. 76 Figure 4.7 Inputs and outputs of the uplift algorithm Case 1: No uplift: y1n+1=0 and y2n+1=0 The uplift algorithm starts by assuming that no panel is uplifted: y1n+1=0 and y2n+1=0. Under those assumptions, there are four unknown variables: Yc1 n+1 ,Yc2 n+1, {sn+1 or Fr n+1} and {θ n+1 or Mc n+1} . Eq 4.12 is then solved for those four unknown variables. The assumption of no uplift is checked using the sign of the vertical reactions: • If Yc1n+1<0, it means that the vertical contact reaction on panel 1 is downwards, which is impossible since contact forces have to be upwards. So the assumption that panel 1 is not uplifted was false: the algorithm goes to case 2. • If Yc2n+1<0, then the non-uplift assumption on panel 2 was incorrect: the algorithm goes to case 3. • If Yc1n+1≥0 and Yc2n+1≥0, then no panel is uplifted, hence assumption of case 1 was correct. Therefore the algorithm returns the displacement and the contact force as outputs. Case 2: Uplift of panel 1: Yc1n+1=0 and y2n+1=0 In this case, panel 1 is uplifted, so there is no vertical contact forces on it. However, its vertical displacement y1n+1 is unknown. In case 2, the unknown variables are: y1 n+1 ,Yc2 n+1, {sn+1 or Fr n+1}, and {θ n+1 or Mc n+1}, , Eq 4.12 is solved for those unknown variables. It is assumed that the panels cannot be all uplifted at the same time, so the calculated displacement is not checked in this algorithm. As a result, the algorithm returns the displacement and the contact forces. 77 Case 3: Uplift of panel 2: y1n+1=0 & Yc2n+1=0 This case works exactly as case 2, the vertical displacement of panel 2 being unknown this time. All the steps of the uplift algorithm are summarized on Figure 4.8. A similar procedure is used if the system contains more than two panels. The uplift algorithm is then used inside a second algorithm called the rocking algorithm. Figure 4.8 Functioning of the uplift algorithm on a 2 panel example 78 Rocking algorithm The rocking algorithm can calculate displacement and contact force using only one assumption about the sliding condition: either the frictional force or the sliding displacement has to be known. Successive trials and checks are performed on the rocking conditions. For every assumption, the displacement and contact forces are calculated using the uplift algorithm. The inputs and outputs of the rocking algorithm are summarized on Figure 4.9. Case 1: Rocking - Mcn+1=0 In this first case, the rocking algorithm assumes that panels are rocking, which means that there is no moment due to the contact forces: Mcn+1=0. Using this assumption plus the input sliding displacement or force, displacement and contact forces are computed using the uplift algorithm. Then, the rocking angle is checked: • If the sign of the rocking angle is consistent with the rocking case (whose derivation was explained in section 4.2.1.2), then the rocking assumption was correct, and the algorithm returns the computed displacement and forces. • If the sign of the rocking angle is opposite to the rocking case, then rocking does not occur. The algorithm proceeds to case 2. Case 2: No rocking - θ n+1=0 In this case, displacement and contact forces are computed using the uplift algorithm with the condition θn+1=0. These displacement and forces are finally return as outputs by the rocking algorithm. These steps are summarized on Figure 4.10. 79 Figure 4.9 Inputs and outputs of the rocking algorithm Figure 4.10 Functioning of the rocking algorithm 80 Contact algorithm The contact algorithm makes successive sliding assumptions in order to determine whether sliding occurs or not. The rocking algorithm is used in each step of this process, and the input does not necessitate any assumption, as shown on Figure 4.6. Case 1: Assume no friction - Frn+1 =0 This first step is aimed at determining in which way the panels tend to slide. To do so, the sliding displacement is computed under the assumption that there is no frictional force. As a result, the frictional force Frn+1 is known in case 1 (Frn+1 =0), but the sliding displacement sn+1 is unknown. The sign of the sliding displacement under this assumption is computed using the rocking algorithm. Case 2: Assume sliding - Frn+1 is known In this case, panels are assumed to slide. Under this assumption, the frictional force is known. Its direction is the opposite direction of the sliding in the no-friction case (from the case 1), and its value is computed using Eq 4.13: Fr= )( 1 1+nssign *µ*neg(Y1n+1+Y1n+1) (4. 13) Where : • 1 1+ns is the sliding displacement given by the first step. • neg(Y1n+1+Y1n+1) is the negative part of the sum of the vertical loads applied on panel 1 and panel 2 at time step n+1. This sum will typically contain panel weight and the vertical seismic forces. The rocking algorithm is then used to derive the sliding displacement, written 2 1+ns . The sliding assumption is checked as follows: • If 0* 2 1 1 1 ≥++ nn ss , then the sliding assumption was correct. Indeed, this equations says that the panels slide in one direction whether there is a frictional force or not. The contact algorithm is then finished, it can return the displacement and contact forces of the sliding case. • If 0* 2 1 1 1 <++ nn ss , then panels are not sliding. Then the algorithm goes to case 3. 81 Case 3: Panels are not sliding - sn+1 = sn In this case, the rocking algorithm is called again with the input sn+1 = sn. The outputs of the rocking algorithm - displacement and contact forces, are then returned as the outputs of the contact algorithm. Figure 4.11 summarizes all the successive steps of the contact algorithm. Figure 4.11 Functioning of the contact algorithm 82 4.2.1.6 Equilibrium check At any time step n+1, Eq 4.12 is solved using the tangent stiffness at the previous displacement xn. This assumption is not exact when some nonlinearities occur in the resisting force connectorsf . To solve this problem, the external force is set equal to the resisting force at the calculated displacement xn+1 at the end of the time step. Contrary to typical algorithms, the displacement is not recalculated at the current time step. The following example explains th method using a one degree-of-freedom example. Figure 4.12 Modification of external force for equilibrium check Figure 4.12 shows the resisting force of the system as a function of the displacement. At step one, the displacement is x1, and the external force is Xext1. At step two, the applied force is Xext2. Since the displacement is computed using the tangent stiffness at the displacement x1, the equilibrium equation Eq 4.12 returns the displacement x2, whereas the real displacement for this external force is x2target, which is higher. Since, the internal force is lower than Xext2, the external force Xext2 is modified by the program to force equilibrium. Thus, the external force at time step 2 is not the force that was really applied. At step 3, the stiffness at the displacement x2 is used. Since the starting force is the "Modified Xext2", the right displacement is obtained at time step 3. As a result, even though the force 83 computed at one time step is not exact, the error does not build up, and the correct forces and displacement are found at the next time step. A time step sufficiently small should be selected to avoid having two stiffness changes at two successive time steps. It is important to note that this method would not work if the system behaviour was not piecewise linear. The system behaviour is piecewise linear because all the material behaviours are piecewise linear. 4.2.2 Dynamic model 4.2.2.1 Mass matrix As the center of mass of the panels is located at their centroid, it is easier use the degrees-of- freedom shown on Figure 4.13 for dynamic analysis. All the terms that are calculated using the middle of the panels as a reference have an m superscript. Figure 4.13 Degrees of freedom defined at the middle of the panels 84 4.2.2.2 Damping matrix The damping matrix is assembled using a 3% Rayleigh damping assigned to the first two modes. The periods of the first two modes are assumed to be Tr and 0.3*Tr where Tr is the roof fundamental period. Eq 4.14 shows the general Rayleigh damping equation. C = a0.M +a1.K (4. 14) Where a0 and a1 are constant coefficients which are computed using Eq 4.15 and 4.16, respectively: 21 21 0 2 ωω ωωζ + =a (4. 15) 21 1 2 ωω ζ + =a (4. 16) Where ζ is the damping ratio and ω1 and ω1 are the pulsation corresponding to Tr and 0.3*Tr, respectively. 4.2.2.3 Dynamic algorithm The Newmark β-γ method is used to calculate the displacements in dynamic analysis. This method assumes that the displacement and the velocity at the step n+1 can be approximated using the following formulas: 1 22 1 11 ))5.0(( ))1(( ++ ++ +−++= +−+= nnnnn nnnn xhxhxhxx xhxhxx &&&&& &&&&&& ββ γγ (4. 17) The motion of the system is ruled by the equation of motion: m n m n m n m n mm n m CFXfxCxM 11111 +++++ +=++ &&& (4. 18) Using the Newmark equations in the equation of motion, Eq 4.17 is derived. 85 m n m n m n m n m n m n m n m m n m n m n m nm CFX fxhxxx h C x h xhxxM 11 11 2 1 ) 2 1()1())(( )1 2 1( ++ ++ + += + −+−+−+ −− −− &&& && & β γ β γ β γ ββ (4. 19) Using the fact that: m n m n m n m n mm n m CFXfxCxM +=++ &&& and mnmnmnmnmn ffxxK −=− ++ 11 )(* The following equation is derived: )()))1 2 (() 2 (( ).(1 11 12 m n m n m n m n m n m n m nm n m n m n m n mm CFXCFxhxCx h xMX xxKC h M h +−+−++−+ =− ++ ++ + &&& &&& β γ β γ ββ β γ β (4. 20) Contact is easier to deal with using corner coordinates. In order to use those coordinates, the conversion matrices T and R were defined as follows: • T converts force at middle into force at the corner: mXTX .= (4. 21) • R converts displacement of the middle into displacement of the corner: xRxm .= (4. 22) Eq 4.20 was multiplied by T on the left side, and corner displacement was introduced: )..()))1 2 (() 2 (.( )...(1. 11 12 m n m n m n m n m n m n m nmm n nn m n mm CFTXTCFxhxCx h xMXT xRxRKC h M h T +−+−++−+ =− ++ ++ + &&& &&& β γ β γ ββ β γ β (4. 23) This equation can be written: )())1 2 (() 2 (. ).(....1 11 12 nnn m n m n m n m nmm n nnn mm CFXCFxhxCx h xMXT xxKRCT h RMT h +−+ −++−+ =− ++ ++ + &&& &&& β γ β γ ββ β γ β (4. 24) 86 This last equation has the same form as Eq (4.11) which is the equation being solved by the contact algorithm. As a result, the contact algorithm can be used in dynamics too. The proper input has to be used, as shown on Figure 4.14: Figure 4.14 Use of contact algorithm in dynamics 4.2.2.4 Panel pounding In dynamics, collision between the panels and the slab occurs every time the panels rock from one way to the other. This phenomenon, called panel pounding, results in a damping of the rocking motion. The real amount of dissipated energy depends on the materials, the foundation and the soil type. A precise rocking model would require a specific attention on the pounding model. The tilt-up model by Olund (2009) described in section 2.4.1 assumed that this pounding is elastic: no energy is dissipated when the panels pound the ground. In this study, the functioning of the contact algorithm compels to make the following assumption: when the panels rock from one way to the other, their vertical and rotational kinetic energy are lost, but the horizontal kinetic energy is conserved. 87 4.3 Model calibration 4.3.1 Typical building geometry As written in part 2.4.1, Olund (2009) conducted a survey amongst tilt-up designers in both Canada and the U.S in 2009 in order to get the most common characteristics of a concrete tilt- up building. He used this study to design two typical tilt-up buildings: one using solid panels, and one using panels with openings. The design he got with solid panels is reused in the analytical part of this study. Here are the characteristics of this design: • The building is 30.48m wide, 60.96m long, and 9.144m (30ft) high. The perimeter of the building is made of two rows of four panels in the short direction, and two rows of eight panels in the long direction. Panel thickness is 0.184m. This thickness is chosen for two reasons: Firstly it allows the slenderness ration of panels to be right below 50, and secondly 0.184m is the actual width of 2"x8" piece of lumber, thus making easy formwork manufacturing. Panel width is 7.62m. This makes the weight of the panel just below 30,000kg, which is the panel weight limitation in many areas in order to limit the sizes of cranes. • Regular connectors EM1~5 are used to connect concrete panels to the roof diaphragm, to the slab, and to each other. Three panel-slab connection and two panel-panel connections per panel are used in the short direction. Two panel-slab connections per panel and no panel-panel connections are used in the long direction. The calculation of the number of required connectors in the short direction is explained in Appendex A. • A steel roof diaphragm is used. It is made of the four following components: - A steel angle L102x102x13 going all along the perimeter 88 - A W610x195 beam spreading at midwidth through the whole length of the building. Its midspan is supported by an HSS254x254x9.5 column - 62 - 1050mm deep OSWJ with L/L chord members, which cross the short direction from the panels to the middle beam. - A steel deck Type B . It is 38mm deep, and its thicknesses is 18ga or 20ga depending on the area. Side lap fasteners are screws, and deck-to-frame fasteners are nails. Figure 4.15 shows the general layout of half of the building: Figure 4.15 General layout of the typical building used for the analysis (Olund, 2009) 89 The weakest direction of this building is the short direction. So this direction was modeled. A value 0.5 was chosen for the concrete-concrete friction coefficient (as recommended by the Concrete Design Handbook - Third Edition). The paired connection models described in section 3.3.5 were used to model the building connections. Panel-panel connections were modeled using the EM5-EM5 model, and panel- slab connections were modeled using the EM5-EM4 model. In practice, the EM5-EM3 connection is much more commonly used to link the panels to the slab. However, using the EM5-EM4 model enables to study the behaviour of buildings where the whole capacity of the EM5 connector is used. The seismic intensity at which the bottom connectors would fail if EM5-EM3 connectors were used can still be seen in the results since those two models have the same elastic stiffness. 4.3.2 Roof diaphragm model calibration To calibrate the stiffness of the roof diaphragm, a 2D finite element model was developed using the software OpenSees (version 2.2.2, McKenna et al., 2011). This model is shown in Figure 4.16. 90 Figure 4.16 OpenSees 2D roof model All the truss elements are pin-connected to each other. The girder and the angle are modeled by beam-column elements. A 200GPa steel was used for each element. In order to model panel restraint, the left-end and right-end node columns are restrained vertically, and the top and bottom node rows are restrained horizontally. Half of the mass of the out-of-plane panels is distributed along the top and bottom nodes rows. The diaphragm shear behaviour is modeled using diagonal bars. The stiffness of those bars is set to fit a shear stiffness of 3.0kN/mm, following the experimental test results from Essa et al. (2003) (Tests 4 and 7 on Table 2.3 and 2.4). This calibration was done by simulating the 91 experimental test shown on Figure 2.13 in the finite element software, as shown on Figure 4.17. In order to model experimental conditions, bottom nodes are clamped, and the top layer of nodes has a rigid body constraint. Figure 4.17 Simulation of the test using Opensees Table 4.1 shows the mass, period and stiffness of the roof model that were obtained from modal analysis. It can be noted that the roof period is different from the one obtained in Olund (2009) where a similar roof diaphragm was modeled and a period of 0.58s was estimated. The difference is due to the stiffness of the diagonal elements. In the proposed model, the stiffness of the diagonal elements was set by reproducing the test results of Essa et al. (2003) with the roof model, whereas in the model of Olund (2009), an equation obtained from the Steel Deck Institute (SDI) was used. Table 4.1 Roof model parameters Period (s) 1.04 Mass (metric ton) 205 Stiffness (kN/m) 7500 92 4.4 Model Validation The model was checked using hand calculation, comparison with previous tilt-up numerical models, and a numerical model developed on another finite-element software. 4.4.1 Hand calculation i) Static behaviour Static checks were performed on two configurations: • Configuration 1: one single panel with two panel-slab connectors, at the top of which were applied a positive and a negative loading cycle using the following force sequence: 0 →155kN → -155kN → 0 • Configuration 2: a two-panel system with two panel-panel connectors and one panel- slab connector. The following force sequence was applied at the top of the panels: 0 →575kN → 0 For each configuration, panel dimensions were 9.14m x 7.28m x 0.18m, and the panel weight was 288kN. Figure 4.18 shows the sliding and the rocking displacement of configuration 1, and Figure 4.19 shows the sliding, rocking and uplift (of the left panel) displacements of configuration 2. Points that were checked by hand calculation are described by the black squares. In each of these plots, the x-axis represents the force applied at the top of the panels. The rocking drift is defined as the displacement of the top of the panels due to the rocking motion. 93 (a) Sliding; (b) Rocking Figure 4.18 Displacement of configuration 1 Figure 4.19 Displacement of configuration 2 94 Sliding behaviour Using Coulomb frictional model, a minimum force is needed to trigger sliding because of frictional resistance. This force is given by: µ** wpnpFslidehc = (4. 25) Using a value of 0.5 for the friction coefficient, the minimum force needed to trigger slilding is, for each configuration: kNFslide hcConfig 1441 = kNFslide hcConfig 2882 = As shown on Figure 4.18 (a) and 4.19, the numerical model gives: kNFslide elConfig 144 mod 1 = kNFslide elConfig 288 mod 2 = If the force is increased beyond the onset of sliding, then the stiffness is the sum of the shear stiffness of all the panel-slab connectors. In the linear range of those connectors, the stiffness is: sh hc knpKslide *2*= (4. 26) Where np is the number of panels and ksh is the shear stiffness of the bottom connectors. Computing Eq 4.27, the following sliding stiffness are derived for configurations 1 and 2: mmkNKslide hcConfig /3041 = mmkNKslide hcConfig /6082 = As shown on Figure 4.18 (a) and 4.19, the sliding stiffness in the numerical models are: mmkNKslide elConfig /304 mod 1 = mmkNKslide elConfig /609 mod 2 = Then, when the force is decreased and applied in the opposite direction, the panel should not move in a first time. The frictional force that was resisting positive displacement is now 95 resisting negative displacement. The force required to trigger sliding in the negative direction is: µ***2max wpnpFFslide neghc −= (4. 27) Where Fmax is the maximum positive force. In configuration 1, the maximum applied force was 155kN. So Eq 4.27 gives: kNFslide neghc 133−= On Figure 4.18 (a), sliding in the negative direction starts at: kNFslide negel 133mod −= As a conclusion, the sliding curves from Figure 4.18 (a) and 4.19 are accurate: the sliding behaviour is well described by the model. Rocking behaviour In this section, the rocking behaviour of the model, as shown on Figure 4.18 (b) and 4.19, is checked. First, the force needed to overcome the moment due to the weight and thus trigger rocking is given by: H d wpnpFrock hc *2 **= (4. 28) For each configuration, this force is: kNFrock hcConfig 1141 = kNFrock hc Config 2292 = The proposed model gives: kNFrock elConfig 114 mod 1 = kNFrock el Config 229 mod 2 = This demonstrated that rocking motion starts at the correct force. In order to get the rotational stiffness of configuration 2 using hand calculation, the following free body diagram was used: 96 Figure 4.20 Two panel model for hand calculation The fact that both panels are compelled to have the same sliding and the same rotation is modeled by the rigid sliding frame in the upper part. The applied force F is applied at the top of this frame, and thus is redistributed to the two panels - namely F1 on panel 1 and F2 on panel 2. This can be written: FFF =+ 21 (4. 29) Since the frame is rigid, panel 2 applies a moment on panel 1, which is called M1. Respectively, panel 1 applies a moment -M1 to panel 2 because there is no shear force in the horizontal beam of the frame. The moment equilibrium at the bottom right corner gives: For panel 1: 11 22 *)) 4 1() 4 3((2/* MHFddkdw axp +=++ θθ (4. 30) For panel 2: 212 22 2*)) 4 1() 4 3((2/* dkMHFddkdw ppaxp θθθ −−=++ (4. 31) Where: • kax is the axial stiffness of the bottom connections; • kpp is the shear stiffness of the panel-panel connections; • θ is the rocking angle; 97 • H and d are the panel height and width, respectively. By adding Eq 4.30 and 4.31 and using Eq 4.29 to get rid of F1+F2, Eq 4.32 is derived: dwHFkkd pppv **)24 5(2 −=+θ (4. 32) So, if a force of 350kN is applied, the rocking angle is: mradhc 121.0=θ In the proposed model, the rotation at a force of 350kN is: mradHmmel 123.0/12.1mod ==θ The rotational stiffness is thus correct. As a result, the rotational behaviour is well described by the model, both qualitatively and quantitatively. Uplift behaviour In configuration 2, the vertical equilibrium of panel 1 is written: 0_2)) 4 1() 4 3(( =+++−− forceContactdkddkw ppvp θθθ (4. 33) The panel is lifted when the contact force is down to zero. This condition is fulfilled when: 02 =+−− dkdkw ppvp θθ (4. 34) Therefore, the angle at which the uplift of panel 1 occurs is: mrad dkdk w ppv phc uplift 93.22 = +− =θ As shown in Figure 4.19, uplift occurs at the following angle in the model: mradHmmeluplift 90.2/65.2mod ==θ Therefore uplift occurs at the correct step. 98 As a conclusion, the nonlinear cyclic behaviour of the building model is validated in the cases where the EM5 connectors are in their linear range. The validation of the building model in the nonlinear range of the connectors is shown in section 4.4.2. ii) Vibrations Period of vibrations of a pure sliding motion and a pure rocking motion were checked. Those two checks were performed on a 4-0-2 configuration without roof. Pure sliding The theoretical sliding period is: s k Mn Ts sh p 077.0 8 2 == pi (4. 35) Where: • M is the mass of one panel • ksh is the shear stiffness of the bottom connectors. Modal analysis on the MATLAB model gives a period of 0.076s, which is very close to the theoretical period. Pure rocking The theoretical rocking period is: s ddk JTr ax 30.0 4 3 4 4 42 22 = + = pi (4. 36) Where: • J is the second moment of inertia of mass of one panel, calculated at the rocking corner; • kax is the axial stiffness of bottom connectors. 99 Modal analysis on the MATLAB model gives a period of 0.31s. As a conclusion, the period of the rocking and sliding motions is accurate. Hence, the response rate of those two motions is realistic. 4.4.2 Comparison with model of Devine (2009) The pushover behaviour of the tilt-up building model was compared with the one obtained by Devine (2009). In the pushover analysis from Devine, a load was applied at the top of a panel row. Various panel and connector configurations were studied. In the current study, only four panel systems are studied. In Figure 4.21, some pushover plots from Devine are compared with pushover plots obtained from the proposed model in a similar configuration. The connector configuration is described by three digits: the first one is the number of panels, the second one is the number of panel-panel connectors, and the third one is the number of panel-slab connectors. The x-axis is the drift, which is the displacement of the top of the panel divided by the panel height. 100 Figure 4.21 Comparison of pushover behaviour of the model with F.Devine's model (2009) for five configurations 101 The pushover curves match well with the previous ones obtained from Devine (2009). The minor differences may be due to the differences in the connection model. The most noticeable difference occurs in the 4-3-2 configuration. The model of Devine failed in a sliding mechanism, whereas in the proposed model, no sliding failure is experienced because of the higher shear strength of the bottom connectors. Those differences are sufficiently small to consider this comparison as a validation of the static monotonic pushover behaviour of the model. In particular, the implementation of the nonlinear EM5 connection model in the building model is working properly. 4.4.3 Comparison with another numerical model As explained in 4.2.2, a significant amount of rocking kinetic energy is lost every time the panels pound the foundation. This might affect the rocking demand on the panels. In order to check the dynamic rocking behaviour of the panels, a simple numerical model was developed using the software Perform 3D, Version 4.0.3 (Computers and Structures, 2007). The geometry of the model is shown on Figure 4.22. It consists of the following elements: • One single rigid panel, with a lumped mass at its center. • A mass-spring system attached at the top of the panel to represent the roof. The mass and stiffness were chosen as a quarter of the ones of the roof model described in section 4.3.2 in order to have the same panel loading as in the typical building where the roof load is evenly distributed to four panels. • A spring and dashpot system at each bottom corner of the panel. The aim of those elements is to represent the contact with the foundation slab: a 'compression only' stiff spring is attached to the corner, the bottom of this spring is attached to a spring-dashpot system which models the energy loss during collisions. The bottom of the spring- dashpot system is fixed. • A horizontal restraint of the panel bottom, aimed at preventing sliding. • A 3% two-mode Rayleigh damping was implemented at the periods of 1.04s and 0.31s. 102 Figure 4.22 Perform 3D panel model There is no tilt-up connector in this model. Two values of the damping coefficient C were tried: C=0 and C=180 kN.m-1.s. C=0 corresponds to an undamped foundation. C=180 kN.m-1.s represents a foundation with very high damping where nearly all the rocking energy is lost when the panel pounds the slab. This Perform 3D model was compared to the proposed model with four panels, a full roof and no connector. The two models were subjected to a ground motion scaled to four different amplitudes. The chosen ground motion record is 'LOS000' from the Northridge earthquake (1994, Magnitude 6.7, recorded at Canyon County) , it was taken on the PEER-NGA database (http://peer.berkeley.edu/nga/); scaling factors were 1, 2, 5 and 10. 103 For each amplitude, the rocking displacement with and without foundation damping was compared to the one given by the proposed model. Figure 4.23 shows the rocking drift time-history for: • The proposed model; • The Perform 3D model without damping; • The Perform 3D model with connection damping. The reduction of the peak rocking drift of the undamped model is compared for the proposed model and the Perform 3D damped model in Table 4.2. Table 4.2 Peak rocking drift reduction Peak rocking drift reduction compared to the undamped model (in %) Ground motion scaling 1 2 5 10 Average Proposed model 61 28 57 25 43 Perform model with damping 72 65 72 41 62 104 Figure 4.23 Comparison of the proposed model with the Perform model using time history analysis with one ground motion record scaled at four amplitudes 105 Based on the time history analysis presented in Figure 4.23 and on the peak drift reduction presented in Table 4.2, it can be concluded that: • The rocking drift given by the proposed model is always lower than the undamped displacement given by the Perform 3D model. The peak drift reduction compared to the undamped Perform 3D model is 43% on average. No trend was observed regarding whether small or large intensities result in a high peak drift reduction ratio. • The rocking drift given by the model is always higher than the drift of the Perform 3D damped system. • Depending on scaling, the response of the proposed model is either closer to the damped model or to the undamped model. However, on average, the peak drift is closer to one given by the damped model. As a conclusion, the dynamic rocking behaviour of the proposed model is validated. It results in a rocking demand intermediate between an undamped model and a model with very high energy dissipation in the foundation. The reduction of the rocking drift demand compared to an undamped model is 43% on average, according to the four nonlinear time history analysis performed in this section. 106 5 - NONLINEAR RESPONSE OF TILT-UP BUILDINGS In this chapter, the nonlinear response of the typical building described in part 4.3 is studied. The building model developed in Chapter 4 is calibrated to this building, and the coupled connection behaviour developed in Chapter 3 is incorporated. The seismic performance of this prototype building is studied using Incremental Dynamic Analysis. In addition, a detailed parametric study is conducted to investigate the structural response of the building with different connection configurations, roof diaphragm stiffness and strength, and friction coefficient between the panels and foundation slab. As shown in sections 4.2.2.4 and 4.4.3, the dynamic rocking displacement given by the model is significantly decreased by the model used for the panel-slab contact, which dissipates a lot of rocking energy. Therefore, using the results presented in this chapter for design purpose might be unconservative, unless the soil and foundation of the building have a very high damping. 107 5.1 Seismic performance evaluation of the prototype building The seismic performance of the prototype building described in section 4.3 is analyzed in this section. The seismic performance was conducted using the Incremental Dynamic Analysis procedure proposed by Vamvatsikos and Cornell (2002). 5.1.1 Seismic hazard The prototype building is located on a Class C soil (2010 NBCC) in Surrey, British Columbia, Canada. The target spectrum is selected to represent the 5% damped response spectrum with 2% probability of exceedance in 50 years. A set of 10 ground motion records was selected from the PEER-NGA database (http://peer.berkeley.edu/nga/), such that the median spectrum of those records fits best the target spectrum in the period range between 0.2s and 1.5s. Table 5.1 shows the characteristics of the selected ground motions, and Figure 5.1 shows the 5% damped response spectra for the selected ground motions along with the target spectrum. Table 5.1 Summary of the selected ground motions NGA# Component Event Year Station Mag (Mw) Mechanism Rjb (km) Rrup (km) Vs30 (m/s) 265 FN Victoria- Mexico 1980 Cerro Prieto 6.3 Strike-Slip 13.8 14.4 659 1010 FP Northridge-01 1994 LA - Wadsworth VA Hospital South 6.7 Reverse 14.6 23.6 414 3473 FP Chi-Chi- Taiwan-06 1999 TCU078 6.3 Reverse 5.7 11.5 443 3473 FN Chi-Chi- Taiwan-06 1999 TCU078 6.3 Reverse 5.7 11.5 443 2734 FP Chi-Chi- Taiwan-04 1999 CHY074 6.2 Strike-Slip 6.0 6.2 553 1111 FN Kobe- Japan 1995 Nishi-Akashi 6.9 Strike-Slip 7.1 7.1 609 727 FN Superstition Hills-02 1987 Superstition Mtn Camera 6.5 Strike-Slip 5.6 5.6 362 125 FN Friuli- Italy-01 1976 Tolmezzo 6.5 Reverse 15 15.8 425 125 FP Friuli- Italy-01 1976 Tolmezzo 6.5 Reverse 15 15.8 425 739 FN Loma Prieta 1989 Anderson Dam (Downstream) 6.9 Reverse-Oblique 19.9 20.3 489 108 Figure 5.1 Design spectrum and spectra of the selected ground motions with 5% damping 5.1.2 Structural response The building response was studied using the Incremental Dynamic Analysis where the ground motion presented in Table 5.1 were scaled at an increasing seismic intensity. In this study, the intensity parameter is chosen as the spectral acceleration at the fundamental period of the building Sa(Tr), where Tr is 1.04s. Three engineering demand parameters have been used to study the building response: • The sliding drift: defined as the peak sliding displacement divided by the building height, which is 9.14 m. This drift is directly related to the shear demand on the bottom connectors: all the panel-slab connectors fail in shear when the sliding drift exceeds 0.18%. It is important to note that, if an EM5-EM3 connection is used to link the panels to the slab, shear failure of the bottom connections occurs at a drift of 0.02% only. • The rocking drift: defined as the peak rocking angle. This drift is directly related to the shear demand on the panel-panel connectors and the axial demand on the bottom connectors. At a rocking drift of 0.43%, most panel-panel connectors fail in shear, at 109 1.75% some panel-slab connectors fail in uplift, and failure of all the panel-slab connectors occurs between 2.6% and 5.2%. • The normalized roof shear demand: defined as the peak roof shear demand divided by the roof self weight (2010kN). The complete failure of the building is defined when all the panel to slab connections have failed. Though it is possible for the tilt-up building to resist gravity and lateral loads after all the panel to slab connections have failed, the only resistance left is provided by the frictional force and the self weight resisting overturning, this is not a desired mechanism. Figure 5.2 shows the sliding drift, rocking drift, and normalized roof shear obtained from the Incremental Dynamic Analysis. Figure 5.2 IDA results for the 10 selected ground motion records 110 As shown in Figure 5.2, the median sliding drift demand can be broken down into three parts: • If Sa(Tr) is lower than 0.18g, then panels do not slide. • If the value of Sa(Tr) is between 0.18g and 0.7g, the median response is linear and sliding displacement is very limited for the majority of the ground motion records. • If the value of Sa(Tr) is increased past 0.7g, then the sliding drift increases a lot. Shear failure of all the panel-slab connectors is reached at Sa(Tr) = 1.05g. The rocking drift follows the same pattern: • If the value of Sa(Tr) is below 0.12g, panels do not rock. • For values of Sa(Tr) between 0.12g and 0.45g, the rocking drift increases quasi linearly at low rate. • When the value of Sa(Tr) increases beyond 0.45g, rocking drift increases at a much higher rate Failure of the building occurred at a median spectral acceleration of 1.05g. At this intensity, the median rocking drift is 0.32%, which is much lower than 1.75% at which axial failure of some bottom connectors would occur. For two ground motions, some panel-panel connectors experienced failure at a lower intensity than the building sliding failure. If the EM5-EM3 model was used to model the bottom connectors, sliding failure would occur at a drift of 0.02%, which is reached at an intensity of 0.52g. This intensity is 50% lower than the intensity at which the bottom connectors fail with the EM5-EM4 model. Hence, it is concluded that the anchorage of the bottom connection in the panels is a critical aspect of the design. The shear demand on the roof diaphragm increases nearly linearly with the spectral acceleration. This shows that panel drift does not affect significantly the roof behaviour. This is due to the fact that the sliding and rocking drifts are very low. If the contact damping was reduced, higher overturning drifts would be reached, and the roof IDA curves would be less linear. 111 The design spectral acceleration at the fundamental building period is 0.34g. At this intensity, amongst the ten considered ground motions, only one results in yielding of the connectors - yielding being experienced in shear by the panel-panel connectors and in the axial direction by some panel-slab connectors. The same result is observed if EM5-EM3 connections are used between the slab and the panels. Therefore, the connection design is adequate according to this model. Four phenomena can explain this very low level of damage: • The rocking demand is significantly reduced in the model because of the contact model. Furthermore, since the rocking demand is low, the shear strength of the bottom connectors is not significantly affected by coupling. • The roof diaphragm is very flexible. This results to a small shear force, and thus a low sliding and rocking demand. If the roof period is reduced to 0.6s, the sliding demand will be 0.05% based on a panel height of 9.14 m, which is 2.5 times the yielding drift. • The shear strength of the bottom connectors is much higher in the model than the value usually chosen in design. 112 5.2 Parametric Study In this part, a parametric study is carried out using the building designed in part 4.3. The importance of the following parameters is investigated: • Connection configuration • Friction coefficient • Roof strength and stiffness • Connectors shear behaviour. 5.2.1 Connection configuration - Pushover analysis The nonlinear behaviour of the typical tilt-up building with various connection configurations was studied using static pushover analysis. A similar study was previously done by Devine (2009). In this study, the roof was replaced by a static load, as shown on Figure 5.3. Figure 5.3 Pushover loading 113 i) Failure sequence This first sub-section is made to provide readers that are not accustomed to tilt-up systems with a better understanding of the failure mechanisms of this kind of building. The sequence of failure of all the connections during a monotonic pushover load is sketched for configurations using four panels and two panel-slab connectors along with the corresponding pushover curve on Figure 5.4. In the figure, "nf" means that the connector did not fail during the analysis. In the configuration without panel-panel connector, all the panels behave similarly. At a drift of 1.75%, all the left bottom connectors fail in the axial direction because of the rocking motion. In the configuration with one panel-panel connector per panel, the strength is greatly increased compared to the previous configuration. Indeed, the panel-panel connectors make the rocking of the panels much stiffer. At a drift of 0.4%, all those connectors fail, and the system behave then exactly as the previous on. In the third configuration, the two panel-panel connections are strong enough to uplift the panel at the left end. As a result, at a drift of 0.4%, the panel-panel connectors in the middle and on the right fail, but the panel-panel connectors between the two left panels do not. The uplift panel increases significantly the overturning resistance of the system. As a result, the pushover force is higher than with the previous configurations at all drift. In the fourth configuration, there are four panel-panel connectors per panel. The overturning is increased so much that the overturning resistance of the system is higher than the sliding resistance. As a result, all the bottom connectors fail in shear simultaneously at a very small drift. 114 Figure 5.4 Connection failure sequence 115 ii) Pushover analysis on multiple connection configurations Figures 5.5 to 5.7 shows the results of the monotonic pushover analysis with 2, 3 and 4 panel to slab connectors, respectively. In each figure the number of the panel to panel connectors was varied from 0 to 4. Figure 5.5 Pushover curves for configurations using two panel-slab connectors per panel 116 Figure 5.6 Pushover curves for configurations using three panel-slab connectors per panel Figure 5.7 Pushover curves for configurations using four panel-slab connectors per panel 117 Based on the results presented in Figures 5.5 to 5.7, it is concluded that: • Increasing the number of panel-panel connectors increases significantly the system strength. However, most of this strength gained by adding panel-panel connectors is lost at very small drift (around 0.4%). • All the systems having the same number of panel-slab connectors have the same ductility, except configuration 4-4-2. This is due to the fact that nearly all the systems fail because of the overturning demand. • Coupling does not affect significantly these pushover curves. Indeed, coupling plays a significant role mainly at high rocking drifts. When rocking drift is high, the bottom connectors that have not failed yet are close to the contact corner, thus create a very limited resisting moment. As a result, the applied force becomes very low at high drifts, and is not sufficient to induce a shear failure of the bottom connectors, even though they have a reduced strength. Effect of coupling would be higher if: - Bottom connectors were located closer to the middle of the panel, thus requiring a higher force to induce overturning failure. - Only two or three panels were used. Indeed, in systems using two panels and a lot of panel-panel connectors, one panel gets uplifted. The increase of the overturning resistance due to this uplifted panel is the same as in a four panel configurations; however, the sliding strength is twice smaller, so the shear strength degradation is much more critical. • In systems having sufficient panel-panel connectors, one or two panels are uplifted. Uplifted panels are responsible for a significant increase of the overturning resistance of the system, which is a desirable characteristic. However, when panels are uplifted, their connection to the slab fails at a lower rocking drift. In monotonic pushover, this phenomenon does not affect the ductility since no configuration fail in sliding. However, in the case of a cyclic loading, configurations where uplift occur would probably have a lower ductility. • Pushover analysis on configuration 4-3-2 predicts an overturning failure. However, the building was found to fail in a sliding mechanism in the Incremental Dynamic 118 Analysis. Three phenomena can explain the difference between static and dynamic behaviour: - First, dynamic load is applied both at mid height and at panel top in the dynamic model, whereas in pushover analysis, the static load is applied only at the top of the panels. In this building, the weight of the in-plane panel row is only 1.6 times lower than the weight of the roof. As a result, a significant part of the load is applied at mid height, resulting in a much reduced moment. - Second, in the dynamic model, the earthquake forces are resisted by both the base (by the panel-slab connectors and friction) and the inertial forces. Contrary to base forces, inertial forces resist the lateral load at panel mid height. When the lateral load is resisted at mid height, a lower overturning moment is created. This dynamic phenomenon can reduce the overturning moment. However, the amplitude of this effect was not studied. - Third, the contact condition used in this study results in a large dissipation of the rocking kinetic energy. A similar pushover study was performed in Devine (2009). As shown in section 4.4.2, pushover plots from the model are very close to the ones obtained by Devine. Most of the conclusions from Devine are in agreement with the proposed model. The main difference lies in the failure mode of some configurations. In Devine (2009), configurations 4-3-2, 4-3-3, 4-4- 4 and 4-4-3 failed in a sliding mechanism, whereas those configurations failed in an overturning mechanism in the proposed model. The reason is that the proposed connection model has a higher shear strength than the model from Devine. 119 5.2.2 Connection configuration - Nonlinear dynamic analysis Using the set of ground motions presented in Table 5.1, the connection configuration for the prototype model are studied using incremental dynamic analysis. To limit the computational analysis, only the following configurations were studied: configurations with 3 panel-slab connections per panel and 0 to 4 panel-panel connections, and configurations with 2 panel-slab connectors per panel and 2 to 4 panel-slab connectors. Figure 5.8 shows the results of the Incremental Dynamic Analysis. The vertical axis represents the shaking intensity described by the spectral acceleration at the fundamental period of the building, Sa(Tr). The horizontal axis represents the median peak response obtained from the ten earthquake records. Figure 5.8 Demand parameters for various connector configurations 120 The following observations can be made: • Increasing the number of panel-panel connectors reduces significantly the overturning demand. However, it also increases the roof shear by making the panels-connections system stiffer. This increase in the roof shear is responsible for an increase in the sliding demand. As a result, in this model, adding panel-panel connectors decreases the seismic intensity at which failure occurs, even if it is responsible for a yielding in uplift of the bottom connectors at a lower seismic intensity. • As expected, increasing the number of panel-slab connectors reduces significantly the shear demand on each bottom connector, resulting in a sliding failure at a much higher intensity. In conclusion, the connection configuration can be improved in two different ways: either adding more panel-slab connectors, or removing the panel-panel connectors. However, the second option is likely to be incorrect if the foundation damping is not as high as assumed in the model. 121 5.2.3 Coefficient of friction The Concrete Design Handbook suggests using a coefficient of friction of 0.5 between the slab and the panels. However, panel motion generates some dust that can lower this friction coefficient. In order to understand how this phenomenon could affect the system behaviour, damage on the designed building was studied using different coefficients of friction µ. Figure 5.9 shows the sliding and rocking displacement demand for each five values of the friction coefficient µ. Figure 5.9 Sliding and rocking demand for various values of the coefficient of friction µ The results show that the coefficient of friction does not affect significantly the rocking response. As expected, as the friction coefficient increases, the system has larger capacity against sliding. In the case where friction coefficient equals 0.5 (the bench mark model), the sliding failure 122 occurred at Sa(Tr) = 1.05 g. When the coefficient of friction is reduced to 0.3 and 0, the shear failure occurs at Sa(Tr) = 0.90g and 0.72g respectively. This is a 31% reduction in the intensity at which failure occurs. It can be noted that friction between the slab and the panels accounts for only 19% of the static sliding resistance of the building, whereas it accounts for 31% of the dynamic sliding resistance. This is probably due to the fact that friction has a very high potential to dissipate energy due to its hysteretic behaviour. 5.2.4 Roof diaphragm stiffness and strength This part studies how the roof diaphragm characteristics can affect the response of the system. Two roof parameters were studied: its stiffness and its strength. Stiffness values of 70%, 100%, 150%, 200%, 300% and 1000% of the reference stiffness were studied, with unlimited strength. Corresponding roof periods are 1.24s, 1.04s, 0.85s, 0.73s, 0.60s and 0.32s. Strength of 1000kN, 2000kN, and 3000kN were studied, using the reference stiffness. To get the building response with a roof whose strength is limited, a bilinear elastic plastic material was affected to the roof spring, with a post-yielding stiffness ratio of 0.01. Results of the numerical analysis are shown on Figure 5.10. Behaviour using a strength of 3000kN is the same as the behaviour of the reference building, therefore it was not shown on the figure. 123 Figure 5.10 Demand parameters for various roof diaphragm characteristics Based on the results presented in Figure 5.10, it can be concluded that: • The roof stiffness has a huge impact on all the demand parameters. Maximum damage is observed when the stiffness is 300%, corresponding to a roof period of 0.60s. With this stiffness, shear failure occurs at a spectral acceleration 2.1 times lower than the failure spectral acceleration for the reference configuration. • Damages are not highest when the diaphragm is infinitely stiff. A roof period of 0.60s results in the higher damage than with a period of 0.32s • In all the configurations tried, sliding failure always occurs before rocking failure. 124 • As expected, if a stiffer roof diaphragm is used, the roof will experience a higher shear demand and a lower displacement. However, a rigid diaphragm will experience a lower shear demand than a diaphragm having a period of 0.6s because of resonance in the latter. • When the strength of the roof diaphragm is limited, sliding and overturning demands are considerably reduced. As a result, if the roof diaphragm is not strong enough, the connections will not experience significant damage. • Whereas those eight roof diaphragms result in very different sliding and rocking demand, they are all treated in the same way by the current design static procedure by NBCC2010. 125 6 - SUMMARY AND CONCLUSIONS 6.1 Summary The principal problem currently encountered in modelling the seismic response of tilt-up buildings is the complex behaviour of the EM5 connector. Previous experimental studies were performed on this connector, focusing on the coupling between the axial and the shear response. However, those tests had not been used to develop a cyclic numerical model of the EM5 connection. This research uses a collection of experimental data from previous studies to develop a numerical model of the EM5 connection. Another numerical model aimed at studying the nonlinear dynamic behaviour of tilt-up buildings with solid panels was then implemented. The building model was calibrated to a typical tilt-up building, and the coupled EM5 model was incorporated. This model was studied using incremental nonlinear dynamic analyses in order to assess the seismic behaviour and performance of a typical building. A parametric study was conducted to compare the design configuration with alternative connector configurations. In addition the effects of changing roof parameters (strength and stiffness), friction coefficient at the panel-slab interface, and EM5 envelope parameters was studied. 126 6.2 Discussion of results A numerical model of the EM5 connector was developed based on previous experimental studies. It takes into account the coupling between the shear and the axial behaviour. The model response showed a good agreement with the experimental test data. A two dimensional tilt-up building model was then developed. It assumes that panels are rigid and that the sliding displacement and rocking angle are identical for each panel. Contact and friction between the foundation slab and the panels were modeled. A rigid contact condition was used. The static and dynamic building response was validated against hand calculation and two other tilt-up models. Analysis of the building response showed that this contact model resulted in a high energy dissipation every time the panels pound the ground, thus lowering the rocking demand during seismic loading. This finding highlighted that contact condition between the foundation slab and the panels plays a significant role in the rocking demand on the panels. The following conclusions were made regarding the seismic performance of the prototype building: • The building is found to fail in sliding mechanism, whereas static analysis suggests it should fail in an overturning mechanism. • Based on the results of the Incremental Dynamic Analyses, the connection configuration designed following the recommendations of Weiler (2006) in the Concrete Design Handbook is conservative. When the ten selected ground motion records are scaled to a probability of exceedance of 2% in 50 years, connection yielding is observed for only one record. The high performance observed in the model is due to the following reasons: 1) the flexibility of the roof diaphragm is not taken into account in the design procedure; 2) the shear capacity of the panel-slab connector in the model is much higher than the one used in design. 3) Some of the earthquake energy was dissipated through rocking which was not accounted in the conventional equivalent force based approach. 127 • The combination of the EM5-EM3 does not allow to use the all the shear ductility capacity of the EM5 connector because of the low strength of the EM3 connector. When the tilt-up building is equipped with EM5-EM4 connectors, the building is likely to withstand twice the shaking intensity as the same building with the EM5-EM3 connector. Here are the conclusions of the parametric study: • Adding panel-panel connectors changes the vibration mode from rocking to sliding. This increases the stiffness of the tilt-up building, which results to higher shear demand on the roof and the panels. • Increasing the number of panel-slab connectors is the best way to improve the seismic resistance of the building, even though static analysis predicts an overturning final failure mechanism. • The horizontal location affects the rocking response of the tilt-up panels. Placing the connectors closer to the middle of the panel increases the cyclic ductility of the system. • Roof strength plays an important role in the damage of the tilt-up building. When the roof strength is low, it is likely that the damage will concentrate in the roof diaphragm, instead of the panel-to-panel or panel-to-slab connections. • Friction coefficient plays an important role in the sliding demand on the building. When friction coefficient is increase from 0 to 0.5, the tilt-up building is able to resist a shaking intensity which is 44% greater, whereas the static strength is increased by only 23%. This phenomenon is probably due to the high energy dissipation due to friction. • The sliding and rocking demand of the tilt-up building is greatly affected by the roof stiffness. When the period of the roof is increased from 0.6s to 1.04s by decreasing the roof stiffness, the seismic intensity at failure is increased by 118%. The roof stiffness is not taken into account in the equivalent static force procedure by 2010 NBCC: roof diaphragm is assumed to be rigid. This hypothesis results in a very conservative connection design for the prototype building. 128 6.3 Conclusion An empirical EM5 connection model that takes the axial and shear coupling into account was successfully implemented. A 2D finite-element software was developed in this study to analyse the nonlinear response of tilt-up buildings and incorporate the EM5 connection model. The software includes a detailed friction and contact algorithm. The software offers the ability to simulate the nonlinear static and dynamic response of the tilt-up buildings. It was found that foundation damping is a key part of the modeling of tilt-up structures with rigid wall panels. The seismic overturning demand can be greatly over or underestimated depending on the damping condition. In this study, the rocking response was probably underestimated because of the way contact is dealt with in the finite element software. However, the developed contact model would be realistic in the case of a highly damping foundation. The model was used to confirm the seismic performance of a prototype model whose connections are designed using the procedure described by Weiler (2006) in the Cement Association of Canada Concrete Design Handbook. The result shows the current design of the tilt-up building is conservative. A parametric study was performed using the model. It provided valuable information about how the connector configuration affects the seismic response of the structure, and it revealed the critical importance of the roof flexibility and strength in the building behaviour. 6.4 Recommendation for future research The EM5 numerical model developed in this study could be implemented in a 3D finite element software in order to study more complicated building geometries. All the mathematical equations that are needed to do this work are explained in Chapter 3. 129 Foundation damping was found to affect significantly the rocking demand of tilt-up buildings. As a result, research about the modeling of damping of tilt-up foundations would provide valuable information for calibrating the future tilt-up building models with solid panels. The result of such research might enable to use a higher force reduction factor for the rocking motion in order to take the damping due to pounding into account. The developed model is based on experimental test results on the connections. Validation of the system response could be provided by experimental testing of entire panels. 130 REFERENCES Baker, J.W. and Cornell, C.A., (2006), “Spectral Shape, Epsilon and Record Selection”, Earthquake Engineering & Structural Dynamics, 34 (10) 1193-1217. Chopra, A.R., 2000, “Dynamics of Structures: Theory and Applications to Earthquake Engineering, Second Edition”, Prentice Hall. Canadian Standard Association (CSA), 2004, “Design of Concrete Structures”, CSA-A23.3-94, CSA, Canada. Devine, F. (2008) “Investigation of Concrete Tilt-up Wall to Base Slab Connections,” M.A.Sc.Thesis, Department of Civil Engineering, University of British Columbia, Vancouver, Canada. Essa, H.S., Tremblay, R., and Rogers, C.A. (2003), “Behaviour of Roof Deck Diaphragms under Quasi-Static Cyclic Loading,” Journal of Structural ngineering, 129(12), 1658-1666. Lemieux, K., Sexsmith, R., and Weiler, G. (1998) “Behaviour of Embedded Steel Connectors in Concrete Tilt-up Panels,” ACI Structural Journal, V. 95, No. 4, July-August, pp. 400- 411. 131 National Research Council of Canada (NRCC), 2010, “National Building Code of Canada (NBCC), Part 4”, NBCC 2010, NRCC, Canada. Olund, O. (2009) “Analytical study to investigate the seismic Performance of single story tilt- up Structures” M.A.Sc.Thesis, Department of Civil Engineering, University of British Columbia, Vancouver, Canada. PEER, 2006, PEER NGA Database, Pacific Earthquake Engineering Research Centre, University of California, Berkeley, California. Vamvatsikos, D. and Cornell, C.A. (2002) “Incremental Dynamic Analysis”, Earthquake Engineering and Structural Dynamics, 31: 491-514. Wallace, J (2006) “One story Tilt-up Concrete Building” in NBCC 2005 Seismic design seminar Weiler, G. "Chapter 13 - Tilt-up Concrete Wall Panels," Cement Association of Canada Concrete Design Handbook, 2006. 132 APPENDIX A. DESIGN NOTES ON TYPICAL BUILDING Connection design for in-plane loading of the short direction of the typical building References: 1. CAN/CSA A23.3 2. "Concrete Design Handbook - Third Edition" by the Cement Association of Canada Building geometry: • 2 rows of 4 panels in the short direction, 2 rows of 8 panels in the long direction • Panel dimensions = 9.144m x 7.62m x 0.184m Seismic acceleration for computation of the equivalent static force: • Building period following 2010 NBCC : shT na 26.005.0 4 3 == • S(0.2)=1.0 • Fa=1 (soil class C) • Ie=1.0 • Higher mode factor = 1.0 • Rd=1.0 Ro=1.0 (Will use different Reduction factors for sliding and rocking) ⇒ Seismic acceleration = E=2/3 * Fa*S*Ie*Mv/RdRo=0.67g Short direction Shear due to roof and out-of-plane panels on one in-plane panel row: • Number of in-plane panels per row = nip = 4 133 • Number of out-of-plane panels per row = nop = 8 • Concrete volumique weight = 24kN/m3 • Panel weight = Wp =303kN • Roof dead load = qroof =1.0kN/m2 • Snow load = qsnow = 1.6kN/m2 in Surrey • Shear force due to the roof : Vroof=((qroof + 0.25 qsnow)*Aroof + 0.5*Wp*nop)*e*0.5=1670kN • Shear due to in-plane panels weight: Vpanels=E*nip*Wp*0.5=806kN • Total base shear: Vf= 2476kN Sliding strength: Use Rd=1.0 and Ro=1.3 • Shear capacity of connectors: VrEM3=110kN VrEM5=125kN • Shear demand on connectors: VDs=Vf/(RdRo)-0.5*nip*µ*Wp=1300kN • Number of panel-slab connections needed per panel: Nps = VDs /(4* VrEM3)=2.95 ⇒ Need 3 panel-slab connectors per panel Overturning strength: Use Rd=1.5 and Ro=1.3 • Overturning moment demand: Mf=(Vroof*h+Vpanel*h/2)/(RdRo)=9720kN.m • Resisting moment due to panel weight : Mr=Wp*d/2*nip/2=4590kN.m • Required panel-panel connection strength per panel: Vfr=(Mf-Mr)/(d*(nip/2-1))=224kN • Shear resistance of one EM5-EM5 connection: VrEM5=125kN 134 • Number of panel-panel connectors needed per panel: Npp = Vfr / VrEM5 = 1.79 ⇒ Need 2 panel-panel connectors per panel Connection design : 4-2-3 135 APPENDIX B. TIME-HISTORY ANALYSIS RESULTS This appendix shows some of the time-history analysis results that were used in Figure 5.2 to perform Incremental Dynamic Analysis. The sliding, rocking and roof displacements are displayed. NGA 3437 Sa = 0.4g 136 Sa = 0.6g Sa = 0.8g 137 Sa = 1.0g NGA 727 FN Sa = 0.4g 138 Sa = 0.6g Sa = 0.8g 139 NGA 737FN Sa = 0.4g Sa = 0.6g 140 Sa = 0.8g Sa = 1.0g 141 APPENDIX C. MODEL CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Nonlinear simulation of tilt-up system % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % clean start clear all; close all; clc; clear; %% Geometry - units = kN & m % Panel Height : H = 9.144; % Panel width : d = 7.62; % Panel thickness : t = 0.184; % Panel weight: wp = H*d*t*24; % Friction coefficient: mu = 0.50; % Damping coef: ksi=0.03; % Number of panels : np = 4; % Number of connectors per panel (bottom-pp-0) : % Add columns to study multiple configurations nc1List = [ 3 ; 2 ; 0 ]; % Roof mass: mRoof=205; %Unit= metric ton ( m.a=f => ton * m * s^-2 = kN ) % Roof stiffness: kRoof0=7500; % Roof strength: kyRoof=50000; %Put high value for elastic roof % Newmark parameters g=0.5; %Gamma b=0.25; %Beta %Coupling (1) or not coupling (0) : degr=0; %% GM loading %GM number choice for analysis: %GMnumList=[1 2 3 4 5 6 7 8 9 10 11 12]; 142 GMnumList=[2]; %PGA in g units : PGAglist=[0.75]; %Period used to scale Sa --------------------------------------------- Tsa=1.04 %Sa list for ground motion scaling (unit = m/s2) : SaTrList0=[0.66;1;1.5;2;2.5;3;3.5;4;4.5;5;5.5;6;6.5;7;7.5;8;9;10;11;12;13]; % Convert PGA in m/s2 units : PGAlist = 9.81*PGAglist; scList = SaTrList0; %Linear interpolation between GM steps>> %Number of added time step at every time step ____ plus one : nSteps =2; %Max number of time step ntGMmax = 5000; ntsW=10; HX=H+0.01; %% Initialization of variables % Friction Coef / Roof stiffness / Roof strength % Add columns to analyse different configurations MuKrFy=[ mu ; kRoof0 ; kyRoof ] ; % Distance from middle to Contact points: dc = d/2 - 0.075; dtList=[1]; %Connection parameter list: % [F1 D1 F2 Wref] for bottom shear, uplift, and panel-panel shear F1D1F2Wr_hvs= [ 207 0.002 254 21.1 34 0.002 80 16.7 207 0.0026 254 21.1]; % failure displacemetn for bottom shear, uplift, and panel-panel shear df_hvs = [ 0.0165 0.1 0.033 ] ; % Damage storage variables: po = size(nc1List); p1 = size(MuKrFy); ndp = p1(1,2); NppF = zeros(length(scList),po(1,2),ndp,length(GMnumList)); NppR = zeros(length(scList),po(1,2),ndp,length(GMnumList)); NpsF = zeros(length(scList),po(1,2),ndp,length(GMnumList)); NpsR = zeros(length(scList),po(1,2),ndp,length(GMnumList)); MaxRoofDriftList = zeros(length(scList),po(1,2),ndp,length(GMnumList)); MaxRoofDDList = zeros(length(scList),po(1,2),ndp,length(GMnumList)); ResSlide = zeros(length(scList),po(1,2),ndp,length(GMnumList)); MaxVt_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxVr_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxVp_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxMt_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxMr_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxMp_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxSl_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) MaxRD_l = zeros(length(scList),po(1,2),ndp,length(GMnumList)) SaTrList = zeros(length(scList),length(GMnumList)); Sa02List = zeros(length(scList),length(GMnumList)); %% START ITERATE ON ALL THE GROUND MOTIONS ******************************** for Zgm = 1:length(GMnumList) 143 GMnum=GMnumList(1,Zgm) %% START ITERATE ON PARAMETERS ******************************************** for Zdp = 1:ndp 'New set of parameters ***** ' mu = MuKrFy(1,Zdp); kRoof = MuKrFy(2,Zdp); Fy = MuKrFy(3,Zdp); F1D1R2W_axial=F1D1F2Wr_hvs(Zdp,5:8) dfH = df_hvs(Zdp,1); dfV = df_hvs(Zdp,2); dfS = df_hvs(Zdp,3); DTM=dtList(1,Zdp); %% START ITERATE ON CONFIGURATIONS **************************************** for Zcf = 1:po(1,2) 'New Connection configuration ---- ' nc1=nc1List(1:3,Zcf) % Location of connectors : [ xc,xc1 ] = getConnLocation( nc1,d,H,np ); % Number of local DOF per panel: n1 = 2*(nc1(1)+2*nc1(2)); % Number of total local DOF : n = n1*np; %% START ITERATE ON ALL THE GM SCALINGS *********************************** for Zsc = 1:length(scList) clear K K0 M C XListAPPLIED_dyn XListAPPLIED xListC xList Param SaTrT=SaTrList0(Zsc,1) % Create matrix of EM5 parameters : [ Param ] = getParamMatrix( np,nc1,F1D1F2Wr_hvs(Zdp,:) ); RoofP = [0;0;kRoof;Fy;0.1*kRoof;0]; kax=1000; %% APPLIED FORCE %Triangular load: %[ lF1p ] = getLoadList_Triangular(Xmax,HX,H,nts,ncyc,wp,np,CoefRoof); %Trapezoidal load: %[ lF1p ] = getLoadList_Trapeze(Xmax,HX,H,nts,ncyc,wp,np,CoefRoof,ntsTRAPEZE); %Load and sudden drop - for pounding study: %[ lF1p ] = getLoadList_LoadAndDrop(Xmax,HX,H,nts,ncyc,wp,np,CoefRoof,ntsTRAPEZE); %Seismic loading scaled at PGA: %[ lF1p,dt,accList,accList1,PGA,PGAusedNow,PGAv,PGAbrut ] = getLoadList_GMpga(PGA,ntsW,wp,np,mRoof,nSteps,ntGMmax,GMnum); %Seismic loading scaled at Sa(Tsa): [ lF1p,dt,accList,accList1,PGA,PGAbrut ] = getLoadList_GMsr(SaTrT,Tsa,ntsW,wp,np,mRoof,nSteps,ntGMmax,GMnum); %Seismic loading scaled with unscaled ground motions: %[ lF1p,dt,accList,accList1,PGA ] = getLoadList_GM_noScaling(wp,np,mRoof,nSteps,ntGMmax,GMnum); %In case of seismic loading, get the spectrum 144 if Zcf==1 && Zdp==1 spectra = ResponseSpectraElastic(accList,dt,1,0.1,[0.01;0.2;0.5;1.15;2]); %spectra.acc Sa=spectra.totAcc; %spectra.psdAcc %<< bad because damping is high and Troof is high SaTrList(Zsc,Zgm)=Sa(4); Sa02List(Zsc,Zgm)=Sa(2); Sa001=Sa(1); Sa05=Sa(3); SaTr=Sa(4); spectra.pga; end dt=DTM*dt; h=dt; % XListAPPLIED = %[ total applied lateral force ; Total Moment about center ; panelweights ] XListAPPLIED = lF1p; %% Initialize lists %update nts for loading and unloading nts=length(lF1p); %2*nts-1; %initialize list of force that will be filled with external forces with %moment at corner (weight is included in moment): XList2DOF = zeros(3+np,nts); XListAPPLIED_dyn=zeros(3+np,nts); XdynList=zeros(3+np,nts); TrueXdyn=zeros(3+np,nts); %list of displacements, velocity, and acceleration (list of zeros that will be filled): xList = zeros(3+np,nts); xListC = zeros(3+np,nts); vList = zeros(3+np,nts); aList = zeros(3+np,nts); xEQ=zeros(3+np,nts); Xv=zeros(3+np,nts); TiltList = zeros(1,nts); %Initialize Sum of moments at bottom middle point PrevSum=0; %Initialize Sum of moments at bottom middle point SumM=0; %Inililize list of Rock/slide RS=zeros(4,nts); FrSignList=zeros(1,nts); TrialDispListC=zeros(2*(3+np),nts); TrialDispListFR=zeros(2*(3+np),nts); %Precision list: PrecList=zeros(3+np,nts); List1=zeros(8,nts); List2=zeros(8,nts); Mcon=0; ll4=zeros(4,nts); Fl2=zeros(3+np,nts); ListD=zeros(1,nts); %Listk=zeros(n,n,nts); Fconn=zeros(3+np,1); %% Recorders RecBotVERTd = zeros(np*nc1(1),nts); RecBotVERTf = zeros(np*nc1(1),nts); RecBotHORd = zeros(np*nc1(1),nts); RecBotHORf = zeros(np*nc1(1),nts); 145 RecPpSHEARd = zeros(nc1(2),nts); RecPpSHEARf = zeros(nc1(2),nts); RecFW = zeros(3*nc1(1),nts); RecRoof = zeros(3,nts); % Damage/failure parameters MaxRoofDrift = 0; MaxBaseShear = 0; %<< not relevant MaxV=0; MaxTopDisp=0; MaxVt = 0; MaxVr = 0; MaxVp = 0; MaxMt = 0; MaxMr = 0; MaxMp = 0; MaxSl = 0; MaxRD = 0; %% get the mass matix: [ M ] = getMassMatrixROOF( np,d,H,wp,mRoof ); % Damping matrix: [ K0,K0p,~,~,~,~ ] = getStiffnessMatrixROOF( - 1,n,np,nc1,n1,xc,xc1,d,H,kax,RoofP,0,Param,zeros(1,n),2,dfH,dfV,dfS ); w1=2*pi/1.04; w2=2*pi/0.5; C = 2*ksi/(w1+w2)*(2*w1*w2*M+1*K0); %1st AND 2nd mode damping %% ******* Start time iterations ****************************************** save('Yoi') for t=2:nts '***************** New time step *****************************'; t; XList2DOF(:,t)=zeros(3+np,1); xList(:,t)=zeros(3+np,1); xListC(:,t)=zeros(3+np,1); %Resest list of displacement directions: ways=zeros(1,n); wayR=0; %% What stiffness matrix should we use ? panel tilted // or \\ ?? XListAPPLIED_dyn(:,t)=XListAPPLIED(:,t)+1/b*M*(1/h*vList(:,t-1)+0.5*aList(:,t- 1))+C*(g/b*vList(:,t-1)+(0.5*g/b-1)*h*aList(:,t-1)); %[ Tilt,TiltList,PrevSum ] = getTiltDyn( XListAPPLIED_dyn,xListC,TiltList,HX,H,kRoof,np,PrevSum,t ) ; [ Tilt,TiltList,PrevSum ] = getTiltDyn2( XListAPPLIED_dyn,XListAPPLIED(:,t),xListC,TiltList,HX,H,kRoof,np,PrevSum,t ) ; %% Conversion matrices % Convert forcese and displacement at corner : R = eye(3+np); R(1,2)=-H/2; for j=3:2+np R(j,2)=Tilt*dc; end T=transpose(R); for Round=1:2 %% GET THE GLOBAL STIFFNESS MATRIX for xList(:,t-1)!!!!!! if t>=2 [ K,KnoRoof,B,k,A,S ] = getStiffnessMatrixROOF( Tilt,n,np,nc1,n1,xc,xc1,d,H,kax,RoofP,wayR,Param,ways,t,dfH,dfV,dfS ); 146 end %% CALCULATE THE NEW DISPLACEMENT : %%% Dynamic algo with contact-friction %%% XListAPPLIED_dyn=XListAPPLIED; XListAPPLIED_dyn(:,t)=XListAPPLIED(:,t)+1/b*M*(1/h*vList(:,t-1)+0.5*aList(:,t- 1))+C*(g/b*vList(:,t-1)+(0.5*g/b-1)*h*aList(:,t-1)); % Applied force is at middle XdynList(:,t)=T*XListAPPLIED_dyn(:,t); %<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<DEC03-Avoid instability for i=1:np XListAPPLIED_dyn(2+i,t)=min(XListAPPLIED_dyn(2+i,t),0); end %<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<DEC03-Avoid instability pseudoK = K+T*(1/(b*h^2)*M+g/(b*h)*C)*R; %<<< at the corner [ XList2DOFdyn,xListC,TrialDispListC,TrialDispListFR,RS,L,C1,C2,dX_FR2,dXFR,Cd ] = ContactAlgorithm4D( XList2DOF,XListAPPLIED_dyn,xListC,pseudoK,t,Tilt,TiltList,Round,np,H,d,mu,wp,kRoof,T rialDispListC,TrialDispListFR,RS); %XList2DOF(:,t)=XList2DOFdyn(:,t)-T*(1/b*M*(1/h*vList(:,t-1)+0.5*aList(:,t- 1))+C*(g/b*vList(:,t-1)+(0.5*g/b-1)*h*aList(:,t-1))); XList2DOFdynNOTDYN=XList2DOFdyn(:,t)-T*(1/b*M*(1/h*vList(:,t-1)+0.5*aList(:,t- 1))+C*(g/b*vList(:,t-1)+(0.5*g/b-1)*h*aList(:,t-1))); for j=1:3+np if Cd(j,1)==0 XList2DOF(j,t)=XList2DOFdyn(j,t) ; else ROW = zeros(1,3+np); ROW(1,j) = 1; XList2DOF(j,t)=ROW*XList2DOFdynNOTDYN; end end XList2DOF(:,t-1)=XList2DOFdyn(:,t-1); %just to deal with Tilt change if Round == 1 && n>0 xCur=A*xListC(1:2+np,t); xPrev=A*xListC(1:2+np,t-1); %Calculate and save the displacement directions for j=0:np-1 for i=1:nc1(1)*2 ways(1,i+j*n1)=sign(xCur(i+j*n1,1)-xPrev(i+j*n1,1)); end for i=nc1(1)*2+1:(nc1(1)+nc1(2))*2 % left side if j==0 else ways(1,i+j*n1)=sign((xCur(i+j*n1,1)-xCur(i+2*nc1(2)+(j-1)*n1,1))- (xPrev(i+2*nc1(2)+(j-1)*n1,1)-xPrev(i+j*n1,1))); end end for i=(nc1(1)+nc1(2))*2+1:(nc1(1)+2*nc1(2))*2 % right side if j==np-1 else ways(1,i+j*n1)=sign((xCur(i+j*n1,1)-xCur(i-2*nc1(2)+(j+1)*n1,1))- (xPrev(i-2*nc1(2)+(j+1)*n1,1)-xPrev(i+j*n1,1))); end end end 147 wayR=sign(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t)-(xListC(3+np,t-1)-xListC(1,t- 1)+H*xListC(2,t-1))); RecRoof(3,t)=wayR; ways; %Reset Xlist and xListC XList2DOF(:,t)=zeros(3+np,1); xList(:,t)=zeros(3+np,1); xListC(:,t)=zeros(3+np,1); end if wayR>-1000 else break end end %of Rounds 1 and 2 % if wayR>-1000 % else % break % end xList(:,t)=R*xListC(:,t); % Newmark equations: aList(:,t)=(xList(:,t)-xList(:,t-1)-h*vList(:,t-1))/(b*h^2)-(0.5/b-1)*aList(:,t-1); vList(:,t)=(g/(b*h))*(xList(:,t)-xList(:,t-1))+(1-g/b)*vList(:,t-1)+(1- 0.5*g/b)*h*aList(:,t-1); if t>=159 %xListC(2,t)*H>0.04 end %<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<DEC06-Avoid instability if aList(1,t)*(XListAPPLIED(1,t)+kRoof*(xListC(3+np,t-1)-xListC(1,t-1)+H*xListC(2,t- 1))-Fconn(1,1))<0 aList(1,t)=0;%aList(1,t)/2;%0% end %<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<DEC06-Avoid instability %Newest contact law !! if xList(2,t-1)~=0 && xList(2,t)==0 vList(2:2+np,t)=zeros(1+np,1); %DOnt change Vx and Ax now aList(2:2+np,t)=zeros(1+np,1); xListC(3:2+np,t)=zeros(np,1); xList(3:2+np,t) =zeros(np,1); end % If one component has the same value twice in a row, then a and v should % be zero : for i=1:2+np %roof does not need that since it's linear if xList(i,t)==xList(i,t-1) aList(i,t)=0; vList(i,t)=0; end end for i=3:2+np %roof does not need that since it's linear %if xListC(i,t)==0 && xListC(i,t-1)~=0 %<< new equation 148 if xListC(i,t)==0 %&& xListC(i,t-1)~=0 %<< new equation aList(i,t)=d/2*Tilt*aList(2,t)-H*Tilt*vList(2,t);%0; vList(i,t)=d/2*Tilt*vList(2,t);%0; end end if n>0 %No equilibrium check if there are no connectors. %% Equilibrium check [ Fconn,~,Fl ] = GetConForces( Param,xListC(:,t),nc1,np,RoofP,H,A,B,S,t,degr,[207 0.002 254 21.1 34 0.002 80 16.7 207 0.0026 254 21.1],dfH,dfV,dfS ); %Internal and external forces: ForceDueToConnectors_ExternalForcesAndDynamicOnes_AppliedForce=[Fconn XList2DOF(:,t)-T*(M*aList(:,t)+C*vList(:,t)) XList2DOF(:,t)]; ConnDisp=A*xListC(1:2+np,t); %Force correction: XList2DOF(:,t) = -Fconn+T*(M*aList(:,t)+C*vList(:,t)); %X2dof is at corner, acc forces are at middle PrecList(:,t)=(Fconn+XList2DOF(:,t)-M*aList(:,t)- C*vList(:,t))./(abs(Fconn)+0.000001*ones(3+np,1)); %% Update connection and roof parameters: [ Param,RecBotVERTd,RecBotVERTf,RecBotHORd,RecBotHORf,RecPpSHEARd,RecPpSHEARf,RecFW ] = updateParam( Param,xListC(:,t),nc1,np,A,t,k,RecBotVERTd,RecBotVERTf,RecBotHORd,RecBotHORf,RecPpSH EARd,RecPpSHEARf,degr,RecFW,[207 0.002 254 21.1 34 0.002 80 16.7 207 0.0026 254 21.1],dfH,dfV,dfS ); [ RoofP,RecRoof ] = updateRoofP( RoofP,xListC(:,t),RecRoof,np,H,t); %% Stop iterations if failure is reached: if KnoRoof(1,1)==0 || KnoRoof(2,2)==0 'Failure reached' break end %% Save system state: List1(1,t)=t; TrueXdyn(:,t) = XListAPPLIED(:,t)-M*aList(:,t)-C*vList(:,t); xEQ(:,t)=(T\K0/R)\XListAPPLIED(:,t); %<<< at middle Xv(:,t)=C*vList(:,t); else [ RoofP,RecRoof ] = updateRoofP( RoofP,xListC(:,t),RecRoof,np,H,t); end if abs(xListC(3+np,t))>MaxRoofDrift MaxRoofDrift=abs(xListC(3+np,t)); MaxRoofDD=abs(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t)); end if abs(xListC(1,t)-H*xListC(2,t))>MaxTopDisp MaxTopDisp=abs(xListC(3+np,t)); end L1=zeros(1,3+np); L2=zeros(1,3+np); L1(1,1)=1; L2(1,2)=1; for i=np L2(1,2+i)=Tilt*d/2; end 149 %Vt=abs(XListAPPLIED(1,t)-L1*(M*aList(:,t)+C*vList(:,t))+kRoof*(xListC(3+np,t)- xListC(1,t)+H*xListC(2,t))); if Fconn(1,1)==0 Vt=0; else Vt=abs(Fconn(1,1))+mu*np*wp; end if Vt>MaxVt MaxVt=Vt; end Vp=abs(XListAPPLIED(1,t)-L1*(M*aList(:,t)+C*vList(:,t))); %Vp=abs(XListAPPLIED(1,t)); if Vp>MaxVp MaxVp=Vp; end %Vr=abs(kRoof*(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t))); Vr=abs(kRoof*(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t))); if Vr>MaxVr MaxVr=Vr; end %Mt=abs(-H/2*(XListAPPLIED(1,t)-L1*(M*aList(:,t)+C*vList(:,t)))- H*kRoof*(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t))-L2*(M*aList(:,t)+C*vList(:,t))); Mt=abs(Fconn(2,1))+np*wp*d/2; if Mt>MaxMt MaxMt=Mt; end %Mr=abs(H*kRoof*(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t))); Mr=abs(H*kRoof*(xListC(3+np,t)-xListC(1,t)+H*xListC(2,t))); if Mr>MaxMr MaxMr=Mr; end %Mp=Mt-Mr; Mp=abs(H/2*(XListAPPLIED(1,t)-L1*(M*aList(:,t)+C*vList(:,t)))); %Mp=abs(XListAPPLIED(1,t))*H/2; if Mp>MaxMp MaxMp=Mp; end %Update max sliding and max drift: if abs(xListC(1,t))>MaxSl MaxSl=abs(xListC(1,t)); end if abs(xListC(2,t))*H>MaxRD MaxRD=abs(xListC(2,t))*H; end end %% Save damage state and failure [ NppF(Zsc,Zcf,Zdp,Zgm),NppR(Zsc,Zcf,Zdp,Zgm),NpsF(Zsc,Zcf,Zdp,Zgm),NpsR(Zsc,Zcf,Zdp,Z gm) ] = getDamageState( Param,np,nc1 ); MaxRoofDriftList(Zsc,Zcf,Zdp,Zgm)=MaxRoofDrift; MaxRoofDDList(Zsc,Zcf,Zdp,Zgm)=MaxRoofDD; ResSlide(Zsc,Zcf,Zdp,Zgm)=xListC(1,t); %xCmegaMatrix(:,:,Zsc,Zcf) = xListC; MaxVt_l(Zsc,Zcf,Zdp,Zgm)=MaxVt; MaxVr_l(Zsc,Zcf,Zdp,Zgm)=MaxVr; MaxVp_l(Zsc,Zcf,Zdp,Zgm)=MaxVp; MaxMt_l(Zsc,Zcf,Zdp,Zgm)=MaxMt; MaxMr_l(Zsc,Zcf,Zdp,Zgm)=MaxMr; 150 MaxMp_l(Zsc,Zcf,Zdp,Zgm)=MaxMp; MaxSl_l(Zsc,Zcf,Zdp,Zgm)=MaxSl; MaxRD_l(Zsc,Zcf,Zdp,Zgm)=MaxRD; MaxRoofDD MaxRoofDrift psf=NpsF(Zsc,Zcf,Zdp,Zgm) psr=NpsR(Zsc,Zcf,Zdp,Zgm) %Plot results for this time-history analysis %Total applied horizontal force: TotalLat = XListAPPLIED(1,:) ; %Time list ListT=List1(1,:)*dt; %all connectors behaviour figure; for j=1:np for i=1:nc1(1) subplot(5,8,8*(i-1)+2*(j-1)+1);plot(RecBotHORd(i+(j- 1)*nc1(1),:),RecBotHORf(i+(j-1)*nc1(1),:)); grid subplot(5,8,8*(i-1)+2*(j-1)+2);plot(RecBotVERTd(i+(j- 1)*nc1(1),:),RecBotVERTf(i+(j-1)*nc1(1),:)); grid end if j<np && nc1(2)>1 subplot(5,4,16+j);plot(RecPpSHEARd(j,:),RecPpSHEARf(j,:)); grid end end % subplot(5,4,20);plot(RecRoof(1,:),RecRoof(2,:)); % xlabel('Roof force(kN) VS disp(m)') % grid %Applied force : figure; subplot(2,1,1); plot(ListT,XListAPPLIED(1,:),ListT,XListAPPLIED(3+np,:)); xlabel('Applied Force (kN) (blue=panel, green=roof) VS Time (s)') grid subplot(2,1,2); % All the displacemnts for a given time: %plot(ListT,xListC(1,:),ListT,- H*xListC(2,:),ListT,xListC(3,:),ListT,xListC(4,:),ListT,xListC(5,:),ListT,xListC(6,: ),ListT,xListC(7,:)); plot(List1(1,:),xListC(1,:),List1(1,:),- H*xListC(2,:),List1(1,:),xListC(3,:),List1(1,:),xListC(4,:),List1(1,:),xListC(5,:),L ist1(1,:),xListC(6,:),List1(1,:),xListC(7,:)); xlabel('Sliding(blue) ; -Theta*H(green) ; Uplft panel 1(red), 2(cyan), 3(purple), 4(gold) ; Translation of roof () VS Time (s)') ylabel('(m)') grid legend('Sliding','-Theta*H','Uplft panel 1','Uplft panel 2','Uplft panel 3','Uplft panel 4','Translation of roof') % end end end save('Save_intermediate') end %% Save results in file ! % Save the whole environment : 151 save('Results') %% Print some results '************ Max disp & damages ************' scLiist = transpose(scList) for Zgm=1:length(GMnumList) Zgm for Zdp=1:ndp Zdp MuKrFy(:,Zdp) %ppF_ppR_psF_psR_config1=[ NppF(:,1,Zdp,Zgm),NppR(:,1,Zdp,Zgm),NpsF(:,1,Zdp,Zgm),NpsR(:,1,Zdp,Zgm) ] end end MaxRoofDriftList MaxRoofDDList ResSlide %% Multiple analysis plots % nn=1; % figure; % for Zcf=1:nn % MATpsF=[]; % for Zgm=1:length(GMnumList) % MATpsF=[MATpsF (1-0.002*Zgm)*NpsF(:,Zcf,1,Zgm)]; % end % subplot(nn,1,Zcf); % plot(PGAlist,MATpsF ) % xlabel('PGA (m.s-2)') % ylabel('Number of failed PS') % grid; % end % % MATroofD=[] % for Zgm=1:length(GMnumList) % MATroofD=[MATroofD MaxRoofDriftList(:,1,1,Zgm)]; % end % figure; % plot(PGAlist,MATroofD) % xlabel('PGA (m.s-2)') % ylabel('Roof drift in config1') % grid; % % MATroofD=[] % for Zgm=1:length(GMnumList) % MATroofD=[MATroofD MaxRoofDriftList(:,2,1,Zgm)]; % end % figure; % plot(PGAlist,MATroofD) % xlabel('PGA (m.s-2)') % ylabel('Roof drift in config2') % grid; % % % Variable = Sa(Troof) % % gmm=11; % SaL=[0.2:0.1:2.5]; % NppFsa=zeros(gmm,length(SaL)); 152 % for gm=1:gmm % NppFsa(gm,:)=ConvConn(SaTrList(:,gm),NpsF(:,1,1,gm),SaL); % % end % % % nn=4; % figure; % for Zcf=1:nn % MATpsF=[]; % for Zgm=1:length(GMnumList) % MATpsF=[MATpsF (1-0.001*Zgm)*NpsF(:,Zcf,1,Zgm)]; % end % subplot(nn,1,Zcf); % for Zgm=1:length(GMnumList) % plot(SaTrList(:,Zgm),MATpsF(:,Zgm)) % hold on % xlabel('Sa(T=Troof) (m.s-2)') % ylabel('Number of failed PS connectors') % grid; % end % end % title('Number of failed PS connectors for some configurations') % % MATroofD=[] % for Zgm=1:length(GMnumList) % MATroofD=[MATroofD MaxRoofDriftList(:,2,1,Zgm)]; % end % figure; % plot(SaTrList(:,1),MATroofD(:,1),SaTrList(:,2),MATroofD(:,2),SaTrList(:,3),MATroofD( :,3),SaTrList(:,4),MATroofD(:,4),SaTrList(:,5),MATroofD(:,5)) % xlabel('Sa(Troof) (m.s-2)') % ylabel('Roof drift in config2') % grid; % % % % Variable = Sa(T=0.2) % figure; % for Zcf=1:nn % MATpsF=[]; % for Zgm=1:length(GMnumList) % MATpsF=[MATpsF (1-0.001*Zgm)*NpsF(:,Zcf,1,Zgm)]; % end % subplot(nn,1,Zcf); % for Zgm=1:length(GMnumList) % plot(Sa02List(:,Zgm),MATpsF(:,Zgm)) % hold on % xlabel('Sa(T=0.2) (m.s-2)') % ylabel('Number of failed PS connectors') % %legend('GM1','GM2','GM3','GM4','GM5'); % grid; % end % end % title('Number of failed PS connectors for some configurations') % % MATroofD=[] % for Zgm=1:length(GMnumList) % MATroofD=[MATroofD MaxRoofDriftList(:,2,1,Zgm)]; % end % figure; % plot(Sa02List(:,1),MATroofD(:,1),Sa02List(:,2),MATroofD(:,2),Sa02List(:,3),MATroofD( :,3),Sa02List(:,4),MATroofD(:,4),Sa02List(:,5),MATroofD(:,5)) 153 % xlabel('Sa(T=0.2) (m.s-2)') % ylabel('Roof drift in config2') % grid; %% Single analysis Plots %Total applied horizontal force: TotalLat = XListAPPLIED(1,:) ; %Time list ListT=List1(1,:)*dt; figure; plot(List1(1,:),xListC(1,:),List1(1,:),- H*xListC(2,:),List1(1,:),xListC(3,:),List1(1,:),xListC(4,:),List1(1,:),xListC(5,:),L ist1(1,:),xListC(6,:),List1(1,:),xListC(7,:)); xlabel('Sliding(blue) ; -Theta*H(green) ; Uplft panel 1(red), 2(cyan), 3(purple), 4(gold) ; Translation of roof () VS Time (s)') ylabel('(m)') grid legend('Sliding','-Theta*H','Uplft panel 1','Uplft panel 2','Uplft panel 3','Uplft panel 4','Translation of roof') 154 function [ Kr,K,B,k,A,S ] = getStiffnessMatrixROOF( Tilt,n,np,nc1,n1,xc,xc1,d,H,kax,RoofP,wayR,Param,ways,t,dfH,dfV,dfS ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Calculate the stiffness matrix % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CALCULATION OF STIFFNESS MATRIX !! % % X x % M Theta % Y1 = K. y1 % Y2 y2 % Y3 y3 % % Moment at corner and (x,y) of corner dc=d/2-0.075; if n1==0 K=zeros(2+np,2+np); B=0; A=0; k=0; S=0; else %% A MATRIX : panel displacement -> node displacement ---------------------- % Number of DOF considering 3DOF per panel: ndof = 2+np; % Build matrix for 1st panel : A1 = zeros(n1,2); j=1; %number of the panel for i=1:n1 if i-2*floor(i/2)==1 %x displacements are on odd lines A1(i,3*(j-1)+1)=1; A1(i,3*(j-1)+2)=-xc(i+1); %xc(i+1) is the y disp of the node end if i-2*floor(i/2)==0 %y displacements are on even lines %A1(i,3*(j-1)+3)=1; A1(i,3*(j-1)+2)=xc(i-1)+Tilt*dc; %xc(i-1) is the x disp of the node end end A1; % Assemble matrices for multiple panels : c0 = zeros(n1,1); c1 = zeros(n1,1); for i=1:n1 if i-2*floor(i/2)==0 c1(i,1)=1; end end A = []; for i=1:np A2=A1; for j=1:np; if j==i 155 A2 = [A2 c1]; else A2 = [A2 c0]; end end A = [A ; A2]; end A; %check %% k MATRIX : STIFFNESS MATRIX in node coordinates ----------------------- % Uncoupled stiffness matrix : k1 = zeros(n,n); for j=0:np-1 for i=j*n1+1:j*n1+2*nc1(1) % Bottom connectors %1st->calc stiffnesses in case of increased/decreased displacement : if i-2*floor(i/2)==1 %x direction [ kp,kn ] = getTangents_H( Param(:,i),t,dfH ); end if i-2*floor(i/2)==0 %y direction [ kp,kn ] = getTangents_V( Param(:,i),dfV ); % kp=17000; % kn=17000; % 'linear bottom connectors' end %2nd->depending on what way the connection goes, use the correct %combination of kp,kn : if ways(1,i)==0 k1(i,i)=(kp+kn)/2; end if ways(1,i)==1 k1(i,i)=kp; end if ways(1,i)==-1 k1(i,i)=kn; end end for i=j*n1+2*nc1(1)+1:j*n1+2*nc1(1)+2*nc1(2) % Side left connectors if j>0 % No pp connector on the edges if i-2*floor(i/2)==1 %x direction if j==0 k1(i,i)=0; %non-modeled panel can rock, thus there is no stiffness (simple, but ok) else k1(i,i)=kax; %could be zero - useless end end if i-2*floor(i/2)==0 %y direction [ kp,kn ] = getTangents_H( Param(:,i),t,dfS ); if ways(1,i)==0 k1(i,i)=(kp+kn)/2; end if ways(1,i)==1 k1(i,i)=kp; end if ways(1,i)==-1 k1(i,i)=kn; end %k1(i,i)=ksh;%---- end end 156 end for i=j*n1+2*nc1(1)+2*nc1(2)+1:j*n1+n1 % Side right connectors if j<np-1 % No pp connector on the edges if i-2*floor(i/2)==1 %x direction if j==np-1 k1(i,i)=0; %non-modeled panel can rock, thus there is no stiffness (simple, but ok) else k1(i,i)=kax; end end if i-2*floor(i/2)==0 %y direction [ kp,kn ] = getTangents_H( Param(:,i),t,dfS ); if ways(1,i)==0 k1(i,i)=(kp+kn)/2; end if ways(1,i)==1 k1(i,i)=kp; end if ways(1,i)==-1 k1(i,i)=kn; end %k1(i,i)=ksh;%---------- end end end % if j==0 % k0=k1; % else % k0=[k0 zeros(j*n,n) ; zeros(n,j*n) k1]; % end end % Add coupling terms (for side connectors): % ------------------------------------------------ % | Abottom | % | Aleft | | k1 kcpl 0 ...| % | Aright -Aright | = | kcpl k1 kcpl ...| % | Abott | | 0 kcpl k1 ...| % | -Aright Aleft | | 0 ... % | Aright | % ------------------------------------------------ % <---ntrans----------> for j=0:np-2 for i=j*n1+2*nc1(1)+2*nc1(2)+1:j*n1+n1 k1(i,i+2*nc1(2)+2*nc1(1))=-k1(i,i); k1(i+2*nc1(2)+2*nc1(1),i)=-k1(i,i); end end %Don't be so negative (k is the opposite of the stiffness matrix - it %translates the equilibrium): k=-k1; k; %% B MATRIX : Get global forces from node forces ------------------------- %Transformation matrix for panel 1: B1=zeros(2,n1); for i=1:n1 if i-2*floor(i/2)==1 %x direction = odd row of force vector B1(1,i)=-1; 157 B1(2,i)=xc1(i+1)-H/2; end if i-2*floor(i/2)==0 %y direction = even row of force vector l1(1,i)=-1; B1(2,i)=-xc1(i-1); end end l0=zeros(1,n1); B = []; for i=1:np B2=B1; for j=1:np; if j==i B2 = [B2 ; l1]; else B2 = [B2 ; l0]; end end B = [B B2]; end %% S matrix : from moment at center to moment at contact point %S = [1 0 0 ; -H/2 -d/2*(-Tilt) 1]; S=[]; for i=1:ndof if i==2 ll=[-H/2 1]; for j=1:np ll=[ll -dc*(-Tilt)]; end else ll=zeros(1,ndof); ll(1,i)=1; end S=[S;ll]; end S; %% Merge all the matrices : K=S*B*k*A; end %% Add roof dof [ kp,kn ] = getTangent_Roof( RoofP ); if wayR==0 kRoof=(kp+kn)/2; end if wayR==1 kRoof=kp; end if wayR==-1 kRoof=kn; end % Panels only: K1 = [K zeros(2+np,1) ; zeros(1,2+np) 0]; % Roof effect : K2 = zeros(3+np,3+np); K2(1,1)=kRoof; K2(1,2)=-H*kRoof; 158 K2(1,3+np)=-kRoof; K2(2,1)=-H*kRoof; K2(2,2)=H^2*kRoof; K2(2,3+np)=H*kRoof; K2(3+np,1)=-kRoof; K2(3+np,2)=H*kRoof; K2(3+np,3+np)=kRoof; Kr=K1+K2; end 159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Calculate mass matrix % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [ M ] = getMassMatrixROOF( np,d,H,wp,mRoof ) %% CALCULATION OF MASS MATRIX !! M=zeros(np+3,np+3); M(1,1)=np*wp/9.8; %Units of mass is Metric ton M(2,2)=np*wp/9.8*(H^2+d^2)/12; for i=3:2+np M(i,i)=wp/9.8; end M(3+np,3+np)=mRoof; end 160 function [ Param,RecBotVERTd,RecBotVERTf,RecBotHORd,RecBotHORf,RecPpSHEARd,RecPpSHEARf,RecFW ] = updateParam( Param,xList_t,nc1,np,A,t,k,RecBotVERTd,RecBotVERTf,RecBotHORd,RecBotHORf,RecPpSHEARd ,RecPpSHEARf,degr,RecFW,P,dfH,dfV,dfS ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Update connection hysteretic parameters % % and record force and displacement % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Get local displacements: xi=A*xList_t(1:2+np,1); n2=2*(nc1(1)+2*nc1(2)); for j=1:np for i=1:nc1(1) %% Bottom********************************************************************* % Bottom vertical connector ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); % get new displacement of the connector D=xi((j-1)*n2+2*(i-1)+2,1); % Get the new parameters : [~,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W]=EM5_ModelVERT(p(1,1),p(2,1),p(3,1),p(4 ,1),p(5,1),p(6,1),p(7,1),p(8,1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15, 1),p(16,1),p(17,1),p(18,1),p(19,1),p(20,1),p(21,1),p(22,1),D,dfV); Param(1,(j-1)*n2+2*(i-1)+2)=Dlast; Param(2,(j-1)*n2+2*(i-1)+2)=Flast; Param(3,(j-1)*n2+2*(i-1)+2)=Dmax; Param(4,(j-1)*n2+2*(i-1)+2)=Dmin; Param(15,(j-1)*n2+2*(i-1)+2)=rForceMp; Param(16,(j-1)*n2+2*(i-1)+2)=rForceMn; Param(22,(j-1)*n2+2*(i-1)+2)=W; % Record Vertical displ of connector : RecBotVERTd((j-1)*nc1(1)+i,t)=Dlast; % record force RecBotVERTf((j-1)*nc1(1)+i,t)=Flast; % Bottom horizontal connector ------------ %------------------------------------------------Effect of coupling nn=3; if degr==1 && t-nn*floor(t/nn)==0 %newWref=P(1,4)/21.1*subplus(21.1-123*Param(3,(j-1)*n2+2*(i-1)+2)); newWref=P(1,4)/21.1*max((21.1-123*Param(3,(j-1)*n2+2*(i-1)+2)),0); % F1(5) = 207-1075*Dmax_v; %newF1=P(1,1)/207*subplus(207-1075*Param(3,(j-1)*n2+2*(i-1)+2)); newF1=P(1,1)/207*max((207-1075*Param(3,(j-1)*n2+2*(i-1)+2)),0); % F2(7) = 254-1145*Dmax_v; %newF2=P(1,3)/254*subplus(254-1145*Param(3,(j-1)*n2+2*(i-1)+2)); newF2=P(1,3)/254*max((254-1145*Param(3,(j-1)*n2+2*(i-1)+2)),0); % Change kr1 (line 18), kr2(19), ku(20) (which are proportional to F1/D1) %equation used : newkr = oldkr*newF1/oldF1 Param(18,(j-1)*n2+2*(i-1)+1)=Param(18,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); Param(19,(j-1)*n2+2*(i-1)+1)=Param(19,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); Param(20,(j-1)*n2+2*(i-1)+1)=Param(20,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); 161 % Change F1,F2,Wref Param(14,(j-1)*n2+2*(i-1)+1)=newWref; Param(5,(j-1)*n2+2*(i-1)+1)=newF1; Param(7,(j-1)*n2+2*(i-1)+1)=newF2; % axial failure => shear failure !! if Dmax>=dfV Param(3,(j-1)*n2+2*(i-1)+1)=0.5; %<< we say that max axial displ is 0.5m, so failed behaviour will be used at next time step. end end %---------------------------------------------------end of Coupling %Now change all the other parameters : p=Param(:,(j-1)*n2+2*(i-1)+1); % Need to get D of this DOF from xi (calculated from xList) ! D=xi((j-1)*n2+2*(i-1)+1,1); if t==42 [p(1,1);p(2,1);p(3,1);p(4,1);p(5,1);p(6,1);p(7,1);p(8,1);p(9,1);p(10,1);p(11,1);p(12 ,1);p(13,1);p(14,1);p(15,1);p(16,1);p(17,1);p(18,1);p(19,1);p(20,1);p(21,1);p(22,1); D;t]; end [~,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1), p(5,1),p(6,1),p(7,1),p(8,1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p (16,1),p(17,1),p(18,1),p(19,1),p(20,1),p(21,1),p(22,1),D,dfH); % 'True' force used by the program (prog does not calc connection % forces, it just uses stiffness): if i==1 Rec(2*(i-1)+5,t)=Param(2,(j-1)*n2+2*(i-1)+1)+(Dlast-Param(1,(j- 1)*n2+2*(i-1)+1))*(-k(1,1)); %Flast =Param(2,(j-1)*n2+2*(i-1)+1)+(Dlast-Param(1,(j- 1)*n2+2*(i-1)+1))*(-k(1,1)); end Param(1,(j-1)*n2+2*(i-1)+1)=Dlast; Param(2,(j-1)*n2+2*(i-1)+1)=Flast; Param(3,(j-1)*n2+2*(i-1)+1)=Dmax; Param(4,(j-1)*n2+2*(i-1)+1)=Dmin; Param(15,(j-1)*n2+2*(i-1)+1)=rForceMp; Param(16,(j-1)*n2+2*(i-1)+1)=rForceMn; Param(22,(j-1)*n2+2*(i-1)+1)=W; % Record Horizontal force/displ of connector : RecBotHORd((j-1)*nc1(1)+i,t)=Dlast; RecBotHORf((j-1)*nc1(1)+i,t)=Flast; % Record new F1, F2 and W of bottom connectors 1 and 2 RecFW(3*i-3+3,t)=Param(14,(j-1)*n2+2*(i-1)+1); %Wref RecFW(3*i-3+1,t)=Param(5,(j-1)*n2+2*(i-1)+1); %F1 RecFW(3*i-3+2,t)=Param(7,(j-1)*n2+2*(i-1)+1); %F2 % Coupling : shear failure => axial failure !! if Dmax>=dfH Param(3,(j-1)*n2+2*(i-1)+2)=0.5; %<< we say that max axial displ is 0.5m, so failed behaviour will be used at next time step. end end for i=nc1(1)+1:nc1(1)+nc1(2) if j==1 162 %% left edge *********************************************************** % no need - look in older versions if needed else %% panel-panel left (change Param even if an elastic connector is used- >simpler)********** % vertical connector ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); D=-xi((j-1)*n2+2*(i-1)+2,1)+xi((j-2)*n2+2*(i+nc1(2)-1)+2,1); [~,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1), p(5,1),p(6,1),p(7,1),p(8,1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p (16,1),p(17,1),p(18,1),p(19,1),p(20,1),p(21,1),p(22,1),D,dfS); Param(1,(j-1)*n2+2*(i-1)+2)=Dlast; Param(2,(j-1)*n2+2*(i-1)+2)=Flast; Param(3,(j-1)*n2+2*(i-1)+2)=Dmax; Param(4,(j-1)*n2+2*(i-1)+2)=Dmin; Param(15,(j-1)*n2+2*(i-1)+2)=rForceMp; Param(16,(j-1)*n2+2*(i-1)+2)=rForceMn; Param(22,(j-1)*n2+2*(i-1)+2)=W; if j>=2 % Record shear force/displ of connector : RecPpSHEARd(j-1,t)=-RecPpSHEARd(j-1,t)+Dlast; RecPpSHEARf(j-1,t)=Flast; end end end for i=nc1(1)+nc1(2)+1:nc1(1)+2*nc1(2) if j==np %% right edge************************************************** else %% panel-panel right*************************************************** % vertical direction ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); D=-xi((j-1)*n2+2*(i-1)+2,1)+xi(j*n2+2*(i-nc1(2)-1)+2,1); [~,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1), p(5,1),p(6,1),p(7,1),p(8,1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p (16,1),p(17,1),p(18,1),p(19,1),p(20,1),p(21,1),p(22,1),D,dfS); Param(1,(j-1)*n2+2*(i-1)+2)=Dlast; Param(2,(j-1)*n2+2*(i-1)+2)=Flast; Param(3,(j-1)*n2+2*(i-1)+2)=Dmax; Param(4,(j-1)*n2+2*(i-1)+2)=Dmin; Param(15,(j-1)*n2+2*(i-1)+2)=rForceMp; Param(16,(j-1)*n2+2*(i-1)+2)=rForceMn; Param(22,(j-1)*n2+2*(i-1)+2)=W; if j<np % Record 2nd part of shear displ of connector : RecPpSHEARd(j,t)=Dlast; end end end end end 163 function [ Param ] = getParamMatrix( np,nc1,P ) % Builds the matrix containing the pinching parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get material parameters % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Bottom Horizontal behavior : zz = zeros(4,1); env = [ P(1,1); P(1,2); P(1,3); 0.015]; %F1,D1,F2,D2 rDF = [ 0.8; 0.8; 0.55; 0.55]; %rDispN,rDispP,rForceN,rForceP gam= 10; Wref= P(1,4); rfM= [0;0];%RForceMp,rForceMn bkr= [ 1; env(1)/env(2); 0.1*env(1)/env(2) ; env(1)/env(2); 0.05];%beta,kr1,kr2,ku,alpha ColBOT_H = [zz;env;rDF;gam;Wref;rfM;bkr;0]; % Bottom vertical behavior : env = [ P(1,5); P(1,6); P(1,7); 0.06]; %F1,D1,F2,D2 rDF = [0.8; 0.8; 0.5; 0.4]; %rDispN,rDispP,rForceN,rForceP gam= 10; Wref= P(1,8); rfM= [ 0; 0];%RForceMp,rForceMn bkr= [ 30; env(1)/env(2); 0.1*env(1)/env(2); env(1)/env(2); 0.05 ];%beta,kr1,kr2,ku,alpha ColBOT_V = [zz;env;rDF;gam;Wref;rfM;bkr;0]; % Panel-Panel horizontal behavior : % Elastic model for this connection ColPP_H = zeros(22,1); % Panel-Panel vertical behavior (half as stiff as horizontal bottom connectors) : env = [ P(1,9); P(1,10); P(1,11); 0.030]; %F1,D1,F2,D2 rDF = [ 0.8; 0.8; 0.55; 0.55]; %rDispN,rDispP,rForceN,rForceP gam= 10; Wref= P(1,12); rfM= [0;0];%RForceMp,rForceMn bkr= [ 1; env(1)/env(2); 0.1*env(1)/env(2) ; env(1)/env(2); 0.05];%beta,kr1,kr2,ku,alpha ColPP_V = [zz;env;rDF;gam;Wref;rfM;bkr;0]; % Side connectors behavior : if nc1(3)==0 % Case of no edge connector : ColSide_H = zeros(22,1); ColSide_V = zeros(22,1); else % With as many edge connectors as pp connectors : %ColSide_H = ColPP_H; ColSide_H = [zz;env;rDF;gam;Wref;rfM;bkr;0;7]; ColSide_V = ColPP_V; end Param = []; for i=1:np for j=1:nc1(1) Param=[Param ColBOT_H ColBOT_V]; end for k=1:nc1(2) if i==1 164 Param=[Param ColSide_H ColSide_V]; else Param=[Param ColPP_H ColPP_V]; end end for l=1:nc1(2) if i==np Param=[Param ColSide_H ColSide_V]; else Param=[Param ColPP_H ColPP_V]; end end end end function [ lF1p,dt,accList,accList1,PGA,PGAbrut ] = getLoadList_GMsr(SaTrT,Tr,ntsW,wp,np,mRoof,nSteps,ntMAX,GMnum) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get scaled seismic laoding % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %ntsW=floor(nts/10); if ntsW<=4 ntsW=5; end lh = zeros(1,ntsW-1); lv1 = [0:(-wp/(ntsW-2)):-wp]; lv=[]; for i=1:np lv = [lv ; lv1 ] ; end lm = zeros(1,ntsW-1); lr = zeros(1,ntsW-1); Lweight = [lh ; lm ; lv ; lr]; switch GMnum case 1 GMmat = xlsread('NGA_125FRIULI_FN'); GMmatV = xlsread('NGA_125FRIULI_vert'); case 2 GMmat = xlsread('NGA_125FRIULI_FP'); GMmatV = xlsread('NGA_125FRIULI_vert'); case 3 GMmat = xlsread('NGA_265VICT_FN'); GMmatV = xlsread('NGA_265VICT_vert'); case 4 GMmat = xlsread('NGA_727SUPERST_FN'); GMmatV = zeros(1500,5); case 5 GMmat = xlsread('NGA_739LOMAP_FN'); GMmatV = xlsread('NGA_739LOMAP_vert'); case 6 GMmat = xlsread('NGA_802LOMAP_FP'); GMmatV = xlsread('NGA_802LOMAP_vert'); case 7 GMmat = xlsread('NGA_952NORTHR_FP'); GMmatV = xlsread('NGA_952NORTHR_vert'); 165 case 8 GMmat = xlsread('NGA_1010NORTHR_FP'); GMmatV = xlsread('NGA_1010NORTHR_vert'); case 9 GMmat = xlsread('NGA_1111KOBE_FN'); GMmatV = xlsread('NGA_1111KOBE_vert'); case 10 GMmat = xlsread('NGA_2734CHICHI04_FP'); GMmatV = xlsread('NGA_2734CHICHI04_vert'); case 11 GMmat = xlsread('NGA_3473CHICHI06_FN'); GMmatV = xlsread('NGA_3473CHICHI06_vert'); case 12 GMmat = xlsread('NGA_3473CHICHI06_FP'); GMmatV = xlsread('NGA_3473CHICHI06_vert'); end dt=GMmat(1,2); nt=GMmat(1,1); nt=nt-(nt-5*floor(nt/5)); accList=zeros(1,nt); %% Horizontal motion PGA = 0; size(GMmat) for i=1:nt/5 accList(1,5*i-4:5*i)=GMmat(i+1,1:5); if abs(accList(1,i))>PGA PGA=abs(accList(1,i)); end if abs(accList(1,i+1))>PGA PGA=abs(accList(1,i+1)); end if abs(accList(1,i+2))>PGA PGA=abs(accList(1,i+2)); end if abs(accList(1,i+3))>PGA PGA=abs(accList(1,i+3)); end if abs(accList(1,i+4))>PGA PGA=abs(accList(1,i+4)); end end spectra = ResponseSpectraElastic(accList,dt,1,0.1,[Tr]); %spectra.acc Sa=spectra.totAcc; %spectra.psdAcc %<< bad because damping is high and Troof is high SaTr=Sa(1); PGAbrut=PGA; accList=accList*SaTrT/SaTr; accList=transpose(accList); accList1=transpose(accList); [ accList,dt ] = InterpolL( accList,dt,nSteps,ntMAX ); %% Vertical motion PGAv = 0; for i=1:nt/5 accListV(1,5*i-4:5*i)=GMmatV(i+1,1:5); if abs(accListV(1,i))>PGAv PGAv=accListV(1,i); 166 end if abs(accListV(1,i+1))>PGAv PGAv=accListV(1,i+1); end if abs(accListV(1,i+2))>PGAv PGAv=accListV(1,i+2); end if abs(accListV(1,i+3))>PGAv PGAv=accListV(1,i+3); end if abs(accListV(1,i+4))>PGAv PGAv=accListV(1,i+4); end end accListV=accListV*SaTrT/SaTr; % The same coefficient is used to scale both horizontal and vertical GMs accListV=transpose(accListV); [ accListV,dt0 ] = InterpolL( accListV,dt,nSteps,ntMAX ); %% Suite et fin... % size(accListV) % size(accList) % size(ones(np,nt )) nt = length(accList); lFeq = [-np*wp*accList/9.8 ; zeros(1,nt) ; -wp*(ones(np,1)*accListV/9.8+ones(np,nt )) ; -mRoof*accList ] ; lF1p = [Lweight lFeq]; %PGAvSC = PGAv*PGAt/PGA; end 167 function [F,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W] = EM5_Model(Dlast,Flast,Dmax,Dmin,F1,D1,F2,D2,rDispN,rDispP,rForceN,rForceP,gamma,Wref ,rForceMp,rForceMn,beta,kr1,kr2,ku,alpha,W,D,Df) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % EM5 shear behaviour % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% F=0; [Flast, EM5_MetaBackBone(F1,D1,F2,D2,Dlast), F1, F2]; Dlast; kr = kr2+(kr1-kr2)/(1+W/(alpha*Wref)); if Dmax<D1 && Dmin>-D1 kr=ku; end % initialisation if Dmax == 0 && Dmin == 0 && D ==0 F=0; end % update rDisp and rForce (note: the value of rDisp and rForce is note modified outside the procedure): rDispP = max(W/Wref*rDispP,0.01); rDispN = rDispP; rForceP = rForceP*(1-rForceMp); rForceN = rForceP*(1-rForceMn); %% If displacement increases ----------------------------- if D>=Dlast % + unloading with negative force if Flast<0 F = Flast+(D-Dlast)*ku; if F>0 %F=0; %<<<<<<<<<<<<<<<<<<<< Removed for improved system behaviour end else % + unloading with positive force - flat part (transition between unloading and reloading) if D<=rDispP*(Dmax-0.1*Dmin) && (Dmax>0 || Dmin<0) && Dmax>D1 && Dmin<-D1 F = Flast + (rForceP*EM5_BackBone(F1,D1,F2,D2,Dmax-0.1*Dmin)-Flast)*(D- Dlast)/(rDispP*(Dmax-0.1*Dmin)-Dlast); % The slope of the flat part should be lower than the slope of reloading : if (F-Flast)/(D-Dlast)>kr F = Flast + (D-Dlast)*kr; end %end % + reloading & + loading (this is the case where D is higher than displ of point "A" (flat to reload transition point)) %if D>rDispP*(Dmax-0.1*Dmin) else F = Flast + (D-Dlast)*kr; 168 if F > EM5_MetaBackBone(F1,D1,F2,D2,D) && D>D1 %&& (Dmax>D1||D>D1) && (Dmin<-D1||D<-D1) F = EM5_MetaBackBone(F1,D1,F2,D2,D); end if Dmax<D1 && Dmin>-D1 && F>kr*D F=F1/D1*D; end if Flast > EM5_MetaBackBone(F1,D1,F2,D2,Dlast) && Dlast>0 F = Flast - (D-Dlast)*kr*0.5;%/4; %we use -1/4 of the reloading slope... there might be some better choice %'Flast>MetaBB !!' end end end end % end of displacement increase------------------------- %% If displacement decreases ----------------------------- if D<Dlast % - unloading with positive force if Flast>0 F = Flast+(D-Dlast)*ku; if F<0 %F=0; %<<<<<<<<<<<<<<<<<<<< Removed for improved system behaviour end else % - unloading with negative force on flat part if D>=rDispN*(Dmin-0.1*Dmax) && Dmax>D1 && Dmin<-D1 %F = Flast + (rForceN*EM5_BackBone(F1,D1,F2,D2,Dmin-0.1*Dmax)-Flast)*(D- Dlast)/(rDispN*Dmin-Dlast); F = Flast + (rForceN*EM5_BackBone(F1,D1,F2,D2,Dmin-0.1*Dmax)-Flast)*(D- Dlast)/(rDispN*(Dmin-0.1*Dmax)-Dlast); %F %D % The slope of the flat part should be lower than the slope of reloading : if (F-Flast)/(D-Dlast)>kr F = Flast + (D-Dlast)*kr; end %end % - reloading & - loading %if D<rDispN*(Dmin-0.1*Dmax) else F = Flast + (D-Dlast)*kr; if F < EM5_MetaBackBone(F1,D1,F2,D2,D) && D<-D1 %&& Dmax>D1 && Dmin<-D1 F = EM5_MetaBackBone(F1,D1,F2,D2,D); end % Case of strong damage due to COUPLING if Flast < EM5_MetaBackBone(F1,D1,F2,D2,Dlast) && (Dmax>0 || Dmin<0) && Dlast<0 % avoid bug at 1st iteration F = Flast - (D-Dlast)*kr*1.2 ;%we use -1/4 of the reloading slope... there might be some better choice %'Flast<MetaBB !!' end end end 169 end % end of displacement decrease------------------------- %% FAILURE CASE ------------------------------------------------------------ %Correct the results if connection has failed if Dmax>Df || D>Df || Dmin<-Df || D<-Df F=0; 'Failure occured in EM5'; if (F-Flast)<-ku*abs(D-Dlast) if D>Dlast F=Flast-1.5*ku*abs(D-Dlast); else F=Flast-0.5*ku*abs(D-Dlast); end end if (F-Flast)>ku*abs(D-Dlast) if D<Dlast F=Flast+1.5*ku*abs(D-Dlast); else F=Flast+0.5*ku*abs(D-Dlast); end end end %% Update internal variables rForceMp = rForceMp +0.5*(F+Flast)*(D-Dlast)* exp(-W/(beta*Wref)) /(Wref*beta); rForceMn = rForceMn +0.5*(F+Flast)*(D-Dlast)* exp(-W/(beta*Wref)) /(Wref*beta); if rForceMp>-1000 'cool'; else 'pb2 in em5_model coco!'; end if rForceMp<0 rForceMp=0; end if rForceMn<0 rForceMn=0; end %update for D (if D>Dmax, rForce increases, thus M decreases): if D>Dmax rForceMp = rForceMp - gamma*rForceMp*(F-Flast)/(EM5_BackBone(F1,D1,F2,D2,D)); end if D<Dmin rForceMn = rForceMn - gamma*rForceMn*(F-Flast)/(EM5_BackBone(F1,D1,F2,D2,D)); end % update dissipated energy (dW=F.dD) W = W+0.5*(F+Flast)*(D-Dlast); % update Dlast, Flast Dlast = D; Flast = F; 170 % update Dmax, Fmax... if D>Dmax Dmax = D; elseif D<Dmin Dmin = D; end function [F,Dmax,Dmin,Dlast,Flast,rForceMp,rForceMn,W] = EM5_ModelVERT(Dlast,Flast,Dmax,Dmin,F1,D1,F2,D2,rDispN,rDispP,rForceN,rForceP,gamma, Wref,rForceMp,rForceMn,beta,kr1,kr2,ku,alpha,W,D,Df) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % EM5 axial behaviour % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Failure displacement : %Df = 0.1; % reloading stiffness %kr = 0.6*F1/D1*(Wconst-(W/Wref)^Wpower); %kr starts at k1, and ends at k2 when W is biiig %kr1 = 1.0*F1/D1; %kr2 = 0.05*F1/D1; %higher alpha => kr tends more slowly towards kr2 %alpha = 0.05; kuN = ku; kr = kr2+(kr1-kr2)/(1+W/(alpha*Wref)); % initialisation if Dmax == 0 && D ==0 F=0; end % update rDisp and rForce (note: the value of rDisp and rForce is not modified outside the procedure): %rDispP = W/Wref*rDispP; rDispP = rDispP; rDispN = rDispP; rForceP = rForceP*(1-rForceMp); %rForceN = rForceP*(1-rForceMn); rForceN = 1. ; %% If displacement increases ----------------------------- if D>=Dlast % + unloading with negative force if Flast<0 F = Flast+(D-Dlast)*ku; if F>0 %F=0; %<<<<<<<<<<<<<<<<< bad in nl end else % + unloading with positive force - flat part (transition between unloading and reloading) if D<=rDispP*Dmax && Dmax>0 && Dmax>1.1*D1 171 F = Flast + (rForceP*EM5_ABackBone(F1,D1,F2,D2,Dmax,rForceP)-Flast)*(D- Dlast)/(rDispP*Dmax-Dlast); % The slope of the flat part should be lower than the slope of reloading... if (F-Flast)/(D-Dlast)>kr F = Flast + (D-Dlast)*kr; end % ... and if the slope is negative, then it shouldn't (this % happens when we try to find the stiffness) if F<Flast F = Flast; end else % + reloading & + loading (this is the case where D is higher than displ of point "A" (flat to reload transition point)) %if D>rDispP*Dmax F = Flast + (D-Dlast)*kr; if F > EM5_MetaBackBone(F1,D1,F2,D2,D) %|| Dmax<D1 F = EM5_MetaBackBone(F1,D1,F2,D2,D); end end end end % end of displacement increase------------------------- %% If displacement decreases ----------------------------- if D<Dlast % - unloading with positive force if Flast>0 F = Flast+(D-Dlast)*kuN; if F<0 %F=0; %<<<<<<<<<<<<<<<<< bad in nl end else % - unloading with negative force on flat part if D>=rDispN*(0-0.1*Dmax) || 1==1 %we need to fall in this category for stiffness check F = Flast + (rForceN*(-30)-Flast)*(D-Dlast)/(0-Dlast); % The slope of the flat part should be lower than the slope of reloading : if (F-Flast)/(D-Dlast)>kr F = Flast + (D-Dlast)*kr; end end end end % end of displacement decrease------------------------- %% If a negative displacement is asked (assessment of tangential stiffness): if D<0 F = Flast + ku*(D-Dlast); end %% FAILURE CASE ------------------------------------------------------------ %Correct the results if connection has failed 172 if Dmax>Df || D>Df || Dmin<-Df || D<-Df F=0; 'Failure occured in EM5_Vert'; if (F-Flast)<-ku*abs(D-Dlast) if D>Dlast F=Flast-2.5*ku*abs(D-Dlast); else F=Flast-0.5*ku*abs(D-Dlast); end end if (F-Flast)>ku*abs(D-Dlast) F=Flast+ku*abs(D-Dlast); end end %% Update internal variables D; Dlast; rForceMp = rForceMp +0.5*(F+Flast)*(D-Dlast)* exp(-W/(beta*Wref)) /(Wref*beta); rForceMn = rForceMn +0.5*(F+Flast)*(D-Dlast)* exp(-W/(beta*Wref)) /(Wref*beta); if rForceMp<0 rForceMp=0; end if rForceMn<0 rForceMn=0; end %update for D (if D>Dmax, rForce increases, thus M decreases): if D>Dmax rForceMp = rForceMp - gamma*rForceMp*(F-Flast)/(EM5_BackBone(F1,D1,F2,D2,D)); end %if D<Dmin % update dissipated energy (dW=F.dD) W = W+0.5*(F+Flast)*(D-Dlast); % update Dlast, Flast Dlast = D; Flast = F; % update Dmax if D>Dmax Dmax = D; end function [ xc,xc1 ] = getConnLocation( nc1,d,H,np ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get location of connectors % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Lm23=[-1/4;1/4;0]; Lm4=[-0.3;-0.1;0.1;0.3]; xc1=[]; for i=1:nc1(1) if nc1(1)==1 xc1=[xc1; 0 ; 0]; elseif nc1(1)<=3 xc1=[xc1; Lm23(i,1)*d ; 0]; elseif nc1(1)==4 xc1=[xc1; Lm4(i,1)*d ; 0]; end end 173 for i=1:nc1(2) if nc1(2)<=3 xc1=[xc1; -d/2 ; Lm23(i,1)*H+H/2 ]; elseif nc1(2)==4 xc1=[xc1; -d/2 ; Lm4(i,1)*H+H/2 ]; end end for i=1:nc1(2) if nc1(2)<=3 xc1=[xc1; d/2 ; Lm23(i,1)*H+H/2 ]; elseif nc1(2)==4 xc1=[xc1; d/2 ; Lm4(i,1)*H+H/2 ]; end end % Note : origin is at the middle of the bottom % ------- % | | % | | % | | % ---O--- xc=xc1; for i=2:np xc = [xc ; xc1]; end end 174 function F = EM5_BackBone(F1,D1,F2,D2,D) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get EM5 backbone % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if D>=0 if D <= D1 F = F1*D/D1; elseif D>D1 F = F1+(F2-F1)*(D-D1)/(D2-D1); end else if D >= -D1 F = F1*D/D1; elseif D<D1 F = -F1+(F2-F1)*(D+D1)/(D2-D1); end end function F = EM5_BackBone(F1,D1,F2,D2,D) if D>=0 F = F1+(F2-F1)*(D-D1)/(D2-D1); else F = -F1+(F2-F1)*(D+D1)/(D2-D1) ; end function F = EM5_ABackBone(F1,D1,F2,D2,D,rF) if D>=-D1 F = (F1/D1)/rF*D; Fb = F1+(F2-F1)*(D-D1)/(D2-D1); if F>Fb F=Fb; end else end 175 function [ XList2,xList,TrialDispListC,TrialDispListFR,RS,L,C1s,C2s,dX_FR2,dXFR,Cd ] = ContactAlgorithm4Db( XList2,XListApplied,xList,K,t,Tilt,TiltList,Round,np,H,d,mu,wp,kRoof,TrialDispListC, TrialDispListFR,RS) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Contact algorithm % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Reset xList(t) (for 2nd round !) xList(:,t)=zeros(3+np,1); RS(1,t)=-0.1; RS(2,t)=-0.1; %% Modify moment due to connections when a new point is used for equilibium if TiltList(1,t-1)*TiltList(1,t)<0 && Round==1 XList2(2,t-1)=-XList2(2,t-1); if t>2 TrialDispListC(1,t-1)=10; end end %% Update XList2 with moment at the rocking point : % use Mcorner = Mmiddle - X*H/2 + Y*Tilt*d/2 : XList2(1,t)=XListApplied(1,t); XList2(2,t)=XListApplied(2,t)-XListApplied(1,t)*H/2; XList2(3+np,t)=XListApplied(3+np,t); for i=1:np XList2(2,t)=XList2(2,t)+XListApplied(2+i,t)*Tilt*d/2; XList2(2+i,t)=XListApplied(2+i,t); end %% Get the trial disp with rocking, uplift and no friction %(will Use TrialDispC to get the sliding direction ) dX2(:,1)=XList2(:,t)-XList2(:,t-1); % XList2(:,t) is just Xapplied with converted moment so far dX_C=zeros(np+3,1); dX_C(1,1)=XList2(1,t)-XList2(1,t-1); dX_C(2,1)=XList2(2,t)-XList2(2,t-1); % We know M and X, we know neither x nor Theta % Loop for uplift ! for i=0:np-1 % i = number of panels that fly in the sky % make new matrice, exchange 1st and 2nd line, as well as the % lines of lifted panels L = K; [L] = ExchangeVar( 1 , L ); [L] = ExchangeVar( 2 , L ); [L] = ExchangeVar( 3+np , L ); for j=1:i if Tilt==-1 % //// [L] = ExchangeVar( 2+j , L ); else % \\\\ [L] = ExchangeVar( np+3-j , L ); end end % solve 176 %C1=(dX dM 0's) %C2=(0 0 dyi's dYj's) C1 = zeros(np+3,1); C1(1,1)=dX_C(1,1); C1(2,1)=dX_C(2,1); C1(3+np,1)=dX2(3+np,1); for j=1:i %Flying panels if Tilt==-1 C1(2+j,1)=dX2(2+j,1); %<<<<---- dY=Yapplied-Y2DOF(t-1)----=----- Correct value of XC ?????????? else C1(3+np-j,1)=dX2(3+np-j,1); end end for j=i+1:np if Tilt==-1 C1(2+j,1)=-xList(2+j,t-1); %<<<<---- y=0 and not dy=0 for the panels on the ground!!! else C1(3+np-j,1)=-xList(3+np-j,t-1); end end C2 = L*C1; % get dx dX dx_C=zeros(3+np,1); %(dx dTh): dx_C(1,1)=C2(1,1); dx_C(2,1)=C2(2,1); dx_C(3+np)=C2(3+np,1); dX_C(3+np)=C1(3+np,1); %(dyi's) and (dYi's): for j=1:np if Tilt==-1 % //// if j<=i dx_C(2+j,1)=C2(2+j,1); dX_C(2+j,1)=C1(2+j,1); else dX_C(2+j,1)=C2(2+j,1); dx_C(2+j,1)=C1(2+j,1); end else if j<1+np-i dX_C(2+j,1)=C2(2+j,1); dx_C(2+j,1)=C1(2+j,1); else dx_C(2+j,1)=C2(2+j,1); dX_C(2+j,1)=C1(2+j,1); end end end if (XList2((2+np-i),t-1)+dX_C(2+np-i)-XListApplied((2+np-i),t))>=0 && Tilt==1 || (XList2((2+1+i),t-1)+dX_C(2+1+i)-XListApplied((2+1+i),t))>=0 && Tilt==-1 %go to next iteration if (i+1)-th reaction is negative break %NOTE : Reaction = X2DOF - Xapplied since X2dof=Xapplied+reac end end TrialDispC = xList(:,t-1)+dx_C; %FrSignList(1,t)=-FrSign; % if t>=10 % dxC 177 % dx_C % end %% Convert this trial displacement into real displacement (ie takes % nonrocking into account if necessary). % Get TrialDispC2-------------------------------------------------------- if Tilt*TrialDispC(2)>=0 %-----------------------------Rocking case dX_C2=dX_C; dx_C2=dx_C; RS(1,t)=1; 'Rocking Trial dispC'; else %-----------------------------------------------------Non-rocking case 'NO Rocking Trial dispC'; RS(1,t)=0; % In this case, we know X and Theta, but we know neither M nor x % => we need to change the matrix problem % We use the procedure ExchangeVar to exchange xi and yi in the problem % y = A x the procedure returns a matrix B : % (y_with xi instead of yi) = B (x_with yi instead of xi) L=K; L = ExchangeVar(1,L); L = ExchangeVar(3+np,L); % V = [dx ; dM ; dY1 ; dY2] = L * [dX ; dTh ; dy1 ; dy2] = L*U % Assemble the U vector : U=zeros(3+np,1); U(1,1)=dX_C(1,1); %dX U(3+np,1)=dX_C(3+np,1); for i=2:np+2 %U(i,1)=0; U(i,1)=-xList(i,t-1); % <<<<<<<<---- no rocking => theta=0 and not dTheta = 0 !!! idem with uplift end V=L*U; %V=[dx ; dM ; dY1 ; dY2 ; dr] dX_C2=zeros(3+np,1); dx_C2=zeros(3+np,1); dx_C2(1,1)=V(1,1); dx_C2(3+np,1)=V(3+np,1); dX_C2(1,1)=U(1,1); dX_C2(3+np,1)=U(3+np,1); for i=2:np+2 dx_C2(i,1)=U(i,1); dX_C2(i,1)=V(i,1); end end TrialDispListC(:,t)=[dx_C;dx_C2]; %% Get the trial disp with rocking, uplift and friction (where friction has its max possible value) : %(Use TrialDispC to get the sliding direction) Now just use force evolution... good ?? dX2(:,1)=XList2(:,t)-XList2(:,t-1); % XList2(:,t) is just Xapplied with converted moment so far Frix = mu*np*wp; %<<< Appromate ! A precise assessment would have to take into account vertical acceleration and connection forces %FrSign = sign(TrialDispC(1)-xList(1,t-1)); 178 %FrSign = sign(XListApplied(1,t)-XListApplied(1,t-1)+kRoof*(xList(3+np,t-1)- xList(1,t-1))); FrSign = sign(dx_C2(1,1)); dXFR=zeros(np+3,1); dXFR(1,1)=XList2(1,t)-FrSign*Frix-XList2(1,t-1); dXFR(2,1)=XList2(2,t)-XList2(2,t-1); % We know M and X, we know neither x nor Theta % Loop for uplift ! for i=0:np-1 % i = number of panels that fly in the sky % make new matrice, exchange 1st and 2nd line, as well as the % lines of lifted panels L = K; [L] = ExchangeVar( 1 , L ); [L] = ExchangeVar( 2 , L ); [L] = ExchangeVar( 3+np , L ); for j=1:i if Tilt==-1 % //// [L] = ExchangeVar( 2+j , L ); else % \\\\ [L] = ExchangeVar( np+3-j , L ); end end % solve %C1=(dX dM 0's) %C2=(dx dTh dyi's dYj's) C1 = zeros(np+3,1); C1(1,1)=dXFR(1,1); C1(2,1)=dXFR(2,1); C1(3+np,1)=dX2(3+np,1); Cd = ones(np+3,1); for j=1:i %Flying panels if Tilt==-1 C1(2+j,1)=dX2(2+j,1); %<<<<---- dY=Yapplied-Y2DOF(t-1)----=----- Correct value of XC ?????????? else C1(3+np-j,1)=dX2(3+np-j,1); end end for j=i+1:np if Tilt==-1 C1(2+j,1)=-xList(2+j,t-1); %<<<<---- y=0 and not dy=0 for the panels on the ground!!! Cd(2+j,1)=0; else C1(3+np-j,1)=-xList(3+np-j,t-1); Cd(3+np-j,1)=0; end end C2 = L*C1; % get dx dX dxFR=zeros(3+np,1); %(dx dTh): dxFR(1,1)=C2(1,1); dxFR(2,1)=C2(2,1); dxFR(3+np)=C2(3+np,1); dXFR(3+np)=C1(3+np,1); %(dyi's) and (dYi's): for j=1:np if Tilt==-1 % //// if j<=i dxFR(2+j,1)=C2(2+j,1); 179 dXFR(2+j,1)=C1(2+j,1); else dXFR(2+j,1)=C2(2+j,1); dxFR(2+j,1)=C1(2+j,1); end else if j<1+np-i dXFR(2+j,1)=C2(2+j,1); dxFR(2+j,1)=C1(2+j,1); else dxFR(2+j,1)=C2(2+j,1); dXFR(2+j,1)=C1(2+j,1); end end end % panNumbSL=i % dX=dXFR(2+1+i) % dx=dxFR(2+1+i) % XListApplied((2+1+i),t) % XList2((2+1+i),t-1) % (XList2((2+1+i),t-1)+dXFR(2+1+i)-XListApplied((2+1+i),t)) % if (XList2((2+np-i),t-1)+dXFR(2+np-i)-XListApplied((2+np-i),t))>=0 && Tilt==1 || (XList2((2+1+i),t-1)+dXFR(2+1+i)-XListApplied((2+1+i),t))>=0 && Tilt==-1 %go to next iteration if (i+1)-th reaction is negative break %NOTE : Reaction = X2DOF - Xapplied since X2dof=Xapplied+reac end end % if i==0 'Oh sinistre attachement aux apparences, que fais tu de l humanite!' t for z=2:np-1 if (XList2(2+z,t-1)+dXFR(2+z)-XListApplied(2+z))<0 %if panel z wants to fly L = K; [L] = ExchangeVar( 1 , L ); [L] = ExchangeVar( 2 , L ); [L] = ExchangeVar( 3+np , L ); for j=z:z [L] = ExchangeVar( 2+j , L ); end % solve %C1=(dX dM 0's) %C2=(dx dTh dyi's dYj's) C1 = zeros(np+3,1); C1(1,1)=dXFR(1,1); C1(2,1)=dXFR(2,1); C1(3+np,1)=dX2(3+np,1); Cd = ones(np+3,1); for j=z:z %Flying panels C1(2+j,1)=dX2(2+j,1); %<<<<---- dY=Yapplied-Y2DOF(t-1)----=----- Correct value of XC ?????????? end for j=1:np if j==z else C1(2+j,1)=-xList(2+j,t-1); %<<<<---- y=0 and not dy=0 for the panels on the ground!!! Cd(2+j,1)=0; end 180 end C2 = L*C1; % get dx dX dxFR=zeros(3+np,1); %(dx dTh): dxFR(1,1)=C2(1,1); dxFR(2,1)=C2(2,1); dxFR(3+np)=C2(3+np,1); dXFR(3+np)=C1(3+np,1); %(dyi's) and (dYi's): for j=1:np if j==z dxFR(2+j,1)=C2(2+j,1); dXFR(2+j,1)=C1(2+j,1); else dXFR(2+j,1)=C2(2+j,1); dxFR(2+j,1)=C1(2+j,1); end end end end end TrialDispFR = xList(:,t-1)+dxFR; %FrSignList(1,t)=-FrSign; % if t>=10 % dxC % dxFR % end %% Convert this trial displacement into real displacement (ie takes % nonrocking into account if necessary). % Get TrialDispFR2-------------------------------------------------------- if Tilt*TrialDispFR(2)>=0 %-----------------------------Rocking case dX_FR2=dXFR; dx_FR2=dxFR; RS(1,t)=1; 'Rocking Trial dispFR'; else %-----------------------------------------------------Non-rocking case 'NO Rocking Trial dispFR'; RS(1,t)=0; % In this case, we know X and Theta, but we know neither M nor x % => we need to change the matrix problem % We use the procedure ExchangeVar to exchange xi and yi in the problem % y = A x the procedure returns a matrix B : % (y_with xi instead of yi) = B (x_with yi instead of xi) L=K; L = ExchangeVar(1,L); L = ExchangeVar(3+np,L); % V = [dx ; dM ; dY1 ; dY2] = L * [dX ; dTh ; dy1 ; dy2] = L*U % Assemble the U vector : U=zeros(3+np,1); U(1,1)=dXFR(1,1); %dX U(3+np,1)=dXFR(3+np,1); Cd = ones(np+3,1); for i=2:np+2 %U(i,1)=0; 181 U(i,1)=-xList(i,t-1); % <<<<<<<<---- no rocking => theta=0 and not dTheta = 0 !!! idem with uplift Cd(i,1)=0; % <<<<<< OK for i=2 ????????? might need to delete this ! end V=L*U; %V=[dx ; dM ; dY1 ; dY2 ; dr] dX_FR2=zeros(3+np,1); dx_FR2=zeros(3+np,1); dx_FR2(1,1)=V(1,1); dx_FR2(3+np,1)=V(3+np,1); dX_FR2(1,1)=U(1,1); dX_FR2(3+np,1)=U(3+np,1); for i=2:np+2 dx_FR2(i,1)=U(i,1); dX_FR2(i,1)=V(i,1); end end %TrialDispListFR(:,t)=[dx_C2(1,1) ; dxFR(1,1) ; dx_FR2(1,1) ;0; dxFR]; %[dxFR;dX_FR2(1,1);dXFR(1,1);XList2(1,t);XList2(1,t-1)]; %[dxFR;dx_FR2]; C1s=C1; C2s=C2; Ks=K; %% Compare the 2 displacements in order to determine whether sliding occurs % or not % Case 1 : Sliding occurs -> use TrialDispFR2 (in which rocking is already assessed) if dx_FR2(1,1)*FrSign>0 %if dx_FR(1)*dx_C(1)>0 xList(:,t)=xList(:,t-1)+dx_FR2(:,1); XList2(:,t)=XList2(:,t-1)+dX_FR2(:,1); RS(2,t)=1; 'Sliiiiiiide'; % Case 2 : No sliding -> X is unknown else 'No slide'; RS(2,t)=0; for i=0:np-1 % i = number of panels that fly in the sky % make new matrice, exchange 1st and 2nd line, as well as the % lines of lifted panels % C2 = [dX ; dTh ; dY1 ... dy_np] = L * [dx ; dM ; dy1...dY_np] = L*C1 if Tilt=+1 L = K; [L] = ExchangeVar( 2 , L ); [L] = ExchangeVar( 3+np , L ); for j=1:i if Tilt==-1 % //// [L] = ExchangeVar( 2+j , L ); else % \\\\ [L] = ExchangeVar( np+3-j , L ); end end % solve %C1=(dX dM 0's) %C2=(0 0 dyi's dYj's) C1 = zeros(3+np,1); Cd = ones(np+3,1); C1(1,1)=0; Cd(1,1)=0; C1(2,1)=XList2(2,t)-XList2(2,t-1); %------- dXC(2,1)? C1(3+np,1)=dX2(3+np,1); 182 for j=1:i %Flying panels if Tilt==-1 C1(2+j,1)=dX2(2+j,1); %<<<<---- dY=Yapplied-Y2DOF(t-1)----=----- Correct value of XC ?????????? else C1(3+np-j,1)=dX2(3+np-j,1); end end for j=i+1:np if Tilt==-1 C1(2+j,1)=-xList(2+j,t-1); %<<<<---- y=0 and not dy=0 for the panels on the ground!!! Cd(2+j,1)=0; else C1(3+np-j,1)=-xList(3+np-j,t-1); Cd(3+np-j,1)=0; end end C2 = L*C1; % get dx dX dx=zeros(np+3,1); dX=zeros(np+3,1); %(dx dTh): dx(1,1)=0; dx(2,1)=C2(2,1); dX(1,1)=C2(1,1); dX(2,1)=C1(2,1); dx(3+np)=C2(3+np,1); dX(3+np)=C1(3+np,1); %(dyi's) and (dYi's): for j=1:np if Tilt==-1 % //// if j<=i dx(2+j,1)=C2(2+j,1); dX(2+j,1)=C1(2+j,1); else dX(2+j,1)=C2(2+j,1); dx(2+j,1)=C1(2+j,1); end else if j<1+np-i dX(2+j,1)=C2(2+j,1); dx(2+j,1)=C1(2+j,1); else dx(2+j,1)=C2(2+j,1); dX(2+j,1)=C1(2+j,1); end end end % panNumb=i % dX(2+1+i) % XListApplied((2+1+i),t) % XList2((2+1+i),t-1) % (XList2((2+1+i),t-1)+dX(2+1+i)-XListApplied((2+1+i),t)) % if (XList2((2+np-i),t-1)+dX(2+np-i)-XListApplied((2+np-i),t))>=0 && Tilt==1 || (XList2((2+1+i),t-1)+dX(2+1+i)-XListApplied((2+1+i),t))>=0 && Tilt==-1 %go to next iteration if (i+1)-th reaction is negative break %NOTE : Reaction = X2DOF - Xapplied since X2dof=Xapplied+reac end 183 end if i==0 'heyaaa!' t for z=2:np-1 if (XList2(2+z,t-1)+dXFR(2+z)-XListApplied(2+z))<0 %if panel z wants to fly L = K; [L] = ExchangeVar( 2 , L ); [L] = ExchangeVar( 3+np , L ); for j=z:z [L] = ExchangeVar( 2+j , L ); end % solve %C1=(dX dM 0's) %C2=(dx dTh dyi's dYj's) C1 = zeros(3+np,1); Cd = ones(np+3,1); C1(1,1)=0; Cd(1,1)=0; C1(2,1)=XList2(2,t)-XList2(2,t-1); %------- dXC(2,1)? C1(3+np,1)=dX2(3+np,1); for j=z:z %Flying panels C1(2+j,1)=dX2(2+j,1); %<<<<---- dY=Yapplied-Y2DOF(t-1)----=----- Correct value of XC ?????????? end for j=1:np if j==z else C1(2+j,1)=-xList(2+j,t-1); %<<<<---- y=0 and not dy=0 for the panels on the ground!!! Cd(2+j,1)=0; end end C2 = L*C1; % get dx dX dx=zeros(np+3,1); dX=zeros(np+3,1); %(dx dTh): dx(1,1)=0; dx(2,1)=C2(2,1); dX(1,1)=C2(1,1); dX(2,1)=C1(2,1); dx(3+np)=C2(3+np,1); dX(3+np)=C1(3+np,1); %(dyi's) and (dYi's): for j=1:np if j==z dx(2+j,1)=C2(2+j,1); dX(2+j,1)=C1(2+j,1); else dX(2+j,1)=C2(2+j,1); dx(2+j,1)=C1(2+j,1); end end end end end % 184 xList(:,t)=xList(:,t-1)+dx; XList2(:,t)=XList2(:,t-1)+dX; % Then we check Theta. If Theta.Tilt<0, there is NO ROCKING, and we % need to correct the result -> NO MOTION if Tilt*xList(2,t)<0 RS(1,t)=0; RS(2,t)=0; 'NO MOTION !!!'; %Panel are still, but roof can move !! %Get L : modified K ---- Some K has been calculated before : good %enough since we just need the last column, which does not degrade: L = K; [L] = ExchangeVar( 3+np , L ); %Get C1 : vector of known displacements and forces C1 = zeros(3+np,1); Cd = ones(np+3,1); Cd(1,1)=0; for i=2:2+np C1(i,1)=-xList(i,t-1); Cd(i,1)=0; end C1(3+np,1)=dX2(3+np,1); %Solve the matrix equation to get C2=(dX 0 dyi's dYj's) C2 = L*C1; %Extract forces and displacements from C1 and C2 % get dx dX dx=zeros(np+3,1); dX=zeros(np+3,1); %(dx dTh): for i=1:2+np dX(i,1)=C2(i,1); dx(i,1)=C1(i,1); end dx(3+np)=C2(3+np,1); dX(3+np)=C1(3+np,1); XList2(:,t)=XList2(:,t-1)+dX; xList(:,t)=xList(:,t-1)+dx; else RS(1,t)=1; 'ok'; end end end 185 function [ E ] = ExchangeVar( k , D ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Exchange variables of a matrx equation % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Exchange m-th var of matrix problem D (D should be square) sizeD=size(D); E = zeros(sizeD); for i=1:sizeD(1) for j=1:sizeD(1) if i==k if j==k E(i,j)=1/D(k,k); else E(i,j)=-D(k,j)/D(k,k); end else if j==k E(i,j)=D(i,k)/D(k,k); else E(i,j)=D(i,j)-D(k,j)*D(i,k)/D(k,k); end end end end end 186 function [ k_p,k_n ] = getTangents_H( p,t,Df ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Assess tangent stiffness in horizontal direction % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Returns tangents in case of : % - Positive increment (kp) % - Negative increment (kn) % Displacement step: dd=0.0000000001; % Get the force in case of a small displacement increase : [Fp,~,~,~,~,~,~]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8,1),p (9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p(19,1 ),p(20,1),p(21,1),p(22,1),p(1,1)+dd,Df); % Get the force in case of a small displacement decrease : [Fn,~,~,~,~,~,~]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8,1),p (9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p(19,1 ),p(20,1),p(21,1),p(22,1),p(1,1)-dd,Df); % Get the current force in the connection : F = p(2,1); % p(1,1) % dd % Fn % Fp % Calc stiffnesses : k_p=(Fp-F)/dd; k_n=(Fn-F)/-dd; [k_p,k_n]; end function [ k_p,k_n ] = getTangents_V( p,dfV ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Tangent stiffness of EM5 in vertical direction % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Returns tangents in case of : % - Positive increment (kp) % - Negative increment (kn) % Displacement step: dd=0.000001; % Get the force in case of a small displacement increase : [Fp,~,~,~,~,~,~]=EM5_ModelVERT(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8, 1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p( 19,1),p(20,1),p(21,1),p(22,1),p(1,1)+dd,dfV); % Get the force in case of a small displacement decrease : [Fn,~,~,~,~,~,~]=EM5_ModelVERT(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8, 1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p( 19,1),p(20,1),p(21,1),p(22,1),p(1,1)-dd,dfV); 187 % Get the current force in the connection : F = p(2,1); % Calc stiffnesses : k_p=(Fp-F)/dd; k_n=(Fn-F)/-dd; end 188 function [ Tilt,TiltList ] = getTilt( XListAPPLIED,Mcon,TiltList,H,t ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Determine rocking case (static) % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Stiffness of // Theta<0 -> Tilt=-1 %Stiffness of \\ Theta>0 -> Tilt=+1 % Calc moment to due to horizontal force + moment to determine what matrix % we should use ----------------------------------<<<<<<<<<<<<<<<<<<<<<<<<<<<< then we will use dD and prvious Theta SumM = 0; for i=1:1 SumM = SumM-XListAPPLIED(3*(i-1)+1,t)*H/2+XListAPPLIED(3*(i-1)+3,t)+Mcon; % !! Sign of Mcon ?? +, -, both ? [SumM , XListAPPLIED(3*(i-1)+1,t)*H/2 , XListAPPLIED(3*(i-1)+3,t)]; end if XListAPPLIED(2,t-1)>0 %Case Previous theta>0 -> we use \\ Tilt = 1; TiltList(t)=10; elseif XListAPPLIED(2,t-1)<0 %Case Previous theta<0 -> we use // Tilt = -1; TiltList(t)=-10; else %Case Previous theta=0 => use force if SumM >= 0 Tilt = 1; TiltList(t)= 1; else Tilt = -1; TiltList(t)= -1; end end end 189 function [ Tilt,TiltList,SumM ] = getTiltDyn2( XListAPPLIED_dyn,Xapplied,xList,TiltList,HX,H,kRoof,np,PrevSum,t ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get rocking case in dynamics % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Stiffness of // Theta<0 -> Tilt=-1 %Stiffness of \\ Theta>0 -> Tilt=+1 SumM=-XListAPPLIED_dyn(1,t)*HX+XListAPPLIED_dyn(2,t)-(xList(3+np,t-1)-xList(1,t- 1)+H*xList(2,t-1))*kRoof*H; %<< at middle MdynForces=-(XListAPPLIED_dyn(1,t)*H/2-Xapplied(1,1))*H/2 if xList(2,t-1)>0 %Case Previous theta>0 -> we use \\ Tilt = 1; TiltList(t)=10; elseif xList(2,t-1)<0 %Case Previous theta<0 -> we use // Tilt = -1; TiltList(t)=-10; else %Case Previous theta=0 => use force if SumM >= 0 Tilt = 1; TiltList(t)= 1; else Tilt = -1; TiltList(t)= -1; end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Plot IDA results % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %SaTrList=SaTrList0(3:39,1); SaTrList=SaTrList0; Ngm=10; Ncf=7; Ndp=1%length(MuKrFy); %% Plot max sliding disp figure; for cf=3:3 Sll=zeros(Ndp,length(SaTrList)); for dp=1:Ndp for sc=1:length(SaTrList) 190 Sll(dp,sc)=median(MaxSl_l(sc,cf,dp,:)/9.8); end end %subplot(1,Ncf,cf); switch cf case 1 %plot(Sll,SaTrList,'b') plot(Sll,SaTrList/9.8,'-or','LineWidth',2) case 2 plot(Sll,SaTrList,'-+g','MarkerSize',4) case 3 plot(Sll,SaTrList/9.8,'-or','LineWidth',2) case 4 plot(Sll,SaTrList,'c','MarkerSize',4) case 5 plot(Sll,SaTrList,'-or','MarkerSize',4) case 6 plot(Sll,SaTrList,'-.k') case 7 plot(Sll,SaTrList,'--k') end hold on %xlabel(strcat('Sa(Tr) (m.s-2)',Tcf)) ylabel('Sa(Tr) (g)','FontSize',14) xlabel('Max Sliding drift','FontSize',14) %grid; end set(gca,'FontSize',14) %plot(0.0165*ones(2,1),[0;14],'-k','LineWidth',2) %axis([0 0.02 0 13]) hold on %% Plot max rocking disp for cf=3:3 Sll=zeros(Ndp,length(SaTrList)); for dp=1:Ndp for sc=1:length(SaTrList) Sll(dp,sc)=RCtoFC(median(MaxRD_l(sc,cf,dp,:)))/H; end end %subplot(1,Ncf,cf); switch cf case 1 %plot(Sll,SaTrList,'b') plot(Sll,SaTrList/9.8,'-og','LineWidth',2) case 2 plot(Sll,SaTrList,'-+g','MarkerSize',4) case 3 plot(Sll,SaTrList/9.8,'-og','LineWidth',2) case 4 plot(Sll,SaTrList,'c','MarkerSize',4) case 5 plot(Sll,SaTrList,'-or','MarkerSize',4) case 6 plot(Sll,SaTrList,'-.k') case 7 plot(Sll,SaTrList,'--k') end hold on %xlabel(strcat('Sa(Tr) (m.s-2)',Tcf)) ylabel('Sa(Tr) (g)','FontSize',14) xlabel('Max drift','FontSize',14) 191 %legend('Normal','Half Stregth','Double Strength','Double degradation rate','Half degradation rate','120% ductility','150% ductility','Location','NorthWest') %grid; end set(gca,'FontSize',14) axis([0 0.04 0 1.2]) legend('Sliding','Rocking') %% Plot max roof shear figure;%subplot(1,3,3); for cf=3:3 Sll=zeros(Ndp,length(SaTrList)); for dp=1:Ndp for sc=1:length(SaTrList) Sll(dp,sc)=kRoof*median(MaxRoofDDList(sc,cf,dp,:))/2000; end end %subplot(1,Ncf,cf); switch cf case 1 %plot(Sll,SaTrList,'b') plot(Sll,SaTrList/9.8,'-ob','LineWidth',2) case 2 plot(Sll,SaTrList,'-+g','MarkerSize',4) case 3 plot(Sll,SaTrList/9.8,'-ob','LineWidth',2) case 4 plot(Sll,SaTrList,'c','MarkerSize',4) case 5 plot(Sll,SaTrList,'-or','MarkerSize',4) case 6 plot(Sll,SaTrList,'-.k') case 7 plot(Sll,SaTrList,'--k') end hold on %xlabel(strcat('Sa(Tr) (m.s-2)',Tcf)) ylabel('Sa(Tr) (m.s-2)','FontSize',14) xlabel('Max roof shear ratio','FontSize',14) %legend('Normal','Half Stregth','Double Strength','Double degradation rate','Half degradation rate','120% ductility','150% ductility','Location','NorthWest') %grid; end set(gca,'FontSize',14) axis([0 2 0 1.2]) 192 function [ lF1p ] = getLoadList_Trapeze(Xmax,HX,H,nts,ncyc,wp,np,coef,ntsC) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get a trapezoidal load load % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ntsW=floor(nts/10); if ntsW<=4 ntsW=5; end lh = zeros(1,ntsW-1); lv1 = [0:(-wp/(ntsW-2)):-wp]; lv=[]; for i=1:np lv = [lv ; lv1 ] ; end lm = zeros(1,ntsW-1); lr = zeros(1,ntsW-1); Lweight = [lh ; lm ; lv ; lr]; %Max moment (about center of the panel) due to horizontal force: Mmax=-Xmax*(HX-H/2); %list of external forces (x force, y force, moment at center of panel) at each time step lh = [0:(Xmax/(nts-2)):Xmax]; lv = -wp*ones(1,nts-1); lm = [0:(Mmax/(nts-2)):Mmax]; lr = coef*lh; %List of forces : LV = []; LV2 = []; wpNP=[]; for i=1:np LV = [LV ; lv ] ; LV2 = [LV2 ; lv ] ; wpNP=[wpNP ; -wp] end lF1p_1=[ lh ; lm ; LV ; lr ]; lF1p_C=[ Xmax ; Mmax ; wpNP ; coef*Xmax]*ones(1,ntsC); lF1p_2=[ (Xmax*ones(1,nts-1)-lh) ; (Mmax*ones(1,nts-1)-lm) ; LV2 ;coef*(Xmax*ones(1,nts-1)-lh)]; lF1p=[lF1p_1 lF1p_C lF1p_2]; %lF1p=[lF1p_1]; % %Uncomment to add negative loading: Xmax = -Xmax; Mmax=-Xmax*(HX-H/2); lh = [0:(Xmax/(nts-2)):Xmax]; lv = -wp*ones(1,nts-1); lm = [0:(Mmax/(nts-2)):Mmax]; lr = coef*lh; LV = []; for i=1:np LV = [LV ; lv] ; end 193 lF1p_Cn = [ Xmax ; Mmax ; wpNP ; coef*Xmax]*ones(1,ntsC); %List of forces for one panel : lF1p_1=[ lh ; lm ; LV ; lr ]; lF1p_2=[ (Xmax*ones(1,nts-1)-lh) ;(Mmax*ones(1,nts-1)-lm) ; LV ;coef*(Xmax*ones(1,nts-1)-lh)]; lF1p22=[lF1p_1 lF1p_Cn lF1p_2]; lF1p=[lF1p lF1p22]; %Regular positive loading (without zero) : Xmax = -Xmax; Mmax=-Xmax*(HX-H/2); lh = [0:(Xmax/(nts-2)):Xmax]; lv = -wp*ones(1,nts-1); lm = [0:(Mmax/(nts-2)):Mmax]; lr = coef*lh; LV = []; for i=1:np LV = [LV ; lv ] ; end %List of forces for one panel : lF20=[ lh ; lm ; LV ;lr ]; lF21=[ (Xmax*ones(1,nts-1)-lh) ;(Mmax*ones(1,nts-1)-lm) ; LV ;coef*(Xmax*ones(1,nts- 1)-lh) ]; lF22=[lF20 lF1p_C lF21]; %One full +/- loading lFF = [lF22 lF1p22]; for i=2:ncyc lF1p=[lF1p lFF]; end lF1p = [Lweight lF1p]; end 194 function [ lF1p ] = getLoadList_SuddenLoad(Xmax,HX,H,nts,ncyc,wp,np,coef,ntsC) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get a step load % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ntsW=floor(nts/10); if ntsW<=4 ntsW=5; end lh = zeros(1,ntsW-1); lv1 = [0:(-wp/(ntsW-2)):-wp]; lv=[]; for i=1:np lv = [lv ; lv1 ] ; end lm = zeros(1,ntsW-1); lr = zeros(1,ntsW-1); Lweight = [lh ; lm ; lv ; lr]; %Max moment (about center of the panel) due to horizontal force: Mmax=-Xmax*(HX-H/2); %list of external forces (x force, y force, moment at center of panel) at each time step lh = 0; %[0:(Xmax/(nts-2)):Xmax]; lv = -wp; % -wp*ones(1,nts-1); lm = 0; % [0:(Mmax/(nts-2)):Mmax]; lr = 0; % coef*lh; %List of forces : LV = []; LV2 = []; wpNP=[]; for i=1:np LV = [LV ; lv ] ; LV2 = [LV2 ; lv ] ; wpNP=[wpNP ; -wp] end lF1p_1=[ lh ; lm ; LV ; lr ]; lF1p_C=[ Xmax ; Mmax ; wpNP ; coef*Xmax]*ones(1,ntsC); lF1p_2=[ (Xmax*ones(1,nts-1)-lh) ; (Mmax*ones(1,nts-1)-lm) ; LV2 ;coef*(Xmax*ones(1,nts-1)-lh)]; lF1p=[Lweight lF1p_C ]; %lF1p=[lF1p_1]; % end 195 function [ lF1p ] = getLoadList_Sinusoid(T0,Xl,HX,H,wp,np,coef,dt,ncyc) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Get a sinusoidal load % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ntsW=10; if ntsW<=4 ntsW=5; end lh = zeros(1,ntsW-1); lv1 = [0:(-wp/(ntsW-2)):-wp]; lv=[]; for i=1:np lv = [lv ; lv1 ] ; end lm = zeros(1,ntsW-1); lr = zeros(1,ntsW-1); lF1p = [lh ; lm ; lv ; lr]; Tl=[0:dt:ncyc*T0+10*dt]; Tl1=[0:dt:ncyc*T0]; for i=1:length(Xl) Xmax=Xl(i); Rmax=coef*Xmax; %Max moment (about center of the panel) due to horizontal force: Mmax=-Xmax*(HX-H/2); %list of external forces (x force, y force, moment at center of panel) at each time step lh = Xmax*sin(2*pi*Tl1/T0); lv = -wp*ones(1,length(Tl1)); lm = Mmax*sin(2*pi*Tl1/T0); lr = Rmax*sin(2*pi*Tl1/T0); %List of forces : LV = []; LV2 = []; wpNP=[]; for i=1:np LV = [LV ; lv ] ; LV2 = [LV2 ; lv ] ; wpNP=[wpNP ; -wp] ; end lF1p_1=[ lh ; lm ; LV ; lr ]; lF1p = [lF1p lF1p_1]; end end 196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Display Modes % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Run after after having compet stiffness and mass matrix from Main % Time and amplitude of display N0=80; ampl=300; Tilt=-1; R = eye(3+np); R(1,2)=-H/2; for j=3:2+np R(j,2)=Tilt*dc; end T=transpose(R); TiltList=-ones(1000); [VectPPchelou,ValPP]=eig((T\K0/R)/M); NatPeriods = 2*pi*ValPP^-0.5 VectPP=M\VectPPchelou %<<< coord at middle! xMode1=zeros(3+np,N0); Mode1=VectPP(:,3); for t=1:N0 xMode1(:,t)=sin(2*pi*t/(0.5*N0))*Mode1; end x=zeros(5); D66 = (d+0.025)*ones(1,5); H66 = H/2*ones(1,5); figure('color','white'); axis on, axis equal %title('Tilt-up panels Modes','Color',[.3 0 0]) hold on axis([-1.5*d/2 1.5*d/2+(np-1)*d+2 -0.5*H 1.5*H]) xxx = [-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy = [x(3) x(3) H+x(3) H+x(3) x(3)] + x(2)*d/2*([-1 1 1 -1 - 1]+TiltList(1)*ones(1,5)); ht = area(xxx,yyy,'facecolor',[.6 0 0]); set(ht,'XDataSource','xxx') set(ht,'YDataSource','yyy') if np>1 xxx2 = [-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy2 = [x(4) x(4) H+x(4) H+x(4) x(4)] + x(2)*d/2*([-1 1 1 -1 - 1]+TiltList(1)*ones(1,5)); ht2 = area(xxx2,yyy2,'facecolor',[.18 .72 0]); set(ht2,'XDataSource','xxx2') set(ht2,'YDataSource','yyy2') end if np>2 xxx3 = [-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(1)*H*[0 0 -1 -1 0]; 197 yyy3 = [x(5) x(5) H+x(5) H+x(5) x(5)] + x(2)*d/2*([-1 1 1 -1 - 1]+TiltList(1)*ones(1,5)); ht3 = area(xxx3,yyy3,'facecolor',[.2 .2 .9]); set(ht3,'XDataSource','xxx3') set(ht3,'YDataSource','yyy3') end if np>3 xxx4 = [-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(1)*H*[0 0 -1 -1 0]; yyy4 = [x(6) x(6) H+x(6) H+x(6) x(6)] + x(2)*d/2*([-1 1 1 -1 - 1]+TiltList(1)*ones(1,5)); ht4 = area(xxx4,yyy4,'facecolor',[.0 .85 .85]); set(ht4,'XDataSource','xxx4') set(ht4,'YDataSource','yyy4') end xR = -d/2+x(np+3); Dims=[xR,H+0.6,np*d,H/20]; hR = rectangle('Position',Dims); set(hR,'Position',Dims); xq = [0 1]; yq = [0 0.5]; lq=plot(xq,yq,'k','LineWidth',2); for j = 1:N0 x=ampl*xMode1(:,j); xxx = [-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy = [x(3) x(3) H+x(3) H+x(3) x(3)] + x(2)*d/2*([-1 1 1 -1 - 1]+sign(TiltList(j))*ones(1,5)); if np>1 xxx2 = D66 +[-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy2 = [x(4) x(4) H+x(4) H+x(4) x(4)] + x(2)*d/2*([-1 1 1 -1 - 1]+sign(TiltList(j))*ones(1,5)); end if np>2 xxx3 = 2*D66 +[-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy3 = [x(5) x(5) H+x(5) H+x(5) x(5)] + x(2)*d/2*([-1 1 1 -1 - 1]+sign(TiltList(j))*ones(1,5)); end if np>3 xxx4 = 3*D66 +[-d/2+x(1) d/2+x(1) d/2+x(1) -d/2+x(1) -d/2+x(1)] + x(2)*H*[0 0 -1 -1 0]; yyy4 = [x(6) x(6) H+x(6) H+x(6) x(6)] + x(2)*d/2*([-1 1 1 -1 - 1]+sign(TiltList(j))*ones(1,5)); end Dims=[-d/2+x(np+3),H+0.6,np*d,H/20]; set(hR,'Position',Dims); xq1=[-d/4 -d/4+x(1) -d/4 d/4 d/4+x(1) d/4]; yq1=[0 x(3)-3*d/4*x(2) 0 0 x(3)-d/4*x(2) 0]; yq2=[0 x(4)-3*d/4*x(2) 0 0 x(4)-d/4*x(2) 0]; yq3=[0 x(5)-3*d/4*x(2) 0 0 x(5)-d/4*x(2) 0]; yq4=[0 x(6)-3*d/4*x(2) 0 0 x(6)-d/4*x(2) 0]; xq2=xq1+d*ones(1,6); xq3=xq1+2*d*ones(1,6); xq4=xq1+3*d*ones(1,6); set(lq,'x',[xq1 xq2 xq3 xq4],'y',[yq1 yq2 yq3 yq4]) refreshdata 198 drawnow end 199 function [ accInt,dt ] = InterpolL( accList,dt,n,ntMAX ) % Get more points in the ground motion %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Interpolate steps between ground motion steps % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nt = min(n*length(accList),ntMAX); accInt=zeros(1,nt); dt=dt/n; for i=1:length(accList)-1 for j=1:n if (i-1)*n+j<=ntMAX accInt(1,(i-1)*n+j)=accList(i)+(j-1)/(n)*(accList(i+1)-accList(i)); end end end end 200 function [ Ft,Fp,Fl ] = GetConForces( Param,xt,nc1,np,RoofP,H,A,B,S,t,degr,P,dfH,dfV,dfS ) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Compute internal forces % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xi=A*xt(1:2+np,1); n2=2*(nc1(1)+2*nc1(2)); Fl=zeros(np*n2,1); %Fext=zeros(3,1); for j=1:np for i=1:nc1(1) %% Bottom********************************************************************* % Bottom vertical connector ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); [p(1:10,1) ; p(21:22,1)]; % get displacement of the connector for displacement xt : D=xi((j-1)*n2+2*(i-1)+2,1); % Get the force : [Fv,~,~,~,~,~,~,~]=EM5_ModelVERT(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p( 8,1),p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1), p(19,1),p(20,1),p(21,1),p(22,1),D,dfV); Fl((j-1)*n2+2*(i-1)+2,1)=Fv; % Bottom horizontal connector ------------ %------------------------------------------------Effect of coupling % Note : Param is changed but not saved !! nn=3; if degr==1 && t-nn*floor(t/nn)==0 % % Wref(14) = 21.1-123*Dmax_v(3) % newWref=21.1-123*Param(3,(j-1)*n2+2*(i-1)+2); % % F1(5) = 207-1075*Dmax_v; % newF1=subplus(207-1075*Param(3,(j-1)*n2+2*(i-1)+2)); % % F2(7) = 254-1145*Dmax_v; % newF2=subplus(254-1145*Param(3,(j-1)*n2+2*(i-1)+2)); % % % Change kr1 (line 18), kr2(19), ku(20) (which are proportional to F1/D1) % %equation used : newkr = oldkr*newF1/oldF1 % Param(18,(j-1)*n2+2*(i-1)+1)=Param(18,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); % Param(19,(j-1)*n2+2*(i-1)+1)=Param(19,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); % Param(20,(j-1)*n2+2*(i-1)+1)=Param(20,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); % % % Change F1,F2,Wref % Param(14,(j-1)*n2+2*(i-1)+1)=newWref; % Param(5,(j-1)*n2+2*(i-1)+1)=newF1; % Param(7,(j-1)*n2+2*(i-1)+1)=newF2; % Wref(14) = (21.1-123*Dmax_v(3))*WrefTrial/WrefHabituel %newWref=P(1,4)/21.1*subplus(21.1-123*Param(3,(j-1)*n2+2*(i-1)+2)); newWref=P(1,4)/21.1*max((21.1-123*Param(3,(j-1)*n2+2*(i-1)+2)),0); % F1(5) = 207-1075*Dmax_v; %newF1=P(1,1)/207*subplus(207-1075*Param(3,(j-1)*n2+2*(i-1)+2)); newF1=P(1,1)/207*max((207-1075*Param(3,(j-1)*n2+2*(i-1)+2)),0); 201 % F2(7) = 254-1145*Dmax_v; %newF2=P(1,3)/254*subplus(254-1145*Param(3,(j-1)*n2+2*(i-1)+2)); newF2=P(1,3)/254*max((254-1145*Param(3,(j-1)*n2+2*(i-1)+2)),0); % Change kr1 (line 18), kr2(19), ku(20) (which are proportional to F1/D1) %equation used : newkr = oldkr*newF1/oldF1 Param(18,(j-1)*n2+2*(i-1)+1)=Param(18,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); Param(19,(j-1)*n2+2*(i-1)+1)=Param(19,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); Param(20,(j-1)*n2+2*(i-1)+1)=Param(20,(j-1)*n2+2*(i- 1)+1)*newF1/Param(5,(j-1)*n2+2*(i-1)+1); % Change F1,F2,Wref Param(14,(j-1)*n2+2*(i-1)+1)=newWref; Param(5,(j-1)*n2+2*(i-1)+1)=newF1; Param(7,(j-1)*n2+2*(i-1)+1)=newF2; % % axial failure => shear failure !! % if Dmax>=dfV % Param(3,(j-1)*n2+2*(i-1)+1)=0.5; %<< we say that max axial displ is 0.5m, so failed behaviour will be used at next time step. % end end %---------------------------------------------------end of Coupling %Now change all the other parameters : p=Param(:,(j-1)*n2+2*(i-1)+1); % Need to get D of this DOF from xi (calculated from xList) ! D=xi((j-1)*n2+2*(i-1)+1,1); % Get the force : [Fh,~,~,~,~,~,~,~]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8,1) ,p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p(19 ,1),p(20,1),p(21,1),p(22,1),D,dfH); Fl((j-1)*n2+2*(i-1)+1,1)=Fh; end for i=nc1(1)+1:nc1(1)+nc1(2) if j==1 %% left edge *********************************************************** % no need - look in older versions if needed else %% panel-panel left (change Param even if an elastic connector is used- >simpler)********** % vertical connector ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); % get D D=-xi((j-1)*n2+2*(i-1)+2,1)+xi((j-2)*n2+2*(i+nc1(2)-1)+2,1); [F,~,~,~,~,~,~,~]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8,1), p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p(19, 1),p(20,1),p(21,1),p(22,1),D,dfS); Fl((j-1)*n2+2*(i-1)+2,1)=-F; end end for i=nc1(1)+nc1(2)+1:nc1(1)+2*nc1(2) if j==np %% right edge************************************************** % no need - look in older versions if needed else %% panel-panel right*************************************************** 202 % vertical connector ------------- p=Param(:,(j-1)*n2+2*(i-1)+2); D=-xi((j-1)*n2+2*(i-1)+2,1)+xi(j*n2+2*(i-nc1(2)-1)+2,1); [F,~,~,~,~,~,~,~]=EM5_Model(p(1,1),p(2,1),p(3,1),p(4,1),p(5,1),p(6,1),p(7,1),p(8,1), p(9,1),p(10,1),p(11,1),p(12,1),p(13,1),p(14,1),p(15,1),p(16,1),p(17,1),p(18,1),p(19, 1),p(20,1),p(21,1),p(22,1),D,dfS); Fl((j-1)*n2+2*(i-1)+2,1)=-F; end end end %% Calculate the force due to panel connectors: Fp=[S*B*Fl ; 0 ]; %% Forces due to roof : D=xt(3+np,1)-xt(1,1)+H*xt(2,1); [Fspring,~,~] = BilinearHysteretic(RoofP(1,1),RoofP(2,1),RoofP(3,1),RoofP(4,1),RoofP(5,1),D); kRoof=42600; K2 = zeros(3+np,3+np); K2(1,1)=kRoof; K2(1,2)=-H*kRoof; K2(1,3+np)=-kRoof; K2(2,1)=-H*kRoof; K2(2,2)=H^2*kRoof; K2(2,3+np)=H*kRoof; K2(3+np,1)=-kRoof; K2(3+np,2)=H*kRoof; K2(3+np,3+np)=kRoof; Froof1 = -K2*xt; Froof = zeros(3+np,1); Froof(1,1)=-Fspring; Froof(2,1)=H*Fspring; Froof(3+np,1)=Fspring; Froof=-Froof; if Froof==Froof1 'ok'; else 'RoofNL'; end %% TOTAL FORCES DUE TO CONNECTIONS Ft = Fp+Froof; end 203 function [F,Dlast,Flast] = BilinearHysteretic(Dlast,Flast,k0,Fy,k1,D) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Bilinear histeretic material for the roof % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Bilinear histeretic material for the roof %k0 = elastic stiffness %Fy = yielding force %k1 = postyield stiffness F = Flast+k0*(D-Dlast); if F>Fy+k1*(D-Fy/k0) F=Fy+k1*(D-Fy/k0); elseif F<-Fy+k1*(D+Fy/k0) F=-Fy+k1*(D+Fy/k0); end Flast=F; Dlast=D; end
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Numerical modeling for the seismic response of concrete...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Numerical modeling for the seismic response of concrete tilt-up buildings Tellier, Xavier 2013
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Numerical modeling for the seismic response of concrete tilt-up buildings |
Creator |
Tellier, Xavier |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | Concrete tilt-up building is a prevalent construction technique used for industrial and commercial applications in North America. This construction technique offers many significant advantages over conventional cast-in-place construction. This includes the reduction in construction time and the amount of formworks. Despite the large array of buildings that has been constructed using such technique, the nonlinear behaviour of the concrete tilt-up buildings is still not well understood. The nonlinear behaviour of the concrete tilt-up building has been studied in this thesis. The nonlinear response of the concrete tilt-up building is largely affected by the nonlinear behaviour of the connectors between the panels and the slab, and between the panels. Past researches have been conducted to experimentally examine the nonlinear behaviour of the tilt-up panel connectors. The experimental results were used in this thesis to develop an empirical numerical model capable of reproducing the force-deformation response of the tilt-up connectors under combined axial and shear deformation. The numerical model takes the shear strength and stiffness degradation into account after axial cycles of inelastic deformation. A finite-element software was developed specifically to study the nonlinear static and dynamic behaviour of concrete tilt-up buildings. A typical tilt-up building designed in the study of Olund (2009) was modeled. Incremental dynamic analysis was performed using the developed finite element software to assess the seismic performance of the prototype tilt-up building. The results of the incremental dynamic analysis provided valuable information to understand the nonlinear behaviour of the concrete tilt-up buildings under seismic load. Detailed parametric studies were carried out to examine the nonlinear behaviour of tilt-up buildings. Parameters such as connector configurations; variation of the roof stiffness and strength; and coefficient of friction between the panels and slab were studied. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-01-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0073552 |
URI | http://hdl.handle.net/2429/43897 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2013-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2013_spring_tellier_xavier.pdf [ 2.99MB ]
- Metadata
- JSON: 24-1.0073552.json
- JSON-LD: 24-1.0073552-ld.json
- RDF/XML (Pretty): 24-1.0073552-rdf.xml
- RDF/JSON: 24-1.0073552-rdf.json
- Turtle: 24-1.0073552-turtle.txt
- N-Triples: 24-1.0073552-rdf-ntriples.txt
- Original Record: 24-1.0073552-source.json
- Full Text
- 24-1.0073552-fulltext.txt
- Citation
- 24-1.0073552.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0073552/manifest