# INNOVATIVE APPROACHES IN COMPUTATIONAL STRUCTURAL ENGINEERING

EDITED BY : George C. Tsiatas and Vagelis Plevris PUBLISHED IN : Frontiers in Built Environment

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-607-5 DOI 10.3389/978-2-88963-607-5

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# INNOVATIVE APPROACHES IN COMPUTATIONAL STRUCTURAL ENGINEERING

Topic Editors:

George C. Tsiatas, University of Patras, Greece Vagelis Plevris, OsloMet – Oslo Metropolitan University, Norway

Nowadays, numerical computation has become one of the most vigorous tools for scientists, researchers and professional engineers, following the enormous progress made during the last decades in computing technology, in terms of both computer hardware and software development. Although this has led to tremendous achievements in computer-based structural engineering, the increasing necessity of solving complex problems in engineering requires the development of new ideas and innovative methods for providing accurate numerical solutions in affordable computing times.

This collection aims at providing a forum for the presentation and discussion of state-of-the-art innovative developments, concepts, methodologies and approaches in scientific computation applied to structural engineering. It involves a wide coverage of timely issues on computational structural engineering with a broad range of both research and advanced practical applications.

This Research Topic encompasses, but is not restricted to, the following scientific areas: modeling in structural engineering; finite element methods; boundary element methods; static and dynamic analysis of structures; structural stability; structural mechanics; meshless methods; smart structures and systems; fire engineering; blast engineering; structural reliability; structural health monitoring and control; optimization; and composite materials, with application to engineering structures.

Citation: Tsiatas, G. C., Plevris, V., eds. (2020). Innovative Approaches in Computational Structural Engineering. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-607-5

# Table of Contents

Ida Mascolo, Mariano Modano, Ada Amendola and Fernando Fraternali

Christos G. Panagiotopoulos, Vladislav Mantič, Israel G. García and Enrique Graciani

*132 Hybrid Micro-Modeling Approach for the Analysis of the Cyclic Behavior of RC Frames*

Gianni Blasi, Flavia De Luca and Maria Antonietta Aiello

*144 A Comparison of Three Different Methods for the Identification of Hysterically Degrading Structures Using BWBN Model* Ying Zhao, Mohammad Noori, Wael A. Altabey and Taher Awad

*163 Damage Assessment by Numerical Modeling of Sant'Agostino's Sanctuary in Offida During the Central Italy 2016–2017 Seismic Sequence* Ersilia Giordano, Francesco Clementi, Andrea Nespeca and Stefano Lenci

*180 Seismic Assessment of Non-conforming Infilled RC Buildings Using IDA Procedures*

Constantinos Repapis and Christos A. Zeris

*203 FE Modeling of Circular, Elliptic, and Triangular Isolated Slabs With a Continuous Damage Model*

Gelacio Juárez-Luna and Omar Caballero-Garatachea

# Editorial: Innovative Approaches in Computational Structural Engineering

#### George C. Tsiatas 1 and Vagelis Plevris 2 \*

*1 Department of Mathematics, University of Patras, Rio, Greece, 2 Department of Civil Engineering and Energy Technology, OsloMet—Oslo Metropolitan University, Oslo, Norway*

Keywords: structural engineering, innovation, computational, FEM, BEM, meshless methods

#### **Editorial on the Research Topic**

#### **Innovative Approaches in Computational Structural Engineering**

Over the last few decades, tremendous development has been made in the field of computing technology, in terms of both computer hardware and software development. As a result of this progress, numerical computation has now become one of the most effective tools for scientists, researchers, and professional engineers around the world. Although this has led to great achievements in computer-based structural engineering, the increasing necessity to solve complex problems in engineering requires the development of new ideas and innovative methods for providing accurate numerical solutions in affordable computing times (Plevris and Tsiatas, 2018).

The Research Topic "Innovative Approaches in Computational Structural Engineering" aims to provide a forum for the presentation and discussion of state-of-the-art innovative developments, concepts, methodologies, and approaches in scientific computation applied to structural engineering. It involves a wide coverage of timely issues on computational structural engineering with a broad range of research and practical applications.

The collection includes both research and advanced applied topics, with particular emphasis on innovation in computational structural mechanics and engineering. The contributions come from 58 leading researchers and professionals from 10 countries around the world who are actively involved in the field of computational structural engineering. The papers in the collection cover various relevant scientific areas such as modeling in structural engineering; finite element methods; boundary element methods; static and dynamic analysis of structures; performance-based design; structural optimization; meshless methods; modeling of advanced and innovative materials such as shape memory alloys; non-linear structural analysis; system identification methods; structural stability; GPU computing; earthquake engineering; seismic vulnerability; incremental dynamic analysis; structural damage assessment; smart structures and systems; structural reliability; structural health monitoring and control; optimization; and composite materials, with application to engineering structures.

# PAPERS IN THE COLLECTION

In the 1st paper of the collection, Mitropoulou et al. present a methodology to identify the most appropriate damage index, able to provide a reliable description of the structural damage level. A performance-based design framework is formulated based on this damage index, which is to be used as a design tool for achieving safer and more economic designs. For each damage index under consideration, design optimization problems for structural systems are solved by means of a popular metaheuristic search algorithm.

#### Edited and reviewed by:

*Georgios Eleftherios Stavroulakis, Technical University of Crete, Greece*

> \*Correspondence: *Vagelis Plevris<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment*

> Received: *05 March 2020* Accepted: *19 March 2020* Published: *08 April 2020*

#### Citation:

*Tsiatas GC and Plevris V (2020) Editorial: Innovative Approaches in Computational Structural Engineering. Front. Built Environ. 6:39. doi: 10.3389/fbuil.2020.00039*

Next, Charalampakis and Tsiatas examine three uniaxial phenomenological models, i.e., the Graesser-Cozzarelli model, the Wilde-Gardoni-Fujino model, and the Zhang-Zhu model, which are currently used for the description of a Shape Memory Alloy (SMA) behavior. The examination process reveals several limitations and drawbacks, including the large number of parameters and the unclear effect of specific parameters in the overall response. Based on this analysis, powerful metaheuristics are employed for system identification and a new simple rate-independent model is proposed, which addresses all issues in a unified manner producing excellent fit with the experimental data.

The 3rd paper by Yiotis and Katsikadelis deals with the solution to the vibration problem of cylindrical shell panels using the Meshless Analog Equation Method (MAEM). The method is based on the principle of the analog equation, converting the original equations into three substitute ones, two Poisson's and one biharmonic, which are solved using a meshless method. The use of integrated Multiquadric— Radial Basis Functions (MQ-RBFs) to approximate the fictitious sources allows the approximation of the solutions sought by new RBFs, which accurately approximate both the solution and its derivatives. Several shell panels are analyzed using the method, and the numerical results demonstrate its efficiency and accuracy.

Tsalkatidis et al. study the mechanical behavior of bolted hybrid connections that consist of a square hollow steel column (SHS) and a glulam timber beam. A reference model is constructed and verified by comparison to experimental and numerical data from the literature. Several parameters that affect the response of the connection are modified in order to investigate and quantify their effect, resulting in seven different case studies. A proposed optimal configuration of the hybrid connection is investigated and presented.

In the next paper, Tsiatas et al. present a new layered approach to the non-linear analysis of initially straight Euler-Bernoulli rectangular beams by the Boundary Element Method (BEM). The formulation of the problem is based on the displacements, and the equations of equilibrium, derived from the principle of minimum total potential energy, being coupled and non-linear. The beam is discretized in both the longitudinal direction and the cross-sectional plane, and the governing differential equations are solved using the Analog Equation Method (AEM) in conjunction with an iterative numerical process. Various representative examples are examined that not only take geometrical nonlinearity into account but material non-linearity as well. The reliability and effectiveness of the proposed method are validated by comparing the results obtained with those presented in the literature or derived by other Finite Element (FE) models.

Arailopoulos et al. propose an integrated reverse engineering methodology for a large-scale fully operational steam turbine rotor, addressing issues that include developing the CAD and FE model of the structure, as well as the applicability of model updating techniques based on experimental modal analysis procedures. The FE model of the turbine was developed using tetrahedral solid elements resulting in 55 million DOFs. Experimentally identified modal modes and modal frequencies were compared to the FE model predictions, and CMA-ES optimization was used to fine-tune the material parameters, such as modulus of elasticity and density. A simplified FE model was also developed and used for the turbine rotordynamic analysis.

In the 7th paper, Mascolo et al. developed a FE model of the Moderate Rotation Theory (MRT) of laminated composite beams. The proposed model describes laminated composite beams with arbitrary curvature of the beam axis taking into account shear deformation, warping effects, and in-plane and out-of-plane instability. The stability analysis is performed through a path-following procedure and a bordering algorithm. Several numerical results are presented and comparisons are made with classical beam theories and other theories available in the relevant literature.

Jena and Chakraverty investigate the free vibration of the variable cross-section (non-uniform) Single-Layered Graphene Nano-Ribbons (SLGNRs) by using the Differential Quadrature Method (DQM). The Euler–Bernoulli beam theory is considered in conjunction with Eringen's non-local elasticity theory. A convergence study is carried out to illustrate the efficiency of the method and the obtained results are validated with known results.

In the paper by Molina-Moya et al., a parallel iterative solver based on the Generalized Minimum Residual Method (GMRES) with complex-valued coefficients is explored, with applications to BEM. The proposed approach does not require modifications in the main program, which computes the system matrix and the right vector, allowing the use of standard BEM codes, with the solver task being transferred to the Graphic Processing Unit (GPU). The CUDA programming language is used, exploiting the particular architecture of the GPU device for complex-valued systems. The parallelized GMRES solver shows reductions in computing times when compared with its CPU implementation.

In the 10th paper, Panagiotopoulos et al. present a unified methodology to solve problems of frictionless unilateral contact as well as adhesive contact between linear elastic solids. The methodology is based on energetic principles and is applied to a minimization problem of the total potential energy. To solve the quadratic problem, two algorithms are developed, which are both variants of the well-known conjugate gradient algorithm. Both symmetric and nonsymmetric formulations of this matrix are presented and discussed, showing that the non-symmetric formulation provides more accurate results. The framework is found to be very useful in cases where dissipative mechanisms also exist on the boundaries or common interfaces of elastic bodies.

Blasi et al. present a work aimed at developing a hybrid approach to consider the effect of concrete cracking on the hysteretic response of RC frames. The smeared cracking approach is employed, while discrete cracking surfaces are included in the geometrical model where the interface behavior is defined by the combination of contact and cohesive elements. The proposed approach is adopted in ABAQUS to simulate an experimental test on a double cantilever column for the calibration of the numerical model. Numerical and experimental results are compared in terms of hysteretic force-displacement behavior and cumulative dissipated energy, to assess the reliability of the proposed model in simulating the energy dissipation capacity of RC members subjected to lateral cyclic loading.

In the 12th paper by Zhao et al., three different strategies are proposed to identify the general hysteretic behavior of a typical shear structure subjected to external excitations. Different case studies are presented to analyze the dynamic responses of a time-varying shear structural system with the early version of the Bouc-Wen-Baber-Noori (BWBN) hysteresis model. Various computational techniques are utilized. The results show that Intelligent Parameter Varying (IPV) performs superior computational efficiency and system identification accuracy in comparison to Genetic algorithm (GA), and Transitional Markov Chain Monte Carlo (TMCMC) approaches.

The subject of the 13th work of the collection, by Giordano et al., is Sant' Agostino Sanctuary in Offida, central Italy. The village and the sanctuary itself was severely damaged by the seismic events of August 2016. A finite element model using non-linear static analysis (Pushover) is employed to estimate the vulnerability of the sanctuary and its dynamic response, considering the nonlinear behavior of masonry with a homogenized material and a smeared cracking and crushing constitutive law. The structure is shown to be prone to massive damage leading to collapse.

Repapis and Zeris investigate the vulnerability assessment of non-conforming infilled RC structures of the 1960s and 1970s. The seismic performance of typical bare and infilled structures was evaluated. The results from Static Pushover analyses are extended with Incremental Dynamic Analysis predictions using a large number of recorded base excitations from recent destructive earthquakes in Greece and internationally. The analysis shows that fully or partially infilled RC frames can perform well, while frames with an open floor usually have the worst performance due to the formation of an unintentional soft story.

In the 15th paper, Juárez-Luna and Caballero-Garatachea study the nonlinear behavior of RC circular, elliptic, and triangular isolated slabs using computational mechanics. The concrete was modeled with a damage model which includes softening, while the behavior of the reinforcing steel was simulated with a 1D bilinear plasticity model. Validation was done by comparing the computed numerical results with the experimental results reported in the literature. The coefficient method is proposed for its simplicity to calculate design bending moments in slabs with arbitrary geometries.

Soroushian and Farahani apply a technique proposed in 2008 by Soroushian (2008) to the FE analysis of assemblies of beam-columns subjected to continuous static loads. These loads may be computationally expensive especially when they are densely digitized. It is shown that the adapted technique might considerably reduce the computational cost with only a negligible change in accuracy. Furthermore, the good performance of the adapted technique is demonstrated from different points of view and is compared with the performance of the technique in time integration analysis.

In the 17th paper, Marakakis et al. review the current state of the art in shunt piezoelectric systems for noise and vibration control. In fact, resonant shunts have proved to be very efficient and stable for the reduction of vibration on smart piezoelectric structures, such as beams and plates. From the numerous applications which were reviewed it is evident that shunt piezoelectric systems can be very effective for several different purposes such as control of vibrations on hard drives, noise reduction on acoustic applications, improvement of image quality and scan rate of tapping mode Atomic Force Microscopy.

In the last paper of the collection, Papazafeiropoulos et al. propose a novel spectra-matching framework, which employs a linear combination of raw ground motion records to generate artificial acceleration time histories for matching a target spectrum taking not only the acceleration into account, but also the seismic input energy equivalent velocity. The procedure of selection and scaling of the suite of ground motion records to fit a given target spectrum is formulated by means of an optimization problem. The optimization results have shown that there is good agreement between the target and optimized spectra for the different matching scenarios examined, regardless of the nature of the target spectrum.

In closing this editorial, we would like to thank the contributing authors, the reviewers, the staff at Frontiers and all the other people we have worked closely with in trying to make this Research Topic a success. Once again, it has been proven that computational methods are at the forefront of real-world applications. Especially in the structural engineering discipline, the applications are numerous and cover a wide range of topics. Although very powerful capabilities are now available for the simulation and analysis of structures, there are still many exciting research challenges. The field of computational structural engineering will continue to grow and increasingly contribute to technological development in the future, influencing and enriching the engineering profession and our lives in general.

# AUTHOR CONTRIBUTIONS

GT and VP handled manuscripts and edited the Research Topic. Both authors have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

# REFERENCES

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Tsiatas and Plevris. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Damage Index-Based Lower Bound Structural Design

#### Chara Ch. Mitropoulou1 , Giuseppe C. Marano2,3 and Nikos D. Lagaros 4 \*

1 ACE-Hellas, Athens, Greece, 2 College of Civil Engineering, Fuzhou University, Fuzhou, China, 3 Department of Civil Engineering and Architectural Science, Technical University of Bari, Bari, Italy, 4 Veltion Engineering and Applied Sciences Optimization Lab, Department of Structural Engineering, School of Civil Engineering, Institute of Structural Analysis & Antiseismic Research, National Technical University of Athens, Athens, Greece

The scope of the work is to detect the most appropriate damage index, able to provide a reliable description of the structural damage level. This damage index is used to formulate a performance-based design framework to be used as a design tool for achieving safer and more economic designs. This objective is achieved by comparing alternative structural systems that comply with the design demands in the most cost-efficient approach, i.e., those requiring less material volume for concrete and steel reinforcement. In this direction, design optimization problems for structural systems are defined for each damage index under consideration and the structural optimization problems at hand are solved by means of a popular metaheuristic search algorithm.

#### Edited by:

Vagelis Plevris, OsloMet—Oslo Metropolitan University, Norway

#### Reviewed by:

Makoto Ohsaki, Kyoto University, Japan Francesco Tornabene, Università degli Studi di Bologna, Italy Flavia De Luca, University of Bristol, United Kingdom

> \*Correspondence: Nikos D. Lagaros<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment

> Received: 17 April 2018 Accepted: 15 June 2018 Published: 04 July 2018

#### Citation:

Mitropoulou CC, Marano GC and Lagaros ND (2018) Damage Index-Based Lower Bound Structural Design. Front. Built Environ. 4:32. doi: 10.3389/fbuil.2018.00032 Keywords: performance-based design, structural optimization, damage indices, reinforced concrete structures, PSO algorithm

# INTRODUCTION

Performance-based design (PBD) framework was presented throughout the past couple of decades, for designing structural systems exposed to non-deterministic earthquake loading conditions, aiming at increasing structural safety against earthquake hazard. In accordance to performance-based earthquake engineering the structures need to be able resisting various earthquake loading scenarios in a measurable fashion and sustaining potential damages corresponding to desirable performance levels (Vanzi et al., 2015). The contemporary PBD engineering state-of-practice, among others, can be found in US design codes such as FEMA-445 (2006), ASCE-41 (2006), ATC-58 (2009), and FEMA-P-58-1 (2012). Aiming to achieve the PBD engineering goals, these design guidelines suggest that in seismic-prone areas require implementing higher-order analysis methods for performing structural assessment and design of existing or new structural systems, respectively. Structural optimization can provide a highly-effective design framework, yet computationally intensive, exploiting the benefits provided by nonlinear dynamic or static analysis methods. Structural design optimization problems formulated on the basis of PBD frameworks became a topic of growing interest, over the last years, where interesting results have been reached. The progress in structural optimization has facilitated passing from design procedures of trial-and-error basis toward fully automatic ones using advanced design space exploration algorithms (Marano et al., 2009; Quaranta et al., 2014; Fiore et al., 2016a; Greco et al., 2016). Partly, this transition should be credited to metaheuristic search algorithms that last decades developed rapidly. These search algorithms are capable of dealing with optimization problems of high complexity. Furthermore, structural optimization matured from simple academic problems to become the core of contemporary design in case of extremely complicated structural systems (Marano and Greco, 2006; Marano et al., 2006, 2007, 2008; Lagaros, 2014a, 2018; Greco et al., 2015).

So far, many researchers have integrated the PBD concept into design optimization procedures for achieving better designs. Indicatively, Ganzerli et al. (2000) suggested a PBD optimization practice for reinforced concrete (RC) framed structures relying on mathematical programming search algorithms. Esteva et al. (2002) recommended a formulation for life-cycle cost (LCC) aiming to define optimized mechanical properties of structures subjected to earthquake risk. Li and Cheng (2003) presented a design concept based on damage reduction and was implemented into structural optimization problems and exhibited that this idea leads to improved designs with respect to structural performance in terms of both maximum interstory drift and life-cycle cost performance. Chan and Zou (2004) presented an improved approach for RC building structures subject to seismic loading aiming to improve elastic and inelastic interstory drifts into structural optimization problems. Liu et al. (2005) presented a design optimization approach for multipleobjectives using Genetic Algorithm (GA) for the case of steel framed structures taking into account various criteria; for two performance levels the resulting maximum interstory drift, material weight, and design complexity ones. Fragiadakis et al. (2006) presented PBD optimization practice for steel momentresisting frames (MRFs) based on FEMA-350 (2000) probabilistic framework. An overview for the PBD framework state-of-theart was presented by Foley et al. (2007) and implemented such a design framework into multiple-objective structural optimization problem formulations for single and multi-story steel framed structures having fully and partially restrained connections. European seismic code recommendations were assessed by Lagaros and Papadrakakis (2007), when adopted for designing 3D RC building structures, within a PBD multiobjective optimization concept. Aiming to deal with complex optimization problems, confidence levels were integrated with the objective criterion together with the material cost for various optimized design scenarios, the GA search algorithm was adopted by Rojas et al. (2007). Lagaros et al. (2010) presented an automated PBD approach that was implemented for designing RC building structures, examining the effect of infill walls, in this direction PBD approach alternatives were tested relying on non-linear static or dynamic analyses. In the work by Lagaros and Fragiadakis (2011) a comparative study is performed over the state-of-practice on nonlinear static analysis approaches, recommended by the European and US code provisions, were implemented within PBD optimization formulations; while aiming to examine the parameters that affect LCC estimation of structures, multi-objective LCC optimization problems were formulated by Mitropoulou et al. (2011).

The concept of defining the condition of structural loss using an engineering compatible quantity taking values on specific scale by means of a damage index (DI) is appealing due to its ease of implementation. So far several researchers have proposed several DIs for the case of RC or steel structural systems (i.e., Ghobarah, 2004; Estekanchi and Arjomandi, 2007). This work aims at exploring and assessing alternative designs, exhibiting a required performance monitored by a damage index, that comply with design codes according to the utmost economical way, i.e., having lower material volume requirements for concrete and steel reinforcement. In this direction, performance-based structural design optimization problems are formulated for various damage indices and the designs obtained are subsequently assessed

under a set of seismic records properly selected. The main objective is to detect the most appropriate damage index, able to provide a reliable description of the structural damage level. This damage index is used to formulate a performancebased design framework to be used as a design tool for achieving safer and more economic designs. The resulting structural optimization problems are solved using a well-known metaheuristic optimization algorithm, namely the particle swarm optimization method, while the limit states implemented for the formulation of the performance-based design procedures, which describe the status of structural integrity, are based on calibrated values obtained from a former study of the authors (Mitropoulou et al., 2014).

# LOCAL AND GLOBAL DAMAGE INDICES

Damage indices can be classified into two groups (Grigoriu, 1987): (a) DIs based on strength; and (b) those based on response. Strength-based DIs do not require response analyses (Shiga et al., 1982; Khashaee, 2005). However, they need to be calibrated with reference to observed damage based on rather large experimental databases. The performance of structures against seismic loading is generally correlated to the ability to sustain inelastic deformations (Resta et al., 2013), while experimental studies have shown that indices of structural performance (like ductility as well as others), which rely only on the fatigue concept of low-cycle, do not appear to offer an acceptable indicator for structural loss due to seismic loads (Banon and Veneziano, 1982). These observations are in agreement with the notion that brittle systems' failure is developed due to significant deformation demands, while ideal ductile structures' failure begun when inelastic deformations are observed repeatedly (Fiore et al., 2016b). Damage indices used to characterize the performance of structural systems that are not ideal ductile or brittle, need to consider for losses due to repeated and increased inelastic deformations (Park and Ang, 1985; Greco et al., 2017). Thus, more wide-ranging applicable

and reliable indices are required to characterize structural performance.

According to Ghobarah et al. (1999) the DIs based on response are classified in three sets based on what DI accounts for: (a) maximum deformation (Mahin and Bertero, 1974; Saiidi and Sozen, 1981; Toussi and Yao, 1982; Roufaiel and Meyer, 1987; Powell and Allahabadi, 1988); (b) cumulative damage (Banon and Veneziano, 1982; Allahabadi and Powell, 1988); and (c) combination of cumulative damage and maximum deformation. In this work DIs accounting for both cumulative damage and maximum deformation are implemented. In particular, (i) the local DI by Park and Ang (1985), (ii) its alternative presented by Kunnath et al. (1992), (iii) local DI by Chung et al. (1987, 1989), (iv) together with the final and maximum softening DIs presented by DiPasquale and Çakmak (1987, 1988) are evaluated in the current study. For their implementation the local DIs are transformed into global ones by combining the local values of the DIs of the structural elements using weight coefficients.

Several damage indices account for the local damage induced to individual structural elements. Monitoring of local damage indices helps in the identification of the weak or vulnerable members that need to be rehabilitated. Nevertheless, it is not an easy task to obtain clear view for the structure's response provided the input earthquake record out of a list of members DIs. Currently, safety and residual strength of damaged structural systems rely on single or global DIs. The global ones are usually implemented for reliability studies, post-earthquake evaluation of structures, and performance-based design applications. Global DIs could be defined by combining local ones; the weighting scheme is the simplest approach for combining local damage indices (Kunnath et al., 1992). The weighting factors can predict the relative importance and/or the replacement cost of a structural element or part of the structure in preserving safety of the entire structure. However, employing weighted average procedures in order to define global DIs does not justify appropriately the concentration locally of the damages, does not distinguish between different structural elements, and might lead to misrepresentative conclusions. There are cases that for a few structural elements of the structural system to have undergone severe damage deprived of this being reflected in by a DI.

For the reasons described previously DIs expressed using weighting factors have not been considered; however, for purposes of completeness such DIs are briefly outlined below. Direct implementation of the damage model to the overall structural system, a single story frame or a structural member requires determining the corresponding overall structural system, story or member ultimate deformations. Nonlinear behavior is limited to plastic regions close to the structural joints; therefore, it is difficult to establish the relation between member, story or top story structural system deformations, with the rotations of the local plastic zone. Two other global and story damage indices were proposed (Kunnath et al., 1992), which are calculated using weighting factors, these factors rely on the hysteretic energy dissipated at story and component levels,

respectively:

$$DI\_{\text{story}} = \sum\_{i=1}^{n\_{cl}} (\lambda\_i)\_{\text{element}} (DI\_i)\_{\text{element}} \tag{1}$$

$$(\lambda\_i)\_{element} = \frac{E\_i}{\sum\_{\substack{n\_{cl} \\ i=1}}^{n\_{cl}} E\_i} = 1, \dots, n\_{cl} \text{(over the elements)} \quad \text{(2)}$$

and

$$DI\_{overall} = \sum\_{i=1}^{n\_{\rm sf}} \langle \lambda\_i \rangle\_{story} \langle DI\_i \rangle\_{stony} \tag{3}$$

$$(\lambda\_i)\_{\text{story}} = \frac{E\_i}{\sum\_{\substack{n\_{\text{st}}\\i=1}}^n} i = 1, \dots, n\_{\text{st}} \\ (over \text{ the stories}) \qquad \text{(4)}$$

λi denotes the energy weight factor for story i; Ei represents the overall energy absorbed by the ith story or component, while nst and nel denote the total number of stories and components/elements, respectively.

#### LOWER BOUND DESIGN OF STRUCTURES

This work concerns the production of lower-bound optimized designs, which comply with design requirements by a costefficient means (i.e., those requiring less material volume for concrete and steel reinforcement) and subsequently to be compared. These designs are then used as the basis of comparison for identifying the best design procedure. Therefore, structural optimization problems are formulated and optimized designs obtained are subsequently evaluated.

#### Design of Reinforce Concrete Structures

Aiming to evaluate the performance of different designs of RC framed structures; structural optimization problems were formulated as below:

$$\begin{array}{ll}\text{minimize} & \quad \text{CIN}(t, \mathbf{s})\\ \text{subject to} & \quad \mathbf{g}\_{j}^{\text{SERV}}(\mathbf{s}) \le 0 \mathbf{j} = 1, \ldots, k \end{array}$$

$$\mathbf{g}\_{j}^{\text{ULT}}(\mathbf{s}) \le 0\\j = k + 1, \ldots, m \tag{5}$$

The design variables comprising of the dimensions of crosssectional for groups of beams and columns considered for formulating the optimization problem are defined by vector **s**, F represents the feasible part of the design space where all serviceability and ultimate limit state design requirements implemented as constraint functions (i.e., the series of constraint functions g SERV j (s) and gULT j (s)) are satisfied:

$$F = \left\{ \mathbf{s} \in R^D \, \middle| \, \mathbf{g}\_j^{\text{SERV}}(\mathbf{s}) \le 0, j = 1, \dots, k \text{ and } \mathbf{g}\_j^{\text{ULT}}(\mathbf{s}) \le 0 \, j = k \right\} \tag{6}$$

where R D represents the discrete design set where design vectors **s** take values. The structural materials cost CIN of the design is the objective function considered. Aiming to deal with the optimization problem the well-known metaheuristic search algorithm called as particle swarm optimization (PSO) method is implemented.

#### Swarm Intelligence

According to the optimization algorithm particle swarm (Kennedy and Eberhart, 1995), various design vectors collaborate. Every vector of unknowns is labeled as "particle" characterized by its velocity and position, both defined in the D-dimensional design domain, while a sum of particles constitutes the so termed "swarm." During the search procedure where an optimized design vector is on the hunt, a particle "flies" in the design space of the problem. During the search procedure, velocity and position vectors are adjusted for every particle based on personal "experience" as well as that of the others (neighboring particles). Memorizing and tracking the best positions encountered construct particles' experience. PSO algorithm relies on joint search locally (self-experience) with global one (neighboring experience), trying to control exploitation with exploration. Every particle maintains its two basic characteristic vectors, i.e., velocity and position (or location), in the D-dimensional domain which are iterated as below:

$$\mathbf{v}^{j}(t+1) = \boldsymbol{\omega}\mathbf{v}^{j}(t) + c\_{1}\mathbf{r}\_{1} \cdot \left(\mathbf{s}^{Pbj} - \mathbf{s}^{j}(t)\right) + c\_{2}\mathbf{r}\_{2} \cdot \left(\mathbf{s}^{Gb} - \mathbf{s}^{j}(t)\right) \tag{7}$$

$$s^j(t+1) = s^j(t) + \nu^j(t+1) \tag{8}$$

**v** j (t) refers to the jth particle's velocity at time t, **s**j(t) denotes the jth particle's position vector at time t, the self-best location of jth particle is labeled as **s** Pb,j , and the best position obtained

globally is denoted with the vector **s** Gb. Coefficients c1 and c2 denoted as acceleration ones, refer to the level of confidence with the best vector achieved by a specific particle (c1 is labeled as cognitive coefficient) and by entire swarm (c2 is denoted as social coefficient), respectively. Vectors r1 and r2 are composed by arbitrary elements distributed uniformly in [0,1]. Flowchart of **Figure 1** describes the basic steps of PSO algorithm, while particle's motion indicatively for the case of the two-dimensional design space is depicted in **Figure 2**. The lower left dotted circle of **Figure 2** represents present position vector **s** j (t) at time t, and the upper right dotted bold circle denotes new location vector **s** j (t+1) at time t+1. The means that the particle's transition into the D-dimensional search domain is affected by: (i) velocity vector **v** j (t); (ii) self-best found location vector (**s** Pb,j ); and (iii) swarm's global-best location achieved so far (**s** Gb) is presented in **Figure 2**.

# DAMAGE INDEX-BASED DESIGN

In this work the performance degree are denoted in terms of structural deformation and cumulative energy, using different damage indices apart from maximum interstory drift that has been exclusively used in the past in PBD. More specifically Park and Ang DI and its reformed version as presented by Kunnath et al. (1992), as well as the Chung et al. DI, together with the final maximum and softening DIs suggested by DiPasquale and Çakmak (1988), have been considered for the definition of the performance levels. These five DIs account for cumulative damage and maximum deformation. These damage indices were calibrated in a recent study by the authors (Mitropoulou et al., 2014) relying on the crack openings' width. In particular, incremental dynamic analysis was employed for developing an archive of the DIs' values using a coarse modeling of structures employing column-beam 1D elements.

Six optimization problems are pursued within, following the PBD framework based on (i) interstory drift (DIθ ), (ii) Park and Ang (DIPA); (iii) its alternative presented by Kunnath et al. (DIKRL); (iv) Chung et al. (DICMS); (v) maximum softening (DIms) as well as (vi) final softening (DIfs). Moreover, DIs statistical features were calculated implementing horizontal statistics combined with the method of maximum likelihood function together with a properly selected search technique. The values provided in **Table 1** (Mitropoulou et al., 2014) represent calibrated damage states and they are used for defining the performance objective introduced in the proposed design framework.

# Performance-Based Design

Several distinctive features for PBD seismic design procedure with respect to design procedures imposed by the prescriptive codes: (i) Permits structural engineers to select both appropriate seismic hazard level and the analogous performance of the structure, (ii) series of performance objectives are used for designing a structural system. PBD design procedure represents a design procedure relying on displacements where design measures and capacity demands are formulated considering displacements rather than using forces (Sullivan et al., 2003).

PBD implies that structural elements are selected, as well as the assessment; construction and maintenance studies are operated on the construction aiming to meet the purposes imposed by owners/users as well as society (Krawinkler and Miranda, 2004). When designing against seismic hazard, the aim is to construct structural systems having predictable performance capable to withstand seismic loading with measurable metrics. Hence, contemporary conceptual approach for structural design is that structural systems need to achieve performance objectives for several seismic levels varying in the rage of seismic event having low intensity and return period, to more damaging earthquake incidents having larger return periods. The contemporary stateof-practice for PBD structural engineering is offered by the various design codes (ASCE-41, 2006; FEMA-445, 2006; ATC-58, 2009; FEMA-P-58-1, 2012), that theoretically do not vary besides adopt techniques that represent the first significant difference with respect to prescriptive-based guidelines. Several of the contemporary design guidelines applied to design studies for new build structures can be considered as merely PBD ones, as these design guidelines combine all measures into a single level of performance, frequently the collapse prevention or life safety level. The definition of performance objectives represents the most important ingredient of a PBD seismic procedure.

#### DI-Based Design Concept for Safe and Economic Seismic Resistant Structures

The main scope of this step is to incorporate the DIs calibrated in a recent work by the authors into a performance-based design framework and to propose an innovative design concept leading to safe and economic earthquake resistant RC structures. For this purpose, the DIs selected in Mitropoulou et al. (2014) will be used to define different DI-based designs into structural optimization problems. The formulation of the DI-based concept can be stated as follows:

$$DI\_{HL1}^{i}(\mathbf{s}) \le DI\_{HL1}^{\text{allow,i}}$$

$$DI\_{HL2}^{i}(\mathbf{s}) \le DI\_{HL2}^{\text{allow,i}}$$

$$\vdots$$

$$DI\_{HLn}^{i}(\mathbf{s}) \le DI\_{HLn}^{\text{allow,i}}\tag{9}$$

which represent the performance-based design constraints [denoted as g ULT(**s**) in Equation 5]. The DI should be less or equal to the allowable upper bounds for a number of hazard levels (HL1, HL2,. . . ,HLn) while i denotes the DI considered in the formulation (i=1,2,. . . ,m). The boundaries of the feasible part of the design space are defined using the calibrated values of the DIs defined in Mitropoulou et al. (2014) which are given in **Table 2**. Therefore, in order to identify the DI that represents better the structural behavior, various minimum initial material cost DIbased design optimization problems are formulated according to Equation (5), where the PBD checks are implemented according to Equation (9).

#### NUMERICAL TESTS

The objective of this work is to incorporate several damage indices examined in former study by the authors (Mitropoulou et al., 2014) into a performance-based design framework and to identify the best choice to be used for designing reinforced concrete structures.

#### Description of the Test Examples

The two 2D RC moment resisting framed (MRF) building structures shown in **Figure 3** have been taken into account aiming to study the effect of incorporating DI within the design framework of RC building structures. The foremost test example refers to a three-story concrete building while the second one is a six-story one. Class of concrete C20/25 and of steel S500 are assumed. The thickness of the slabs was taken equal to 15 cm for both test examples and their contribution to the moment of inertia of the beams through TABLE 3 | Three-story test example: Concrete volume and steel weight corresponding to the optimized designs based on six DIs.

TABLE 4 | Six-story test example: Concrete volume and steel weight corresponding to the optimized designs based on six DIs.

effective flange width is taken into account. Due to floor finishing-partitions permanent loading of 2 kN/m2 and moving one equal to 1.5 kN/m2 , are applied further to the self-weight of beams and slabs, where an effective zone 10 × 15 m2 is considered for each story. The nominal imposed and permanent loads are combined using the load factors 1.50 and 1.35, respectively.

Columns and beams were simulated adopting the inelastic force-based fiber finite element. The simulation with this type of elements exhibits better accuracy than the plastic hinge formulation (Lagaros, 2014b). The structural analyses required in the framework the study were performed using OpenSEES (McKenna and Fenves, 2009) software. For steel reinforcement pure kinematic hardening bilinear material model was implemented, considering also geometric nonlinearities explicitly. The extended Kent-Park model as modified by Scott et al. (1982) was used for simulating concrete. Despite its relatively simple formulation, the specific model offers for acceptable forecasts of the needs for flexure-dominated RC members. The cyclic inelastic behavior of bars reinforcement was implemented with Menegotto-Pinto model (Menegotto and Pinto, 1973).

#### Analysis Procedure

During seismic structural design and/or assessment a wideranging earthquake events and multiple structural response levels are required for accounting the uncertainties that earthquake hazard introduces into performance-based earthquake engineering (PBEE) assessment or design. Methods implemented for PBEE assessment by means of nonlinear dynamic analyses are classified into single and multiple hazard level ones. Incremental dynamic analysis (IDA) and multi-stripe dynamic analysis (MSDA), their single and multicomponent (Lagaros, 2010) variants, represent the most appropriate ones. The numerical study performed consists of two stages, the optimization and the structural assessment ones, for both multi-stripe dynamic analysis (MSDA) was implemented for calculating the maximum inter-story drift.

Similar to IDA the purpose of MSDA studies is to form the relationship among earthquake intensity level and the correlative maximum structural performance. Earthquake ranking and performance of the structural system are defined through intensity measures (IM) and engineering demand parameters (EDP), correspondingly. MSDA implementation involves multiple inelastic dynamic analyses (stripes) executed in various spectral acceleration stages (see **Figure 4** for the case of the final softening and Park and Ang designs). In every stripe various inelastic dynamic analyses are executed for the seismic records considered, all are scaled to specific spectral acceleration. The group of seismic records adopted for deriving each stripe should preferably be representative of the seismic hazard at specific spectral acceleration; nevertheless, although not necessarily always justified (e.g., Jalayer and Cornell, 2009), it is common to implement the same group of records for all spectral acceleration levels.

#### Formulation of the Optimization Problems

The problem is expressed mathematically as single-objective minimization problem and is shown below:

FIGURE 5 | Three-story test example—performances of the optimized designs (i.e., DIθmax, DIms, DIfs, DIPA, DIKRL, and DICMS) for three hazard levels compared with respect (A) maximum drift, (B) maximum softening, (C) final softening, (D) Park and Ang damage index, (E) Kunnath et al. damage index, and (F) Chung et al. damage index.

$$\begin{array}{ll}\min\_{\mathbf{s}\in F} & C(\mathbf{s})\\\text{where} & C(\mathbf{s}) = C\_b(\mathbf{s}) + C\_{sl}(\mathbf{s}) + C\_{cl}(\mathbf{s}) + C\_{ns}(\mathbf{s})\\\text{subject to } & DI\_{HL1}^i(\mathbf{s}) \le DI\_{HL1}^{\text{allow},i} \\\ & DI\_{HL2}^i(\mathbf{s}) \le DI\_{HL2}^{\text{allow},i} \\\ & \vdots \\\ & DI\_{HLn}^i(\mathbf{s}) \le DI\_{HLn}^{\text{allow},i} \end{array} \tag{10}$$

subject to the constraints of Equation (9) (Quaranta et al., 2014; Fiore et al., 2016a). The total initial structural material cost C(**s**) is the objective criterion considered, and Ccl(**s**), Cb (**s**), Csl(**s**), and Cns(**s**) are the total material cost of columns, slabs, beams ,and nonstructural members, correspondingly. The ith DI for several hazard levels (HL1, HL2, . . . , HLn) should remain below the allowable upper bounds (DIi HLj(s) ≤ DIallow,i HLj ). The limits of the feasible region are given in **Table 2** (Mitropoulou et al., 2014). The PSO algorithm is adopted for dealing with the resulting structural design optimization problems. The boundaries for the maximum interstory drift were obtained by Ghobarah (2004).

The optimized design is computed by applying the PSO method with the following characteristic parameters for both test examples as suggested by the parametric study of Pedersen (2010): Number of particles NP is equal to 100, inertia weight w is equal to −0.6, cognitive and social parameters c1 is equal to −0.65 and c2 is equal to 2.65, respectively. Rectangular crosssectional shape is adopted for both columns and beams, which are separated into different sets. The beams/columns dimensionality together with the longitudinal reinforcing bars represent the design parameters, which are allocated to every set of the beams/columns. The structural members (columns and beams) are divided in two sets for the three-story and three sets for the six-story test examples, resulting into five and seven design variables, respectively.

The optimized designs obtained are labeled as DDI=i associated with the damage index used (i–vi). The optimized

damage index.

designs obtained for different DI considerations with reference to steel and concrete material demands are shown in **Tables 3**, **4**, respectively. According to **Tables 3**, **4** it looks that DDI=ms, DDI=fs, and DDI=CMS are leading to smaller columns for the three-storys but it does not apply to the six-story where DDI=θ is the one leading to smaller columns. With respect to beams, it looks that DDI=θ and DDI=CMS are leading to smaller beams for the three-stories and DDI=θ and DDI=KRL are leading ones for the six-story. Compared to design DDI=θ, it can be seen that, the initial material construction cost of the other five designs is increased by 23–145% for the three-story and by 26–75% for the six-story test example.

# Assessment of the Optimized Designs

Aiming to evaluate the capacity of the optimized designs achieved by means of the six design PBD procedures described previously, a multi-stripe analyses is performed over a properly selected bin 100 natural records from the list of records provided in Mitropoulou et al. (2015). Seismic records were chosen from PEER (2010) database in agreement with the subsequent characteristics: (i) Records chosen in particular geographical zone (latitude 32–41◦ , longitude −124◦ to −115◦ ). (ii) Magnitude (M) to be greater or equal to 5. (iii) Distance from the epicenter (R) is <150 km. **Figures 5**, **6** depict the median values of maximum interstory drifts, maximum softening, final softening, Park and Ang damage index, Kunnath et al. damage index, and Chung et al. damage index. These quantities are obtained for three hazard levels namely 50/50, 10/50, and 2/50 (corresponding to probabilities of exceedance equal to 50, 10, and 2% in 50 years, respectively) and for the optimized designs achieved according to the design frameworks discussed above. Indicatively the multi-stripe analysis results for the three-story with reference to the maximum drift for final softening and Park and Ang designs that were used to derive the corresponding for these two designs median curves of **Figure 5A** are provided in **Figure 4**. In **Figure 5A** is represented the structural performance of the three-stories optimized designs (i.e., DIθmax, DIms, DIfs, DIPA, DIKRL, and DICMS) for the three hazard levels based only on maximum drift; the different lines are representing the different design.

The measure of the performance of one DI over the others is the structural performance achieved by each optimized design compared to the other ones, i.e., the design that has on average the better structural performance for all DIs. For the threestory test model it is the design obtained according to Park and Ang damage index that shows an overall better performance with reference to the six engineering demand parameters. On the other hand, for the six-story test model it is the design achieved according to Chung et al. damage index that exhibits the best performance, with reference to the six engineering demand parameter. Furthermore, for both test examples the maximum interstory drift appear to increase smoothly based on the abscissa escalation for all optimized designs in consideration of the hazard levels, while the structural behavior of the six designs is nonlinear with reference to the other damage indices. In general, it can be said that, for the two test examples considered, it is the Chung et al. damage index that shows an overall good performance with reference to the performance criteria considered. This is because for both test examples the designs obtained through the design formulation with the Chung et al. DI corresponds to the lower DI values, with reference to all six damage indices considered.

#### CONCLUSIONS

In this work we aim to incorporate a number of damage indices previously assessed by Mitropoulou et al. (2014) into a performance-based design framework and to identify the most suitable one for designing reinforced concrete

#### REFERENCES

framed structures. In this direction, structural optimization problems are formulated with performance-based design criteria corresponding to six damage indices for two reinforced concrete buildings.

An investigation was performed on the influence of these six damage indices incorporated into a performance-based design procedure. This was achieved by means of lower bound performance-based design. The ultimate objective was to compare lower-bound designs that satisfying in the most costefficient way the design code requirements, i.e., those requiring less material volume for concrete and steel reinforcement. The optimized designs obtained were subsequently assessed on the basis of their performance on a number of properly selected earthquakes with increasing intensity.

From the results obtained it was observed that, compared to design DDI=θ based on interstory drift, the initial construction cost of the other five designs is enlarged by 23–145% for the three-story test example and by 26–75% for the six-story one. It can also be concluded that for the two test examples considered it is the Chung et al. damage index that shows an overall good performance with reference to the performance criteria considered.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

### ACKNOWLEDGMENTS

This research was supported by the OptArch project: Optimization Driven Architectural Design of Structures (No: 689983) belonging to the Marie Skłodowska-Curie Actions (MSCA) Research and Innovation Staff Exchange (RISE) H2020-MSCA-RISE-2015.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Mitropoulou, Marano and Lagaros. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Simple Rate-Independent Uniaxial Shape Memory Alloy (SMA) Model

Aristotelis E. Charalampakis 1 \* and George C. Tsiatas 2

*1 School of Civil Engineering, National Technical University of Athens, Athens, Greece, 2 Department of Mathematics, University of Patras, Rio, Greece*

In this paper, three specific uniaxial phenomenological models commonly used for the description of a Shape Memory Alloy (SMA) behavior are examined in detail. In particular, the models examined are the Graesser-Cozzarelli model, the Wilde-Gardoni-Fujino model, and the Zhang-Zhu model. The pertinent model parameters are examined with respect to their physical representation, if any. Based on this analysis, a new simple rate-independent model is proposed which addresses all issues in a unified manner. Finally, powerful metaheuristics are employed for system identification, producing excellent fit with experimental data while revealing valuable information regarding the relative sensitivity of the proposed model parameters.

Keywords: Shape Memory Alloy (SMA), identification, hysteresis, rate-independent model, phenomenological

#### Edited by:

*Nikos D. Lagaros, National Technical University of Athens, Greece*

#### Reviewed by:

*Dimitrios Giagopoulos, University of Western Macedonia, Greece Savvas Triantafyllou, University of Nottingham, United Kingdom*

> \*Correspondence: *Aristotelis E. Charalampakis<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment*

> Received: *11 June 2018* Accepted: *30 July 2018* Published: *28 August 2018*

#### Citation:

*Charalampakis AE and Tsiatas GC (2018) A Simple Rate-Independent Uniaxial Shape Memory Alloy (SMA) Model. Front. Built Environ. 4:46. doi: 10.3389/fbuil.2018.00046* INTRODUCTION

model

Throughout history, humanity has sought shelter from natural phenomena. Natural shelters, such as caves, were abandoned for artificial structures made of primitive materials, i.e., wood, stone, and brick. The building blocks were held in place by either sole gravity or by primitive paste. In the classical era, the Romans used opus caementicium, a pozzolanic concrete-like material, to build arches, vaults, and domes. Following a prolonged period of stagnation, new materials were introduced around 1800–1850 (i.e., cast iron and, later on, wrought iron and structural steel, as well as Portland cement), causing a revolution in the construction industry. Nowadays, the extensive use of high-strength steel and concrete, as well as the combination of them (as in composite structures), has led to an explosion in the capacity, performance, and size of structures that are feasible. Nevertheless, there is a consensus that, more or less, these materials have been fully exploited and groundbreaking development is unlikely to be observed.

Seeking the next revolution in the construction industry, researchers drew their attention to the so-called smart materials. These materials exhibit extraordinary properties, ranging from piezoelectricity and pH-sensitivity to magnetostriction and self-healing. A popular class of these materials, commonly known as Shape Memory Alloys (SMAs), exhibit physical and mechanical characteristics that allow their integration into structures. SMAs are capable of sustaining large inelastic strains that can be recovered by heating or unloading, depending on prior loading history. The origin of this unusual behavior is the ability of SMAs to undergo a first-order solid-solid diffusionless, and reversible phase change called martensitic transformation between a parent phase called austenite (A), stable at high temperature and low stress, and a product phase called martensite (M), metastable at low temperature and high stress (Olson and Cohen, 1982; Cisse et al., 2016). The martensitic transformation is at the origin of the two main effects observed in SMAs, namely the shape-memory effect and superelasticity.

Various innovative systems and devices, mainly using NiTi and Cu-based SMAs, have been developed for seismic energy absorption, damping control, structural retrofit. Several prototypes of SMA braces for the seismic protection of structures have been designed, numerically assessed and experimentally tested (Clark et al., 1995; Auricchio et al., 2006; Andrawes and DesRoches, 2007; Zhu and Zhang, 2008; Asgarian and Moradi, 2011). Hybrid devices, comprising of a combination of SMA and other components such as precompressed springs (Ma and Cho, 2008; Miller et al., 2012), or struts and steel tubes (Walter Yang et al., 2010) have been presented. Also, SMA-based isolation devices have been developed finding various applications on buildings or bridge structures (Krumme et al., 1995; Wilde et al., 2000). These devices exploit the hysteretic nature of SMA materials aiming to absorb structural vibrating energy and mitigate seismic effects. Gur and Mishra (2013) proposed a combination of steel-Teflon surface isolation device with superelastic SMA restrainers for energy dissipation and simultaneous minimization of isolator displacements and superstructure accelerations. Ozbulut and Hurlebaus (2010) introduced a sliding-type isolation system, optimally designed with the help of a multi-objective genetic algorithm, which consists of a steel Teflon sliding bearing and a NiTi SMA device, providing vertical load capacity, damping, and self-centering at the same time. Other applications of SMA materials in structural engineering include damping devices for bridge structures (Fujino et al., 1993; Ben Mekki and Auricchio, 2011; Torra et al., 2013), SMA-based structural connections in steel structures (Leon et al., 2001; Speicher et al., 2011; Fang et al., 2014) and reinforced concrete frames (Youssef et al., 2008; Muntasir Billah and Shahria Alam, 2012), structural retrofit of buildings and bridges (DesRoches and Delemont, 2002; Andrawes et al., 2010), use of SMA as reinforcement material in concrete structures (Deng et al., 2006; Zafar and Andrawes, 2013) as well as selfrehabilitation of structural elements (Li et al., 2006, 2008).

Obviously, proper modeling of the extraordinary behavior of SMAs is important and, thus, several models have been proposed in the literature. These can be broadly categorized into microscopic thermodynamic models, based on the Ginzburg-Landau theory or molecular dynamics; micro-macro models, based on micromechanics, micro-planes or micro-spheres; and macroscopic models, based on the theory of plasticity, thermodynamic potentials, finite strains or statistical physics (Cisse et al., 2016). For applications in structural engineering, however, uniaxial phenomenological models are of special interest as integration with existing reliable FEM codes, e.g., OpenSees (Mazzoni et al., 2006), is straightforward. Falling under this category, the Graesser-Cozzarelli model (Graesser and Cozzarelli, 1991), the Wilde-Gardoni-Fujino model (Wilde et al., 2000), and the Zhang-Zhu model (Zhang and Zhu, 2007) are examined herein in detail. Their parameters are examined concerning their physical representation, if any. This process reveals several limitations and drawbacks, most prominent of which is the large number of parameters and the unclear effect of specific parameters in the overall response. These are addressed in a unified manner in a new rate-independent model proposed herein. Finally, metaheuristic algorithms are employed for the identification of unknown parameter values. This process reveals valuable information regarding the relative sensitivity of the parameters as well as the mathematical consistency of the model.

#### BASIC SMA BEHAVIOR

The primary forms of SMA behavior which are pertinent to the applications examined in this study are shown in **Figure 1**. In particular, **Figure 1A** shows the hysteretic loop observed at ambient temperatures T < Mf , where Mf is the temperature at which the microstructure of the material is fully martensitic. Although the loop strongly resembles the one observed in most conventional steels, the hysteretic mechanism is quite different. In steels, hysteresis in cyclic loading is due to dislocation glide but in SMAs it is due to twinning deformation of martensite that occurs by rotation, growth, and shrinkage of individual variants of martensite (Graesser and Cozzarelli, 1991). By application of sufficient strain, only one martensitic variant remains. This variant reverts to the original parent crystal orientation upon the application of heat, thus the term shape-memory effect.

**Figure 1B** shows the hysteretic loop observed at ambient temperatures T > Af , where Af is the temperature at which the microstructure of the material is fully austenitic. This behavior exhibits two very important properties, i.e., energy dissipation and zero residual strain upon unloading, and thus it is termed superelastic. According to this, the stable austenitic

microstructure is loaded elastically up to a threshold stress level, at which a transformation from austenite to martensite is initiated. The transformation is accompanied by the reduced modulus, as compared to the initial elastic loading, which strongly resembles the plastic yielding of steels. Gradually, the volume ratio of martensite is increased within the microstructure until it becomes dominant. Further loading of the now fully martensitic microstructure leads to elastic loading with a modulus smaller than the one of the elastic austenite, yet significantly higher than the one during transformation. Upon unloading, the inverse process is observed but at a significantly lower stress plateau. Ideally, the material returns to its exact original form (Graesser and Cozzarelli, 1991).

# THE GRAESSER-COZZARELLI MODEL

By modifying Ozdemir's model (Ozdemir, 1973), which is a particular case of the well-known Bouc-Wen hysteretic model (Bouc, 1967; Wen, 1976), Graesser and Cozzarelli (1991) proposed a one-dimensional phenomenological law to describe both twinning hysteresis and superelasticity, which is described by the following equations:

$$\dot{\sigma} = E\left[\dot{\varepsilon} - |\dot{\varepsilon}| \left(\frac{\sigma - \beta}{Y}\right)^n\right],\tag{1}$$

$$\beta = E \alpha \left[ \varepsilon\_{in} + f\_T |\varepsilon|^c \text{erf} \left( \bar{a} \varepsilon \right) \mu \left( -\varepsilon \dot{\varepsilon} \right) \right], \tag{2}$$

where, (˙) = ordinary time derivative, σ = stress, ε = strain, β = backstress, E = initial elastic modulus, Y = yield stress, α = Ey/ E − Ey = parameter controlling the post-elastic slope of the curve (Ey = post-elastic modulus), n = parameter controlling the abruptness of transition from the elastic to postelastic branch, εin = ε − σ/E = inelastic strain, fT, a¯, c are model parameters, erf (·) is the error function defined by:

$$\operatorname{erf}\left(\mathbf{x}\right) = \frac{2}{\sqrt{\pi}} \int\_0^\mathbf{x} e^{-t^2} dt,\tag{3}$$

and u (·) is the Heaviside step function defined by:

$$\mu(\mathbf{x}) = \begin{cases} 1, & \mathbf{x} \ge \mathbf{0} \\ 0, & \mathbf{x} < \mathbf{0} \end{cases}. \tag{4}$$

Note that parameter a¯ is used in Equation (2) to avoid confusion between parameters a and α in the original paper. Also, in Equation (1) the term (σ − β) may be negative, which is problematic when n takes real rather than integer values. This problem is addressed in a similar manner as in Zhang and Zhu (2007) in the proposed model which will be presented later.

#### Rate-Independency

The Graesser and Cozzarelli model is rate-independent, as is the Bouc-Wen model (Charalampakis and Koumousis, 2008b; Charalampakis, 2015). This means that in strain-controlled experiments, as the ones shown next, the same hysteretic loop will be traced irrespective of the rate of loading. In fact, the input does not even need to be sinusoidal; any input history with the specific sequence of local maxima and minima, i.e., 0 → εmax → −εmax → εmax → . . ., will produce the same result.

#### Twinning Hysteresis

Twinning hysteresis in the Graesser and Cozzarelli model is observed by setting fT = 0 in Equation (2). The base model shown in **Figure 2A** was used in the original paper (Graesser and Cozzarelli, 1991), where E = 28, 500 ksi; Y = 30ksi; α = 0.0197; n = 3. The input history is sinusoidal with εmax = 0.016.

In **Figure 2B**, the effect of varying parameter E is shown. Apart from the initial stiffness, which is clearly controlled by E, the post-elastic stiffness Ey is also affected since Ey = αE/ (1 + α).

In **Figure 2C** the effect of varying parameter Y is examined. Obviously, Y controls the yield stress. Note that the term "yield stress" should be used rather loosely in our context. In SMA, it corresponds to the stress level at which stress induced phase transformation takes place in the material and has nothing to do with plastic deformation. Nevertheless, it will be used for reasons of simplicity.

In **Figure 2D** the effect of varying parameter α is shown. Clearly, α controls the post-elastic stiffness Ey without affecting the initial stiffness E.

In **Figure 2E** the effect of varying the exponential parameter n is shown. As is the case of the Bouc-Wen model, n controls the abruptness of transition between the initial (elastic) and postelastic branches. Small values of n lead to a smooth transition and vice versa. Note that the sensitivity of the model with respect to this parameter becomes negligible when n ≥ 3. In the figure, the responses for n = 3 and n = 5 are practically the same.

#### Superelasticity

Superelasticity in the Graesser and Cozzarelli model is observed by setting fT > 0 in Equation (2). The base model shown in **Figure 3A** was used in the original paper (Graesser and Cozzarelli, 1991), where E = 28,500 ksi; Y = 30 ksi; α = 0.0197; n = 3; fT = 0.07; a¯ = 2,500; and c = 0.001. The input history is sinusoidal with εmax = 0.016. In **Figure 3B**, the effect of varying parameter fT is shown. It is observed that high values of fT reduce the size of the energy-dissipating loops exhibited during the material transformations. Apparently, the sensitivity of the model with respect to fT is high. In **Figure 3C** the effect of varying parameter a¯ is examined. As is the case with n, high values of a¯ do not alter the response significantly. For small values of a¯, the response becomes smooth around the origin, yielding residual deformations upon unloading. In **Figure 3D** the effect of varying parameter c is shown. In general, the sensitivity of the model with respect to this parameter is very small. **Figure 3E** shows the variation of the quantity |ε| c during the experiment. The responses for c = 0.001 and c = 0.00001 are practically the same, because |ε| c ∼= 1 for ε 6= 0. For c = 0.1, the graph of |ε| c takes a form not radically different from the error function erf (·) used in the same expression, i.e., Equation (2). Thus, the response observed for c = 0.1, can also be obtained by appropriate values of fT, n, and a¯. To demonstrate this, the optimum values of fT, n and a¯ have been obtained which produce a similar curve as the one with c = 0.1 (**Figure 3F**) using metaheuristics. In this regard, the term |ε| c in Equation (2) is deemed redundant and will be removed in the proposed model which will be presented later.

### Error Function

The final note on the Graesser-Cozzarelli model refers to the error function itself. As the evaluation of this function is cumbersome [see Equation (3)], its substitution with the hyperbolic tangent is proposed herein. The hyperbolic tangent is remarkably similar to the error function but more comfortable to evaluate, e.g., by using:

$$\tanh\left(\mathbf{x}\right) = 1 - \frac{2}{1 + e^{2\mathbf{x}}}.\tag{5}$$

In **Figure 4**, the similarity of erf (a¯ε) to tanh aˆε is illustrated. Three levels of strain coefficients, i.e., a¯ = {2, 500; 500; 100} are matched manually to aˆ = {3, 000; 625; 125} to produce similar curves.

# THE WILDE—GARDONI—FUJINO MODEL

The Wilde-Gardoni-Fujino model (Wilde et al., 2000) is an extension of the Graesser-Cozzarelli model which takes into account the elastic behavior of the martensite which is activated after a particular strain level. The following equations describe it:

$$\dot{\sigma} = E \left[ \dot{\varepsilon} - |\dot{\varepsilon}| \left( \frac{\sigma - \beta}{Y} \right)^{n} \right] u\_{I} \left( \varepsilon \right) + E\_{m} \dot{\varepsilon} u\_{II} \left( \varepsilon \right)$$

$$+ \left( 3\alpha\_{1} \dot{\varepsilon} \varepsilon^{2} + 2\alpha\_{2} \operatorname{sgn} \left( \varepsilon \right) \dot{\varepsilon} \varepsilon + \alpha\_{3} \dot{\varepsilon} \right) u\_{III} \left( \varepsilon \right), \qquad \text{(6)}$$

$$\beta = E\alpha \left[ \varepsilon\_{in} + f\_T |\varepsilon|^\varepsilon \text{erf} \left( \bar{a}\varepsilon \right) \mu \left( -\varepsilon \dot{\varepsilon} \right) \right],\tag{7}$$

where the functions uI (ε), uII (ε), and uIII (ε) are given by:

$$
\mu\_I(\varepsilon) = 1 - \mu\_{II}(\varepsilon) - \mu\_{III}(\varepsilon), \tag{8}
$$

0.001 (B) variation in *fT* (C) variation in a¯ (D) variation in *c*, (E) value of *|*ε*| c* for various values of *c*, (F) curve-fit of the response obtained with c = 0.1 by altering *fT* , *n,* and *a¯* only (c = 0).

$$\mu\_{\Pi}(\varepsilon) = \begin{cases} 1, & |\varepsilon| \ge \varepsilon\_m \\ 0, & \text{otherwise} \end{cases},\tag{9}$$

$$\mu\_{\text{III}}(\varepsilon) = \begin{cases} 1, & \text{\$\&\$\varepsilon > 0\$ and \$\varepsilon\_1 < |\varepsilon| < \varepsilon\_m \$}, \\ 0, & \text{otherwise} \end{cases}, \tag{10}$$

In this model, the response is divided into separate regions which are activated or deactivated by the flags described by Equations (8)–(10). The strain εm defines the point at which the transformation from austenite to martensite is completed. Beyond this strain, the response is linear elastic with modulus equal to Em, due to the term Emε˙uII (ε) in Equation (6). A smooth transition is achieved during loading only, due to the term activated by uIII (ε).

Although the model is capable of simulating the martensitic phase of the SMA, its usage is cumbersome. For instance, the coefficients a1, a2, and a3 do not have physical representation and plausible value ranges are not easy to establish. As is the case with the Graesser-Cozzarelli model, the term (σ − β) may be negative, which is problematic when n takes real rather than integer values. Smooth transition between austenitic and martensitic phases can be realized with simpler and more elegant relations, as illustrated in the proposed model described later on.

#### THE ZHANG—ZHU MODEL

The Zhang-Zhu model (Zhang and Zhu, 2007) is a modification of the Wilde-Gardoni-Fujino model which aims to enhance the

stability of numerical simulation and speed up the computation time. The following equations describe it:

$$\begin{split} \dot{\sigma}\_{1} &= E \left[ \dot{\varepsilon} - K\_{\varepsilon}(\varepsilon) \left| \dot{\varepsilon} \right| \text{sgn} \left( \sigma - \beta \right) \left( \frac{|\sigma - \beta|}{Y} \right)^{n} \right] u\_{I} \left( \varepsilon \right) \\ &+ E\_{m} \dot{\varepsilon} u\_{II} \left( \varepsilon \right) + \left( E\_{\mathcal{V}} \frac{\varepsilon\_{m} - \varepsilon}{\varepsilon\_{m} - \varepsilon\_{1}} + E\_{m} \frac{\varepsilon - \varepsilon\_{1}}{\varepsilon\_{m} - \varepsilon\_{1}} \right) \dot{\varepsilon} u\_{III} \left( \varepsilon \right), \text{(11)} \\ &\varepsilon\_{0} - \varepsilon\_{1} \left[ \varepsilon + \varepsilon\_{1} \left( \varepsilon + \varepsilon\_{1} \right) \left( \varepsilon \right)^{-} \right] u\_{I} \left( \varepsilon + \varepsilon\_{1} \left( \varepsilon \right)^{-} \right) \dot{\varepsilon} u\_{I} \left( \varepsilon \right) \right] \end{split}$$

$$\beta = E\alpha \left[ \varepsilon\_{in} + f\_T u \left( -\varepsilon \dot{\varepsilon} \right) \text{sgn} \left( \varepsilon \right) \text{g} \left( \bar{a} \varepsilon\_{in} + \text{sgn} \left( \varepsilon \right) b \right) \right], \tag{12}$$

where the functions uI (ε), uII (ε), uIII (ε), K (ε), g (t) are given by:

$$
\mu\_I(\varepsilon) = \lim\_{\varepsilon \to \varepsilon^+} \mu\_{II}(\varepsilon) - \mu\_{III}(\varepsilon) \,, \tag{13}
$$

$$\mu\_{II}(\varepsilon) = \begin{cases} 1, & \text{\$\dot{\varepsilon}\varepsilon > 0\$ and \$|\varepsilon| \ge \varepsilon\_m} \\ 0, & \text{otherwise} \end{cases},\tag{14}$$

$$\mu\_{III}(\varepsilon) = \begin{cases} 1, & \text{ $\dot{\varepsilon}\varepsilon > 0$  and  $\varepsilon\_1 < |\varepsilon| < \varepsilon\_m $ } \\ 0, & \text{otherwise} \end{cases}, \tag{15}$$

$$K(\varepsilon) = \begin{cases} 1, & \text{ $\dot{\varepsilon}\varepsilon > 0$ } \\ \mu\left(\varepsilon\_{\text{in}\varepsilon}\right), & \text{otherwise} \end{cases},\tag{16}$$

$$\lg(t) = 1 - e^{-t^2},\tag{17}$$

and the signum function is defined as

$$\text{sgn}\,(\mathbf{x}) = \begin{cases} 1, & \mathbf{x} > \mathbf{0} \\ 0, & \mathbf{x} = \mathbf{0} \\ -1, & \mathbf{x} < \mathbf{0} \end{cases}.\tag{18}$$

Note that parameters α, E, and n can take different values for loading and loading phases. Thus, the model contains up to fourteen parameters. The main improvements with respect to the Wilde-Gardoni-Fujino are the following:

The following, however, can be listed as disadvantages of the model:

#### PROPOSED MODEL

Based on the observations above, a simple rate-independent uniaxial phenomenological model is proposed herein which is described by the following terms:

$$\dot{\sigma} = \left(1 - s\left(\varepsilon\right)\right) E \left[\dot{\varepsilon} - \left|\dot{\varepsilon}\right| \text{sgn}\left(\sigma - \beta\right) \left(\frac{|\sigma - \beta|}{Y}\right)^n\right] + s\left(\varepsilon\right) E\_m \dot{\varepsilon}, \tag{19}$$

$$\beta = \operatorname{Ea}\left[\varepsilon - \sigma/E \, \, \middle| \, \, \operatorname{tdh}\left(\tilde{a}\varepsilon\right) \, \mu\left(-\varepsilon\dot{\varepsilon}\right)\right],\tag{20}$$

$$s\left(\varepsilon\right) = \frac{\tanh\left(c\left(|\varepsilon| - \varepsilon\_l\right)\right) + 1}{2},\tag{21}$$

where, s(ε)is a smooth logistic function which yields 0 for|ε|≪εt and 1 for |ε| ≫ εt . Note that s(εt) = 0.5, i.e., at this level of strain the weights of the Graesser-Cozzarelli term and the elastic martensitic term in Equation (19) are equal. The coefficient c controls the abruptness of transition. For high values of c, the transition is abrupt, as evidenced by **Figure 5A**.

In total, the model contains only nine parameters (E, Y, α, n, fT, a¯, Em, εt , c) as opposed to the thirteen parameters of the Wilde-Gardoni-Fujino model (E, Em, Y, α, n, fT, c, a¯, εm, ε1, α1, α2, α3) and the eleven to fourteen parameters of the Zhang-Zhu model (El , Eu, Em, Y, EY, αl , αu, nl , nu, fT, a¯, εm, ε1, b). A brief description of the model parameters and their effect on the overall response, based on the evidence examined earlier, is summarized in **Table 1**. The last three parameters refer to the martensitic phase.

TABLE 1 | Brief description of proposed model parameters and their effect on the overall response.

(B) Zhang and Zhu (2007), (C) Zhuang et al. (2016), (D) Ozbulut et al. (2010), (E) Auricchio et al. (2003), (F) Ozbulut and Hurlebaus (2010).

Despite its simplicity, the proposed model can accurately capture all the pertinent characteristics of the response curve. System identification based on metaheuristics produces excellent fit with experimental data obtained from the literature, as will be demonstrated. The optimum (best) parameter values were evaluated by Differential Evolution, a powerful metaheuristic algorithm (Storn and Price, 1997). The DE/rand/1/bin configuration was used with population size P = 50, F = 0.5 and Cr = 0.9. Ten independent runs were conducted with different random seeds. Each run was terminated after 20000 function evaluations. The sum of squares of the difference between the measured and the predicted time history of stress is used as the objective function to be minimized. Cast in discrete form, this function can be written as

$$OF\left(\mathbf{p}\right) = \sum \left(\sigma\left(t\_i\right) - \widehat{\sigma}\left(t\_i \mid \mathbf{p}\right)\right)^2,\tag{22}$$

where **p** is the parameter vector.

TABLE 2 | Statistical analysis of identified parameters of the proposed model for the experimental data obtained from Zhang and Zhu (2007) shown in Figure 5B.

TABLE 3 | Statistical analysis of identified parameters of the proposed model for the experimental data obtained from Zhuang et al. (2016) shown in Figure 5C.

In **Figure 5B** the case of tensional tests on cold-drawn nitinol wires with a diameter of 0.58 mm obtained from Zhang and Zhu (2007) is shown. A statistical analysis of the identified parameters is summarized in **Table 2**. It is noted that the initial side constraints are quite relaxed, yet the identification process was found to be straightforward. It is also observed that some parameters are very sensitive to the considered experimental data set, i.e., parameters εt , Y, Em, fT, and c, while the others are not so. In general, however, all runs were obviously guided to the same optimum region of the search space, which is a strong indication of a well-behaved model with uniquely defined parameters, each one with a distinct role in the response curve. If needed, the insensitive parameters can be accurately identified with other experiments which focus on the pertinent characteristics of the response curve, e.g., parameter E with a linear elastic experiment with small strains. Improved results can also be obtained by gradual tightening of the side constraints based on trial analyses, as in Charalampakis and Koumousis (2008a), but this lies beyond the scope of the present study and will be pursued in future research.

**Figure 5C** shows excellent fit in the case of the single SMA helical spring with hook-like ends obtained from Zhuang et al. (2016). Being phenomenological, the proposed model can be easily used with force-displacement curves instead of stress-strain curves, with the necessary semantic changes concerning parameters E and Y. By inspection, it is clear that the experimental data does not include any stress-induced martensitic phase. Therefore, the pertinent parameters (Em, εt , c) cannot be identified. This fact is manifested clearly in the results of the identification, where the metaheuristic algorithm selects arbitrary values for these parameters in order to disable them (i.e., εt >0.037, any Em, and large c). In this specific case, we can use s(ε) = 0 instead of Equation (21), or repeat the experiment with a larger displacement to record martensitic phase data. A statistical analysis of the identified parameters is summarized in **Table 3**, where small coefficients of variation are observed.

**Figure 5D** shows an excellent fit for the case of NiTi SMA wires with a diameter of 1 mm. The data refers to the response of different strain levels at a frequency of 1 Hz and has been obtained from Ozbulut et al. (2010). A value of 0.005 has been subtracted from all strains to compensate for the lack of complete tautness. The wires have been trained with a set of 10 load cycles with a strain amplitude of 6% at 0.04 Hz. A statistical analysis of the identified parameters is summarized in **Table 4**, where, again, small coefficients of variation are observed.

In **Figure 5E**, the experimental stress-strain response of an untrained commercial Ni-Ti wire with a diameter of 2.01 mm is shown. The loading cycles are characterized by an increasing value of the maximum elongation that varies from 2 to 8%. The data has been obtained from Auricchio et al. (2003), and an excellent fit can be observed for the second and last cycle, indicatively. Apparently, the values differ significantly, which is the result of the progressive decrement of the initial and final stress thresholds of the phase transformations during the training of a specimen (Auricchio et al., 2003). In similar fashion, in **Figure 5F** very good fit is observed for the experimental data obtained from Ozbulut and Hurlebaus (2010). The data refer to NiTi wires with a diameter of 1.5 mm which are subjected to tensile tests with various temperatures and strain rates. The identified model parameters for the latter cases are summarized in **Table 5**. By appropriate interpolation of model parameters, it may be possible that aspects which are not covered presently, TABLE 4 | Statistical analysis of identified parameters of the proposed model for the experimental data obtained from Ozbulut et al. (2010) shown in Figure 5D.

TABLE 6 | Statistical analysis of identification results for the Wilde-Gardoni-Fujino model with data obtained from Zhang and Zhu (2007).

such as strain-rate and temperature dependency, or training process, are taken into account. This, however, lies beyond the scope of the present study and will be pursued in future research.

In order to demonstrate the advantages of the proposed model, the same identification procedure is applied to the Wilde-Gardoni-Fujino model with data obtained from Zhang and Zhu (2007). The results are summarized in **Table 6**, where great variability in parameter values and quality of solutions can be observed.

# CONCLUSIONS

In this paper, three specific uniaxial phenomenological models commonly used for the description of a Shape Memory Alloy (SMA) behavior were examined in detail, and a new simple rate-independent model was proposed which addresses all issues in a unified manner. From the presented analysis and the numerical results, the following main conclusions can be drawn:

# REFERENCES

# AUTHOR CONTRIBUTIONS

AC had the research idea, drafted the article, and contributed to the derivation of the numerical examples. GT contributed to the conception of the work and interpretation of the results. Both authors contributed to the writing of the manuscript.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a shared affiliation, though no other collaboration, with one of the authors AC at time of review.

Copyright © 2018 Charalampakis and Tsiatas. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# A Meshless Solution to the Vibration Problem of Cylindrical Shell Panels

#### Aristophanes J. Yiotis\* and John T. Katsikadelis

*School of Civil Engineering, National Technical University of Athens, Athens, Greece*

The Meshless Analog Equation Method (MAEM) is a purely mesh-free method for solving partial differential equations (PDEs). In the present study, the method is applied to the dynamic analysis of cylindrical shell structures. Based on the principle of the analog equation, MAEM converts the three governing partial differential equations in terms of displacements into three uncoupled substitute equations, two of 2nd order (Poisson's) and one of 4th order (biharmonic), with fictitious sources. The fictitious sources are represented by series of Radial Basis Functions (RBFs) of multiquadric (MQ) type, and the substitute equations are integrated. The integration allows the representation of the displacements by new RBFs, which approximate the displacements accurately and also their derivatives involved in the governing equations. By inserting the approximate solution in the governing differential equations and taking into account the boundary and initial conditions and collocating at a predefined set of mesh-free nodal points, we obtain a system of ordinary differential equations of motion. The solution of the system gives the unknown time-dependent series coefficients and the solution to the original problem. Several shell panels are analyzed using the method, and the numerical results demonstrate its efficiency and accuracy.

#### Edited by:

*Vagelis Plevris, OsloMet—Oslo Metropolitan University, Norway*

#### Reviewed by:

*Andreas Kampitsis, Imperial College London, United Kingdom Aristotelis E. Charalampakis, National Technical University of Athens, Greece*

#### \*Correspondence:

*Aristophanes J. Yiotis<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Computational Methods in Structural Engineering, a section of the journal Frontiers in Built Environment*

> Received: *13 March 2018* Accepted: *04 July 2018* Published: *12 September 2018*

#### Citation:

*Yiotis AJ and Katsikadelis JT (2018) A Meshless Solution to the Vibration Problem of Cylindrical Shell Panels. Front. Built Environ. 4:40. doi: 10.3389/fbuil.2018.00040* Keywords: MAEM, Meshless Analog Equation Method, cylindrical shells, dynamic analysis, radial basis functions, partial differential equations

#### INTRODUCTION

Thin shell structures have an outstanding efficiency in fully utilizing the structural material and have been extensively used in many engineering applications including aircraft structures, pressure vessels, and others. Static and dynamic analysis is essential for the analysis and design of shell structures. Various numerical methods, such as the Finite Difference Method (FDM) and especially the Finite Element Method (FEM) have been used (Lee and Han, 2001) for the dynamic analysis of linear elastic thin shells characterized by complex geometry, loading and boundary conditions. Both methods have been employed successfully for the solution of a variety of static and dynamic shell problems. The Boundary Element Method (BEM) is an efficient alternative to the domain type methods, especially for thin elastic shallow shells (Beskos, 1991), or combined with the AEM for cylindrical shells (Yiotis and Katsikadelis, 2000).

Such methods require the generation of a mesh which can be an incredibly tedious and timeconsuming process, while their convergence rate is of 2nd order (Cheng et al., 2003). On the other hand, Meshless Methods (MMs) present an attractive alternative to FEM or BEM, especially for shell structures that are complex regarding both the governing equations and the geometry representation. Comprehensive descriptions of different MMs are presented by Liu (2002); Liu and Gu (2005) and in a review paper by Nguyen et al. (2008).

There are several papers on dynamic analysis of shells using MM. Homogeneous shells are studied using various versions of MM in Liu et al. (2006), Ferreira et al. (2006b), and Dinis et al. (2011). Functionally graded (FG) cylindrical thin shells have also been treated by this method (Ferreira et al., 2006a; Zhao et al., 2009; Roque et al., 2010), as well as thick cylindrical shells (Pilafkan et al., 2013) have been analyzed by this method.

The mesh-free multiquadric radial basis functions (MQ-RBFs) method presented in Kansa (2005) has attracted the interest of researchers, due to its exponential convergence and its easiness of implementation. The significant drawbacks of the method are the ill-conditioning of the coefficient matrix and the inability to accurately approximate the derivatives of the sought solution which renders the method inappropriate for a strong formulation of the problem. These drawbacks of the standard MQ-RBF method, are overcome by a new RBF method presented recently by Katsikadelis (2006, 2008a,b, 2009) and Yiotis and Katsikadelis (2008, 2013, 2015a,b). Another critical issue is the implementation of multiple boundary conditions for equations of order higher than 2nd. In this investigation, the δ-technique is employed (Jang et al., 1989) for the 4th order equation. The problem of multiple boundary conditions is not present when the shell is modeled as a 3D body (Katsikadelis and Platanidi, 2007).

In this paper, the MAEM is extended to the dynamic problem of cylindrical shell panels as described by section MAEM Solution. A first approach to this problem was attempted in a previous work (Yiotis and Katsikadelis, 2015b), where some preliminary results only for the eigenfrequency analysis were presented. In section Problem Statement, the statement of the problem is presented, while several example problems are worked out in section Numerical Examples, which illustrate the applicability of the method and demonstrate its efficiency and section Conclusions contains certain conclusions drawn from this investigation.

#### PROBLEM STATEMENT

We consider a thin cylindrical shell with parametric lines x (s = const.) and s (x = const.) which are assumed to be lines of curvature, as well; x is measured along the x lines of the shell and s along the s lines, while z is measured along the normal to the middle surface of the shell, as shown in **Figure 1**. R is the radius of curvature and h is the thickness.

In this investigation we use the Flügge equations for the thin shell theory, based on the following assumptions (Love, 1944):

The first assumption defines the meaning of "thin shells," whereas the second one implies that all calculations refer to the original undeformed configuration and subsequently leads to linear differential equations. Further, the assumption z/R << 1 is adopted in deriving the stress resultants in integrating the stresses through the thickness of the shell. The 4th assumption is known as Kirchhoff's hypothesis yielding

$$
\gamma\_{\pi z} = 0,\tag{1a}
$$

$$
\gamma\_{\varepsilon} = 0,\tag{1b}
$$

$$e\_{\sharp} = 0,\tag{1c}$$

which implies σxz = σsz = 0 (Leissa, 1973).

The equations of motion for the case of the thin cylindrical shell can be derived using Hamilton's principle as follows

$$\delta \left[ \int\_{t\_0}^{t\_1} (\Pi - K) dt - \int\_{t\_0}^{t\_1} W\_{nc} dt \right] = 0,\tag{2}$$

where 5 is the total potential energy given by

$$
\Pi = U\_0 - W\_1 - W\_2,\tag{3}
$$

in which U0, is the strain energy

$$U\_0 = \frac{1}{2} \int\_{\underline{v}} (\sigma\_{\underline{x}} \mathbf{e}\_{\underline{x}} + \sigma\_{\underline{s}} \mathbf{e}\_{\underline{s}} + \sigma\_{\underline{x}\boldsymbol{\varepsilon}} \boldsymbol{\chi}\_{\underline{x}\boldsymbol{\varepsilon}} + \sigma\_{\underline{x}\boldsymbol{\varepsilon}} \boldsymbol{\chi}\_{\underline{x}\boldsymbol{\varepsilon}} + \sigma\_{\underline{s}\boldsymbol{\varepsilon}} \boldsymbol{\chi}\_{\underline{x}\boldsymbol{\varepsilon}}) d\boldsymbol{x} d\boldsymbol{s} d\boldsymbol{z}, \quad \text{(4)}$$

and W1, W2 the works produced by the loading and the boundary forces, i.e.,

$$W\_1 = \int\_S (q\_\mathbf{x}\mu + q\_\mathbf{s}\nu + q\_\mathbf{z}\nu)d\mathbf{x}ds,\tag{5}$$

$$W\_2 = \int (\bar{N}\_\mathbf{x}\mu + \bar{N}\_\mathbf{x}\nu + \bar{Q}\_\mathbf{x}\nu + \bar{M}\_\mathbf{x}\theta\_\mathbf{x} + \bar{M}\_\mathbf{x}\theta\_\mathbf{s})ds$$

$$\int\_{\infty}^{\bar{\zeta}} + \int\_{\infty} (\bar{N}\_{s\mathbf{x}}\mu + \bar{N}\_{s}\nu + \bar{Q}\_{s}\nu + \bar{M}\_{\mathbf{x}\mathbf{x}}\theta\_{\mathbf{x}} + \bar{M}\_{s}\theta\_{s})d\mathbf{x}, \quad \text{(6)}$$

In Equation (4) σx, σs are the normal stresses, σxs the tangential shear stress, σxz, σsz the transverse (in the z direction) shear stresses and ex, es , γxs, γxz, γsz the respective strains at an arbitrary point of shell cross-section.

In Equation (5), qx(t), qs(t), and qz(t) are the three components of the loading in the axial, circumferential and normal to the middle surface directions, respectively, while u, ν, and w represent the axial, circumferential and normal displacements at the middle surface of the shell.

In Equation (6) the quantities N¯ x, N¯ xs, Q¯ x, M¯ x, M¯ xs, and N¯ sx, N¯ s , Q¯ s , M¯ sx, M¯ s denote prescribed boundary forces along an edge (x = const.) and an edge (s = const.) respectively; u, ν and w represent the axial, circumferential and normal displacements at the boundary and θx, θs are the rotations of the normal to the middle surface about the s and x axes respectively.

Furthermore, in Equation (2) the quantity K is the kinetic energy of the body and is given regarding the shell variables as

$$K = \frac{h\rho}{2} \int\_{\chi} \int\_{s} \left[ \dot{u}^2 + \dot{\nu}^2 + \dot{\nu}^2 + \frac{h^2}{12} (\dot{\theta}\_x^2 + \dot{\theta}\_s^2) \right] dxds,\tag{7}$$

where ρ is the mass density of the material of the shell.

The quantity δWnc represents the work of the damping forces, non-conservative forces, due to the virtual displacements and is given by the relation

$$
\delta W\_{nc} = \int\_{\mathcal{X}} \int\_{s} (\eta \dot{u} \delta u + \eta \dot{v} \delta v + \eta \dot{w} \delta w) dx ds,\tag{8}
$$

where η is the damping coefficient.

Neglecting the contribution from the rotatory inertia terms ρh 3 θx/12 and ρh 3 θs/12 in Equation (7), inserting Equation (3) and taking the variation (Katsikadelis, 2016), we obtain the Flügge type differential equations (Flügge, 1962; Kraus, 1967), in terms of the displacements as well as the associated boundary and initial conditions

#### (a) Differential equations

u,xx + 1 − ν 2 u,ss + 1 + ν 2 ν,xs + ν R w,x − h 2 12 1 R w,xxx − 1 − ν 2 R w,xs R + u,s R2 ,s − ηu˙ = − 1 − ν 2 Eh (qx − ρhu¨), (9a) ν,ss + 1 − ν 2 ν,xx + 1 + ν 2 u,xs + w R ,s + h 2 12 1 R2 3(1 − ν) 2 ν,xx − (3 − ν) 2 Rw,xxs − R,s(w,ss + w R2 + R,s R2 ν) − ην˙ = − 1 − ν 2 Eh (qs − ρhν¨), (9b) ∇ 4w + w,ss R2 + w R2 ,ss + w R4 − 1 R u,xxx + 1 − ν 2 u,xs R ,s − 3 − ν 2 ν R ,xxs + R,s R2 ν ,ss + R,s R4 ν + 12 h 2 1 R ν,s + w R + νu,x + ηw˙ = − 12(1 − ν 2 ) Eh3 −qz + ρhw¨ , (9c)

where ∇ 4 = ∂ 4 ∂x 4 + 2∂ 4 ∂x 2∂s 2 + ∂ 4 ∂s 4 is the biharmonic operator E is the modulus of elasticity and ν is Poisson's ratio.

(b) The boundary conditions (Kraus, 1967) On a curved edge (x = 0 or x = l)

$$
\mu = \bar{\mu} \operatorname{or} N\_{\mathfrak{x}} = \bar{N}\_{\mathfrak{x}}.\tag{10a}
$$

$$\nu = \bar{\nu} \text{ or } T\_{\infty} = \bar{T}\_{\infty}, \tag{10b}$$

$$
\bar{w} = \bar{w} \text{ or } V\_{\bar{x}} = \bar{V}\_{\bar{x}}, \tag{10c}
$$

$$\theta\_{\mathbf{x}} = \bar{\theta}\_{\mathbf{x}} \text{ or } M\_{\mathbf{x}} = \bar{M}\_{\mathbf{x}} \left( \theta\_{\mathbf{x}} = -\frac{\partial \mathcal{W}}{\partial \mathbf{x}} \right). \tag{10d}$$

On a straight edge (s = 0 or s = a)

$$
\mu = \bar{u} \text{ or } T\_{s\mathbf{x}} = \bar{T}\_{s\mathbf{x}}.\tag{10e}
$$

$$
\nu = \bar{\nu} \text{ or } \mathcal{N}\_{\mathfrak{s}} = \bar{\mathcal{N}}\_{\mathfrak{s}}.\tag{10f}
$$

$$
\bar{w} = \bar{w} \text{ or } V\_s = \bar{V}\_s,\tag{10g}
$$

$$\theta\_s = \bar{\theta}\_s \text{ or } M\_s = \bar{M}\_s \left( \theta\_s = \frac{\nu}{R} - \frac{\partial w}{\partial s} \right). \tag{10h}$$

Besides, the following corner conditions must be satisfied (Leissa, 1973)

$$\boldsymbol{\nu} = \bar{\boldsymbol{w}} \text{ or } (\boldsymbol{M}\_{\text{xs}} - \boldsymbol{M}\_{\text{sx}})\_k = \bar{\boldsymbol{F}}\_k, \quad k = 1, 2, 3, 4. \tag{10i}$$

(c) The initial conditions

$$\mathcal{w}(\mathbf{x},0) = \mathcal{g}\_{\mathcal{I}}(\mathbf{x}), \quad \dot{\mathcal{w}}(\mathbf{x},0) = h\_{\mathcal{I}}(\mathbf{x}), \tag{11a-b}$$

$$u(\mathbf{x},0) = \mathbf{g}\_1(\mathbf{x}), \quad \dot{u}(\mathbf{x},0) = h\_1(\mathbf{x}), \tag{11c\text{-}d}$$

$$\nu(\mathbf{x},0) = \operatorname{g}\_2(\mathbf{x}), \quad \dot{\nu}(\mathbf{x},0) = h\_2(\mathbf{x}), \tag{11e-f}$$

where gi(x), hi(x) (i = 1, 2, 3) are specified functions and (x) = (x,s).

The stress resultants Nx, Ns , Nxs, Nsx, Mx, Ms , Mxs, Msx, Qx, Qs are expressed in terms of the displacements as

$$N\_{\rm x} = \frac{Eh}{1 - \nu^2} \left[ \mu\_{\rm xx} + \nu \left( \nu\_{\rm xs} + \frac{\nu}{R} \right) - \frac{h^2}{12} \frac{1}{R} \nu\_{\rm xxx} \right],\tag{12a}$$

$$\begin{split} N\_{s} &= \frac{Eh}{1 - \nu^{2}} \Big[ \nu\_{s\text{s}} + \frac{\nu}{R} + \nu u\_{,\text{x}} \\ &+ \frac{h^{2}}{12} \frac{1}{R} \Big( \nu\_{\text{s}\text{s}} + \frac{\nu}{R^{2}} + \frac{R\_{\text{s}\text{s}}}{R^{2}} \nu \Big) \Big], \end{split} \tag{12b}$$

$$N\_{\infty} = \frac{Eh}{2\left(1+\nu\right)} \left[ u\_{\text{ss}} + \nu\_{\text{ss}} - \frac{h^2}{12} \frac{1}{R} \left( \nu\_{\text{ss}} - \frac{1}{R} \nu\_{\text{ss}} \right) \right], \quad \text{(12c)}$$

$$N\_{\rm xx} = \frac{Eh}{2\left(1+\nu\right)} \left[ u\_{\rm zs} + v\_{\rm xx} + \frac{h^2}{12} \frac{1}{R} \left( \nu\_{\rm xx} + \frac{1}{R} u\_{\rm zs} \right) \right], \quad \text{(12d)}$$

$$M\_{\rm x} = -D\{\boldsymbol{\omega}\_{\rm xxx} + \nu \left[\boldsymbol{\omega}\_{\rm xx} - \left(\frac{\nu}{R}\right)\_{\rm x}\right] - \frac{1}{R}\boldsymbol{\omega}\_{\rm x}\},\tag{12e}$$

$$M\_{\rm s} = -D\left(\omega\_{\rm sst} + \frac{\nu}{R^2} + \frac{R\_{\rm s}}{R^2}\nu + \nu w\_{\rm sxx}\right),\tag{12f}$$

$$M\_{\rm xs} = -D\left(1 - \nu\right) \left(\nu\_{\rm ss} - \frac{1}{R}\nu\_{\rm xs}\right),\tag{12g}$$

$$M\_{\rm sx} = -D\frac{(1-\nu)}{2}\left(2\nu\_{\rm sx} - \frac{\nu\_{\rm sx}}{R} + \frac{\mu\_{\rm s}}{R}\right),\tag{12h}$$

$$Q\_{\rm x} = -D\left[\begin{matrix}\boldsymbol{w}\_{\rm xxx} + \boldsymbol{w}\_{\rm xxx} - \frac{\boldsymbol{u}\_{\rm xxx}}{R} + \frac{(1-\nu)}{2}\left(\frac{\mu\_{\rm s}}{R}\right)\_{\rm s} \end{matrix}\right]$$

$$-\frac{1+\nu}{2R}\left(\frac{\nu}{R}\right)\_{\text{xx}}\right],\tag{12i}$$

$$Q\_s = -D\left[\left.\nu\_{\text{sys}} + \left.\nu\_{\text{xxx}} + \left(\frac{\nu}{R}\right)\_{,s} - \frac{(1-\nu)}{R}\nu\_{\text{xxx}}\right]\right]$$

$$
\begin{bmatrix}
\frac{\partial \boldsymbol{\nu}}{\partial \boldsymbol{\nu}} & \frac{\partial \boldsymbol{\nu}}{\partial \boldsymbol{\nu}} & \frac{\partial \boldsymbol{\nu}}{\partial \boldsymbol{\nu}}
\end{bmatrix}
\begin{bmatrix}
\boldsymbol{R}^{2} \\
\boldsymbol{R}^{2}
\end{bmatrix}^{\prime\prime} = \boldsymbol{R}^{\prime}
$$

$$
+ \left(\frac{\boldsymbol{R}\_{\rm s}}{\boldsymbol{R}^{2}} \boldsymbol{\nu}\right)\_{\rm s} \begin{bmatrix}
\boldsymbol{\nu} \\
\boldsymbol{\nu}
\end{bmatrix},
\tag{12j}
$$

where D = E/12(1 − ν 2 ) and, Txs Vx the effective tangential membrane and transverse shear forces at the edges x = 0, l given as

$$T\_{\rm xs} = N\_{\rm xs} + \frac{M\_{\rm xs}}{R} = \frac{Eh}{2(1+\upsilon)} \left[ u\_{\rm ss} + \upsilon\_{\rm xs} + \frac{h^2}{4} \frac{1}{R} \left( \frac{\upsilon\_{\rm xs}}{R} - \upsilon\_{\rm xs} \right) \right],\tag{13a}$$

$$\begin{split} V\_{\rm x} &= Q\_{\rm x} + \frac{\partial M\_{\rm xs}}{\partial s} = -D \left[ \upsilon\_{\rm xxxx} + (2-\upsilon)\upsilon\_{\rm xs} \right. \\ &\quad - \frac{u\_{\rm xx}}{R} + \frac{(1-\upsilon)}{2} \left( \frac{u\_{\rm s}}{R} \right)\_{\rm s} - \frac{3-\upsilon}{2} \left( \frac{\upsilon\_{\rm xs}}{R} \right)\_{\rm s} \end{split} \tag{13b}$$

Similarly, Tsx and Vs represent the effective tangential membrane, and transverse shear force at the edges s = 0, a and are given as

$$T\_{sx} = N\_{sx} = \frac{Eh}{2(1+\nu)}[u\_{ss} + \nu\_{sx} + \frac{h^2}{12}\frac{1}{R}(w\_{\times s} + \frac{1}{R}u\_{s})], \quad \text{(14a)}$$

$$V\_s = Q\_s + \frac{\partial M\_{sx}}{\partial x} = -D[w\_{\times sx} + (2-\nu)w\_{\times x\infty} - \frac{3(1-\nu)}{2}\frac{\nu\_{\times x}}{R}]$$

$$+ \left(\frac{R\_{s\,\nu}}{R^2}\right)\_{\,\nu} + \left(\frac{\nu}{R^2}\right)\_{\,\nu} + \frac{(1-\nu)}{2}\frac{u\_{\times x}}{R}]. \tag{14b)}$$

#### MAEM SOLUTION

MAEM (Katsikadelis, 2002) is used for the solution of the initial boundary problem (9), (10), (11), as shown in the following. Let u, ν and w be the solution to the problem. Since Equations (9) are of the 2nd order with respect to u, ν and of the 4th order with respect to w, the analog equations which are convenient to use are

$$
\nabla^2 u = b\_1(\mathbf{x}, t), \quad \nabla^2 v = b\_2(\mathbf{x}, t), \quad \nabla^4 w = b\_3(\mathbf{x}, t), \tag{15a, \mathbf{b}, \mathbf{c}})
$$

where bi = bi(x, t) (i = 1, 2, 3) are unknown fictitious sources depending on time, which, however, is treated as a parameter, i.e., Equations (15a,b,c) are quasi-static, treated as static at each instant. The fictitious sources can be established as follows.

The fictitious sources are approximated with MQ-RBFs series. Thus we have

$$\begin{aligned} \nabla^2 u &\simeq \sum\_{j=1}^{K+L} a\_j^{(1)}(t) f\_{\hat{j}}(r), & \nabla^2 v &\simeq \sum\_{j=1}^{K+L} a\_j^{(2)}(t) f\_{\hat{j}}(r), \\ &\nabla^4 w &\simeq \sum\_{j=1}^{K+2L} a\_j^{(3)}(t) f\_{\hat{j}}(r) \end{aligned} \tag{16a,b,c}$$

where c is the shape parameter; K, L represent the number of collocating points inside and on Ŵ, respectively; fj(r) = √ r 2 + c 2 , r = x − xj  (see **Figure 2**), and a (1) j , a (2) j , a (3) j timedependent coefficients to be determined. Note that, while the derivatives of the membrane displacements u, ν are collocated at K domain and L boundary points, the derivatives of the normal displacement w according to the δ-technique (Ferreira et al., 2005) are collocated in K domain and 2L boundary nodal points placed on the auxiliary boundary Ŵ˜ at a small distance δ from the actual one.

Equations (16) can be directly integrated to yield

$$u \simeq \sum\_{j=1}^{K+L} a\_j^{(1)}(t)\hat{u}\_j, \quad \nu \simeq \sum\_{j=1}^{K+L} a\_j^{(2)}(t)\hat{v}\_j, \quad w \simeq \sum\_{j=1}^{K+2L} a\_j^{(3)}(t)\hat{w}\_j,\tag{17a,b,c}$$

where uˆj(r), νˆj(r), wˆ j(r) are solutions of the equations

$$
\nabla^2 \hat{u}\_j = f\_j(r), \quad \nabla^2 \hat{v}\_j = f\_j(r), \quad \nabla^4 \hat{w}\_j = f\_j(r). \tag{18a,b,c}
$$

Since the functions fj(r) depend only on the radial distance r, the solution of Equations (18) can be obtained after writing them in polar coordinates. For the 2nd order equations, we have

$$\nabla^2 \hat{u}\_j = \frac{1}{r} \frac{d}{dr} \left( r \frac{d\hat{u}\_j}{dr} \right) = f\_j(r), \tag{19}$$

which after integration gives

$$
\hat{u}\_{\hat{\jmath}} = \frac{1}{9} f\_{\hat{\jmath}}^3 + \frac{1}{3} f\_{\hat{\jmath}} c^2 - \frac{c^3}{3} \ln\left(c + f\_{\hat{\jmath}}\right) + G\_1 \ln(r) + H\_1. \tag{20}
$$

Similarly, we have

$$
\hat{\upsilon}\_{\hat{\jmath}} = \frac{1}{9} f\_{\hat{\jmath}}^3 + \frac{1}{3} f\_{\hat{\jmath}} c^2 - \frac{c^3}{3} \ln\left(c + f\_{\hat{\jmath}}\right) + G\_2 \ln(r) + H\_2. \tag{21}
$$

The regularity condition at r = 0 demands G1 = G2 = 0. The remaining constants H1, H2 together with the shape

TABLE 1 | Eigenfrequency parameters of the shell in Example 1.

parameter c, if not arbitrarily specified, can be determined with an optimization procedure, such as to ensure the regularity of coefficients matrix (control of the condition number) and the error minimization. It has been shown that the coefficient matrix resulting from the new RBFs is always invertible (Sarra, 2006), and as a result, we take in this analysis H1 = H2 = 0 for convenience. Thus only c, the shape parameter is involved in the error minimization procedure (Katsikadelis, 2008a).

For the 4th order equation, one can write

$$
\nabla^4 \hat{\boldsymbol{w}} = \nabla^2 (\nabla^2 \hat{\boldsymbol{w}}) = f\_{\hat{\jmath}}.\tag{22}
$$

Integrating Equation (22), and removing the singular terms and the terms including the arbitrary constants (Yao et al., 2010) yield

$$\hat{\boldsymbol{w}}\_{\circ} = -\frac{7}{60}\boldsymbol{c}^{4}\boldsymbol{f}\_{\circ} + \frac{2}{45}\boldsymbol{c}^{2}\boldsymbol{f}\_{\circ}^{3} + \frac{1}{225}\boldsymbol{f}\_{\circ}^{5} + \frac{2\boldsymbol{c}^{2} - 5\boldsymbol{r}^{2}}{60}\boldsymbol{c}^{3}\ln(\boldsymbol{c} + \boldsymbol{f}\_{\circ}) + \frac{1}{12}\boldsymbol{r}^{2}\boldsymbol{c}^{3}.\tag{23}$$

Direct differentiation of Equations (17) obtains the derivatives of the displacements involved in the governing equations (9a,b,c).

$$\begin{split} u\_{\boldsymbol{i},ikl} &\simeq \sum\_{j=1}^{K+L} a\_j^{(1)}(\boldsymbol{t}) \hat{u}\_{\boldsymbol{j},ikl}, \quad \boldsymbol{\nu}\_{\boldsymbol{:}ikl} \simeq \sum\_{j=1}^{K+L} a\_j^{(2)}(\boldsymbol{t}) \hat{v}\_{\boldsymbol{j},ikl}, \\ \boldsymbol{w}\_{\boldsymbol{i}kl} &\simeq \sum\_{j=1}^{K+2L} a\_j^{(3)}(\boldsymbol{t}) \hat{w}\_{\boldsymbol{j}\boldsymbol{:}ikl}, \end{split} \tag{24a,b,c}$$

where i, k, l stand for x,s.

Furthermore, the derivatives of the displacements with respect to time can also be obtained by direct differentiation of Equations (17). Thus we have

$$\dot{\boldsymbol{u}} \simeq \sum\_{j=1}^{K+L} \dot{\boldsymbol{a}}\_j^{(1)}(t)\hat{\boldsymbol{u}}\_j, \quad \dot{\boldsymbol{v}} \simeq \sum\_{j=1}^{K+L} \dot{\boldsymbol{a}}\_j^{(2)}(t)\hat{\boldsymbol{v}}\_j, \dot{\boldsymbol{w}} \simeq \sum\_{j=1}^{K+2L} \dot{\boldsymbol{a}}\_j^{(3)}(t)\hat{\boldsymbol{w}}\_j,\tag{25a,b,c}$$

$$\ddot{u} \simeq \sum\_{j=1}^{K+L} \ddot{a}\_j^{(1)}(t)\hat{u}\_j, \quad \ddot{v} \simeq \sum\_{j=1}^{K+L} \ddot{a}\_j^{(2)}(t)\hat{v}\_j, \quad \ddot{w} \simeq \sum\_{j=1}^{K+2L} \ddot{a}\_j^{(3)}(t)\hat{w}\_j. \tag{26a,b,c}$$

Collocating Equations (9) at the K nodal points inside and the four boundary conditions, Equations (10), at the L boundary nodal points (**Figure 2**) using the well-known δtechnique for multiple boundaries (Yiotis and Katsikadelis, 2015a), and inserting Equations (17) and (24) to (26) in the resulting expressions, a system of ordinary differential equations is obtained, namely

$$\mathbf{M}\ddot{\mathbf{a}} + \mathbf{C}\dot{\mathbf{a}} + \mathbf{K}\mathbf{a} = \mathbf{g},\tag{27}$$

where **M**, **C**, and **K** are known square matrices having dimension 3K +4L; **g** is a vector including the 3K values of the external load g(**x**, t) and **a** is the vector of the 3K + 4L values of the unknown time-dependent coefficients a (1) j (t), a (2) j (t), a (3) j (t).

Equation (27) is the semi-discretized equation of motion of the cylindrical shell with **M**, **C**, and **K** representing the TABLE 2 | Eigenfrequency parameters of the shell in Example 2.

generalized mass, damping and stiffness matrices, respectively. It can be solved numerically, using any time step integration technique to establish the time-dependent unknown coefficients. Here the method presented in Katsikadelis (2014a,b) is employed. The initial conditions of Equation (27) result from Equations (17) and (25) on the basis of Equations (11) as follows:

$$\mathbf{a}^{(1)}(\mathbf{0}) = \hat{\mathbf{u}}^{-1} \mathbf{g}\_1(\mathbf{x}), \quad \dot{\mathbf{a}}^{-1}(\mathbf{0}) = \hat{\mathbf{u}}^{-1} \mathbf{h}\_1(\mathbf{x}), \tag{28a}$$

$$\mathbf{a}^{(2)}(\mathbf{0}) = \hat{\boldsymbol{\upsilon}}^{-1} \mathbf{g}\_2(\mathbf{x}), \quad \dot{\mathbf{a}}^{(2)}(\mathbf{0}) = \hat{\boldsymbol{\upsilon}}^{-1} \mathbf{h}\_2(\mathbf{x}), \tag{28b}$$

$$\mathbf{a}^{(3)}(\mathbf{0}) = \hat{\mathbf{w}}^{-1} \mathbf{g}\_{\mathfrak{z}}(\mathbf{x}), \quad \dot{\mathbf{a}}^{(3)}(\mathbf{0}) = \hat{\mathbf{w}}^{-1} \mathbf{h}\_{\mathfrak{z}}(\mathbf{x}).\tag{28c}$$

Once the coefficients a (1) j (t), a (2) j (t), a (3) j (t) have been computed, the field functions u, ν,w, their derivatives, and the stress resultants can be evaluated from Equations (17), (24) to (26) and (12) to (14).

TABLE 3 | Eigenfrequency parameters of the shell in Example 3: case (a).

For free vibrations it is **C** = **g**(**x**, t) = **0** and the equation of motion, Equation (27), takes the form

$$\mathbf{M}\ddot{\mathbf{a}} + \mathbf{K}\mathbf{a} = \mathbf{0},\tag{29}$$

while the essential boundary conditions, Equations (10), become homogeneous.

By setting

$$\mathbf{a}(t) = \alpha e^{i\alpha t},\tag{30}$$

Equation (29) results in the eigenvalue problem

$$\left[\mathbf{K} - \omega^2 \mathbf{M}\right] \alpha = \mathbf{0},\tag{31}$$

which gives the eigenfrequencies ωi and the eigenvectors α = [α (1) , α (2) , α (3)] T , where α (1) = [α (1) 1 , α (1) 2 , ..., α (1) K+L ] T , α (2) = [α (2) 1 , α (2) 2 , ..., α (2) K+L ] T , α (3) = [α (3) 1 , α (3) 2 , ..., α (3) K+2L ] T .

TABLE 4 | Eigenfrequency parameters of the shell in Example 3: case (b).

The elements of these vectors are the three sets of coefficients corresponding to the functions u, ν, and w, respectively. Subsequently, the mode shapes are obtained by substituting α = [α (1) , α (2) , α (3)] T in Equations (17).

The accuracy of the approximation (9) depends heavily on c (see Equations 20-21-23). This was also verified in the problem at hand. Thus we come across to the problem of selecting a "good" value for c, that is, a value of the shape parameters that produces results of acceptable accuracy. Several methods have been suggested (Hardy, 1971; Franke, 1982; Foley, 1994; Rippa, 1999; Katsikadelis, 2009) for selecting a good value for c in 2D problems. Katsikadelis (2006, 2008b) proposed the minimization of the functional (total potential) for obtaining

an optimal value for c. For the present problem, the optimal value is obtained by the search method as the value of c which yields the minimum value of the eigenfrequencies ωi . It was observed that the optimum c giving the minimum 1st eigenfrequency differs negligibly from that yielding the higher minimum eigenfrequencies. Therefore the same optimum value of c can be used to avoid the search method for higher eigenfrequencies.

#### NUMERICAL EXAMPLES

On the basis of the above analysis, a Fortran program has been written. The expressions of the derivatives involved in Equations (9) to (11) and Equations (12) to (14) have been obtained using the symbolic language MAPLE. Though the method applies to cylindrical shell of variable radius of curvature, for reasons of simplicity, the efficiency and accuracy of the developed method are demonstrated by studying the forced and free vibrations of circular cylindrical panels, (**Figure 3**), under different sets of boundary conditions. The NASTRAN FEM code and a model with 400 rectangular elements (**Figure 4B**) are used to compare the results. In all examples the employed material constants are: E = 21 × 106 kN/m2 , ν = 0.30. The results have been obtained running the relevant programs on an Intel Core 1.6 GHz with RAM 4 GB computer.

#### Example 1

We study the dimensionless eigenfrequency parameter f = Rω p (1 − ν 2 )ρ/E of a simply supported circular cylindrical shell panel with movable curved edges in the axial direction (Nx = v = w = Mx = 0 along the curved edges and u = v = w = Ms = 0 along the straight edges). The first six eigenfrequency parameters from MAEM are given in **Table 1** and are compared with a FEM solution. The 1st and 2nd vibration modes for the normal displacement are shown in **Figure 5**, respectively. The numerical results have been obtained with the parameters L = 80, K = 361 and δ = 1.e − 6 as shown in **Figure 4A**. The optimal value copt = 0.06 corresponds to the 1st mode and has also been used for the other five modes. This value is close to that obtained by the formula c = 2/ √ (K + 2L) = 0.088 proposed in Ferreira et al. (2006b). The employed geometrical data are: h = 0.10 m, R = 10.00 m, l = 5.00 m, a = 5.00 m. The CPU time for the FEM solution was 10 s, while for the MAEM was 85 s. Note that the employed code has not been optimized.

#### Example 2

In this example the same shell as in Example 1 is analyzed under the following boundary conditions: Nx = ν = w = θx = 0 along the curved edges and u = ν = w = θs = 0 along the straight edges. The same collocation points as in the Example 1 have been used. The first six eigenfrequency parameters are shown in **Table 2** as compared again with those obtained from a FEM solution. The 1st and 2nd vibration modes for the normal displacement are shown in **Figure 6** respectively. The value copt = 0.06 was employed to obtain results for the first six modes. The CPU time for the FEM solution was 10 s, while for the MAEM was 88 s.

#### Example 3

In this example, a cylindrical shell panel with geometrical data R = 9.896 m, l = 4.949 m, a = 5.00 m is analyzed. Two cases with regard to the thickness have been considered: (a) 2R sin(φ/2)/h = 100 and (b)2R sin(φ/2)/h = 20 . All edges are clamped, i.e., u = ν = w = θx = 0 along the curved edges and u = ν = w = θs = 0 along the straight edges. The numerical results have been obtained with L = 80, K = 361 randomly distributed as shown in **Figure 7**, that is using an irregular distribution, and δ = 1.e − 6. In both cases, the search method resulted copt = 0.05. The first eight eigenfrequency parameters are shown in **Tables 3**, **4** as compared with those obtained from an analytical solution (Lim and Liew, 1994). The 1st and 2nd vibration modes for the normal displacement are shown in **Figure 8** for case (a) and in **Figure 9** for case (b), respectively. These figures show that that the vibration modes are influenced by the thickness of the shell, which is verified in (Webster, 1968; Lim and Liew, 1994).

#### Example 4

In this example, the forced vibrations of a simply supported circular cylindrical shell panel (u = ν = w = Mx = 0

#### REFERENCES

Beskos, D. E. (1991). "Static and dynamic analysis of shells," in Boundary Element Analysis of Plates and Shells, ed D. E. Beskos (Berlin: Springer-Verlag), 93–140. along the curved edges and u = ν = w = Ms = 0 along the straight edges) with zero initial conditions (u(x, 0) = u˙(x, 0) = ν(x, 0) = ˙ν(x, 0) = w(x, 0) = ˙w(x, 0) = 0) has been studied. Its geometrical data are those of Example 1. The applied load is the normal pressure given by qz = sin(t) kN/m2 . The mass density is ρ = 2.446kNm−4 sec2 . The numerical results have been obtained with L = 80, K = 361 distributed as shown in **Figure 4A** and δ = 1.e−6. The employed optimal value is copt = 0.06. The time history of the normal displacement w, the bending moment Ms and the membrane force Ns at the center of the shell are shown in **Figure 10** as compared with those obtained by a FEM solution. The CPU time for the FEM solution was 40 s, while for the MAEM was 325 s.

### CONCLUSIONS

The Meshless Analog Equation Method, a truly meshless method, has been applied to the dynamic analysis of thin cylindrical shell panels in the present study. MAEM is based on the principle of the analog equation, converting the original equations into three substitute equations, two Poisson's and one biharmonic, which are solved using a meshless method. The use of integrated MQ-RBFs to approximate the fictitious sources allows the approximation of the sought solutions by new RBFs, which approximate both the solution and its derivatives accurately. This way the strong formulation of the problem avoids the drawbacks inherent in the conventional MQ-RBFs, while maintaining all the advantages of a truly mess-free method. A method is presented to obtain optimum values for the shape parameter, eliminating the uncertainty in its choice. It was observed that the optimum value of the shape parameter for the 1st mode differs negligibly from those of higher modes and therefore the same value can be used to obtain the eigenfrequencies of higher modes. The solution algorithm is straightforward and quite reasonably easy to program. The numerical examples presented demonstrate the efficiency and accuracy of the proposed method and show that MAEM can be used as an efficient solver for challenging problems in engineering analysis.

# AUTHOR CONTRIBUTIONS

AY has implemented the MAEM (Meshless Analog Equation Method) for the solution of problems governing by fourth order differential equations (plates-shells etc.). JK has invented the AEM (Analog Equation Method), which is a general method for the solution of various engineering problems.

Cheng, A. H. D., Golbeg, M. A., Kansa, E. J., and Zammito, G. (2003). Exponential1 convergence and h-c multiquadric collocation method for partial differential equations. Numerical Methods Partial Differential Equations 19, 571–594. doi: 10.1002/num.10062

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer AC declared a shared affiliation, with no collaboration, with the authors to the handling editor at time of review.

Copyright © 2018 Yiotis and Katsikadelis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Numerical Investigation of Bolted Hybrid Steel-Timber Connections

Themistoklis Tsalkatidis 1 \*, Youcef Amara2 , Samuel Embaye2 and Einar Nathan1

1 Faculty of Science and Technology, Norwegian University of Life Sciences, Ås, Norway, 2 UnionConsult, Oslo, Norway