Strojniški vestnik Journal of Mechanical Engineering no. 7-8 year 2016 volume 62 Strojniški vestnik - Journal of Mechanical Engineering (SV-JME) Aim and Scope The international journal publishes original and (mini)review articles covering the concepts of materials science, mechanics, kinematics, thermodynamics, energy and environment, mechatronics and robotics, fluid mechanics, tribology, cybernetics, industrial engineering and structural analysis. The journal follows new trends and progress proven practice in the mechanical engineering and also in the closely related sciences as are electrical, civil and process engineering, medicine, microbiology, ecology, agriculture, transport systems, aviation, and others, thus creating a unique forum for interdisciplinary or multidisciplinary dialogue. The international conferences selected papers are welcome for publishing as a special issue of SV-JME with invited co-editor(s). Editor in Chief Vincenc Butala University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Technical Editor Pika Škraba University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Founding Editor Bojan Kraut University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Editorial Office University of Ljubljana, Faculty of Mechanical Engineering SV-JME, Aškerčeva 6, SI-1000 Ljubljana, Slovenia Phone: 386 (0)1 4771 137 Fax: 386 (0)1 2518 567 info@sv-jme.eu, http://www.sv-jme.eu Print: Grafex, d.o.o., printed in 310 copies Founders and Publishers University of Ljubljana, Faculty of Mechanical Engineering, Slovenia University of Maribor, Faculty of Mechanical Engineering, Slovenia Association of Mechanical Engineers of Slovenia Chamber of Commerce and Industry of Slovenia, Metal Processing Industry Association President of Publishing Council Branko Širok University of Ljubljana, Faculty of Mechanical Engineering, Slovenia International Editorial Board Kamil Arslan, Karabuk University, Turkey Hafiz Muhammad Ali, University of Engineering and Technology, Pakistan Josep M. Bergada, Politechnical University of Catalonia, Spain Anton Bergant, Litostroj Power, Slovenia Miha Boltežar, UL, Faculty of Mechanical Engineering, Slovenia Franci Čuš, UM, Faculty of Mechanical Engineering, Slovenia Anselmo Eduardo Diniz, State University of Campinas, Brazil Igor Emri, UL, Faculty of Mechanical Engineering, Slovenia Imre Felde, Obuda University, Faculty of Informatics, Hungary Janez Grum, UL, Faculty of Mechanical Engineering, Slovenia Imre Horvath, Delft University of Technology, The Netherlands Aleš Hribernik, UM, Faculty of Mechanical Engineering, Slovenia Soichi Ibaraki, Kyoto University, Department of Micro Eng., Japan Julius Kaplunov, Brunel University, West London, UK Iyas Khader, Fraunhofer Institute for Mechanics of Materials, Germany Jernej Klemenc, UL, Faculty of Mechanical Engineering, Slovenia Milan Kljajin, J.J. Strossmayer University of Osijek, Croatia Peter Krajnik, Chalmers University of Technology, Sweden Janez Kušar, UL, Faculty of Mechanical Engineering, Slovenia Gorazd Lojen, UM, Faculty of Mechanical Engineering, Slovenia Thomas Lübben, University of Bremen, Germany Janez Možina, UL, Faculty of Mechanical Engineering, Slovenia George K. Nikas, KADMOS Engineering, UK José L. Ocana, Technical University of Madrid, Spain Miroslav Plančak, University of Novi Sad, Serbia Vladimir Popović, University of Belgrade, Faculty of Mech. Eng., Serbia Franci Pušavec, UL, Faculty of Mechanical Engineering, Slovenia Bernd Sauer, University of Kaiserlautern, Germany Rudolph J. Scavuzzo, University of Akron, USA Arkady Voloshin, Lehigh University, Bethlehem, USA Vice-President of Publishing Council Jože Balič University of Maribor, Faculty of Mechanical Engineering, Slovenia Cover: The image presents an experimental setup for performing accelerated vibration fatigue tests for nonlinear rivet joints under base excitation with random signal. This relatively simple experimental setup along with extensive numerical post-analysis enables the assessment of fatigue parameters of blind-rivet joint under bending load due to the dynamic response of the riveted plate. The methodology behind this experiment introduces physically consistent evaluation of the product failures during standard and widely performed vibration tests with electro-dynamic shakers. Image Courtesy: Laboratory for dynamics of machines and structures, Faculty of Mechanical Engineering, University of Ljubljana ISSN 0039-2480 General information Strojniški vestnik - Journal of Mechanical Engineering is published in 11 issues per year (July and August is a double issue). Institutional prices include print & online access: institutional subscription price and foreign subscription €100,00 (the price of a single issue is €10,00); general public subscription and student subscription €50,00 (the price of a single issue is €5,00). Prices are exclusive of tax. Delivery is included in the price. The recipient is responsible for paying any import duties or taxes. Legal title passes to the customer on dispatch by our distributor. Single issues from current and recent volumes are available at the current single-issue price. To order the journal, please complete the form on our website. For submissions, subscriptions and all other information please visit: http://en.sv-jme.eu/. You can advertise on the inner and outer side of the back cover of the journal. The authors of the published papers are invited to send photos or pictures with short explanation for cover content. We would like to thank the reviewers who have taken part in the peerreview process. © 2016 Strojniški vestnik - Journal of Mechanical Engineering. All rights reserved. SV-JME is indexed / abstracted in: SCI-Expanded, Compendex, Inspec, ProQuest-CSA, SCOPUS, TEMA. The list of the remaining bases, in which SV-JME is indexed, is available on the website. The journal is subsidized by Slovenian Research Agency. Strojniški vestnik - Journal of Mechanical Engineering is available on http://www.sv-jme.eu, where you access also to papers' supplements, such as simulations, etc. Contents Strojniški vestnik - Journal of Mechanical Engineering volume 62, (2016), number 7-8 Ljubljana, July-August 2016 ISSN 0039-2480 Published monthly Editorial 401 Papers Xiaoning Liu, Gengkai Hu: Elastic Metamaterials Making Use of Chirality: A Review 403 Dezhi Wang, John E. Mottershead: Measurement Precision and Spatial Resolution with Kriging Digital Image Correlation 419 Yang Liu, Sheikh Islam, Ekaterina Pavlovskaia, Marian Wiercigroch: Optimization of the Vibro-Impact Capsule System 430 Ioannis P. Mitseas, Ioannis A. Kougioumtzoglou, Pol D. Spanos, Michael Beer: Nonlinear MDOF system Survival Probability Determination Subject to Evolutionary Stochastic Excitation 440 Fabian M. Gruber, Daniel J. Rixen: Evaluation of Substructure Reduction Techniques with Fixed and Free Interfaces 452 Lionel Manin, Stefan Braun, Damien Hugues: Modelling the Belt - Envelope Interactions during the Postal Mail Conveying by a Sorting Machine 463 Martin Česnik, Janko Slavič, Miha Boltežar: Assessment of the Fatigue Parameters from Random Vibration Testing: Application to a Rivet Joint 471 Editorial This special issue of the Journal of Mechanical Engineering is devoted to the International Conference on Engineering Vibration (ICoEV), held in September 2015 in Ljubljana, Slovenia. This successful conference attracted over 200 participants, coming from more than 30 different countries. The Steering and Organising Committees of the conference made great efforts to maintain the high scientific standard of the papers. The conference was based on three and half days of technical presentations, structured into four parallel sessions. On each day of the conference there was an invited lecture in the morning. The oral presentations were organized into 17 mini-symposia and eight other general topics. Based on the conference contributions a selection of 7 extended research manuscripts are published in this special issue. This selection of manuscripts presents the state of the art in the field of engineering vibration and dynamics. A brief summary of each article is presented below. In their review paper, Xiaoning Liu and Gengkai Hu from the School of Aerospace Engineering, Beijing Institute of Technology, China, discuss some interesting properties of metamaterials in Elastic metamaterials making use of chirality: a review. They show that the coupling between the rotation and the bulk deformations of two-dimensional chiral solids is able to provide a unique resonant mechanism when designing elastic metamaterials, as well as presenting recent advances in this area. Dezhi Wang from Institute for Risk and Uncertainty and John Mottershead from the Centre for Engineering Dynamics, University of Liverpool, UK, analyse the performance of a new global Digital Image Correlation (DIC) approach known as Kriging DIC, which is assessed by a comparison with the classical subset-based DIC through a standard evaluation procedure in their paper entitled Measurement Precision and Spatial Resolution with Kriging DIC. The third paper, from the field of vibro-impact dynamics, Optimization of the vibro-impact capsule system, is written by Yang Liu and Sheikh Islam from Robert Gordon University, School of Engineering, UK, and Ekaterina Pavlovskaia and Marian Wiercigroch from the University of Aberdeen, Centre for Applied Dynamics Research, UK. The paper is focused on the choice of the excitation parameters and the shape of the capsule. The fourth paper is entitled Nonlinear MDOF structural system survival probability determination subject to evolutionary stochastic excitation. The group of co-authors is Ioannis P. Mitseas and Michael Beer from the Faculty of Civil Engineering and Geodetic Science, Leibniz University Hannover, Germany; Ioannis A. Kougioumtzoglou from the Department of Civil Engineering and Engineering Mechanics, Columbia University; and Pol D. Spanos, L.B. Ryon Chair in Engineering, Rice University. In their paper an approximate analytical technique for determining the survival probability and first-passage probability density function (PDF) of nonlinear multi-degree-of-freedom (MDOF) structural systems subject to a non-stationary stochastic excitation vector is developed. The fifth paper is co-authored by Fabian M. Gruber and Daniel J. Rixen, both coming from the Technical University of Munich, Institute of Applied Mechanics, Germany, and is entitled Evaluation of substructure reduction techniques with fixed and free interfaces. In this contribution they evaluate the primal (classical) formulation of the Craig-Bampton method, the MacNeal method, the Rubin method and the dual formulation of the Craig-Bampton method. The sixth paper is written by Manin Lionel and Braun Stefan from the University of Lyon, INSA-Lyon, LaMCoS, France, and Hugues Damien from Solystic S.A, France, and discusses the prediction of the trajectory of an envelope when conveyed between two flat belts in a mail-sorting machine. The title is Modelling the belt - envelope interactions during the postal mail conveying by a sorting machine. The last, seventh, paper is a contribution from Martin Česnik, Janko Slavič and Miha Boltežar, all from the University of Ljubljana, Faculty of Mechanical Engineering, Slovenia. It is entitled Assessment of the Fatigue Parameters from Random Vibration Testing: Application to a Rivet Joint and discusses some fatigue issues in the field of vibrational fatigue, an expanding area in recent years. It introduces a new fatigue-parameter assessment method based on random vibration loading and its application to a blind-hole rivet joint that diminishes the need for additional fatigue tests. The variety of the selected papers shows the broadness of the research topics as well as the application of vibrational aspects in modern engineering. Miha Boltežar Janko Slavič guest editors Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 403-418 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3799 Review Scientific Paper Elastic Metamaterials Making Use of Chirality: A Review Xiaoning Liu - Gengkai Hu Key Laboratory of Dynamics and Control of Flight Vehicle, Ministry of Education, School of Aerospace Engineering, Beijing Institute of Technology, China Metamaterials are artificially designed composite materials exhibiting unusual physical properties not easily found in nature. In most cases, these properties appear in a wave environment where local resonance comes into play. Therefore the design principle of these metamaterials relies intimately on creation of appropriate local resonances and their interplay with background waves. In this review, we show that the coupling between rotation and bulk deformations of two-dimensional chiral solids is able to provide a unique resonant mechanism in designing elastic metamaterials, and the recent advances in this area will be reviewed. We begin with a metamaterial with a single-negative parameter by integrating a chiral lattice with resonating inclusions, and demonstrate that this metamaterial not only is suitable for broadband vibration isolation but also provides a mixed-type resonance due to microstructure chirality. This mixed-type resonance is further explored to realize elastic metamaterial with double-negative parameters, which can refract elastic wave with a negative refraction angle. Finally, we present also recent development on micropolar constitutive models, which are potentially suitable for modeling chiral elastic metamaterials. Keywords: elastic metamaterials, chiral solids, negative refraction, micropolar theory Highlights • Recent advances in design of elastic metamaterials making use of chirality are reviewed. • The bulk-rotational coupled mechanism of two dimensional chiral solids provides a distinct principle realizing doubly negative elastic metamaterial which is feasible to be engineered. • The chiral lattice integrated with resonators is also a highly designable vibration isolating candidate. • Chirality demands more sophisticated continuum theory, e.g. micropolar theory, in order to fully characterize its wave behavior. 0 INTRODUCTION Metamaterials are manmade composites with exotic properties not found in nature through deliberate microstructure design. There have been intensive research activities on these materials during the past decade. The breakthrough on metamaterials should be attributed to Pendry and his co-workers, they designed respectively electromagnetic (EM) composites with metallic wires and split-ring resonators [1] and [2] which can exhibit effectively negative permittivity and permeability in certain frequency ranges. In addition, EM metamaterials with simultaneously negative material parameters are termed as left handed materials (LHMs) and they are featured by the property that a wave traveling through such material will display anti-parallel phase and group-velocity directions. For LHM, a number of exotic phenomena such as negative refraction, reversed Doppler-effect and Cherenkov radiation [3] and [4] are demonstrated theoretically or experimentally. Since metamaterials significantly enlarge material space available in designing wave-control devices, extensive investigations were carried out in the corresponding resonant mechanism, dynamic homogenization theory and microstructural design. The rapid advance of metamaterials was also boosted by the advent of transformation theory [5] and [6], which maps a deformed space into stringent distribution of material property and provides a superior methodology in designing wave devices. This and the metamaterial technique lead to a number of potential applications, perhaps the most creative examples are the invisible cloaking and the super lensing [7] and [8]. Similar ideas are naturally extended to conceive mechanical (acoustic/elastic) metamaterials. The first elastic metamaterial (EMM) was proposed by embedding periodically silicon rubber coated lead spheres into a polymer matrix [9]. The resulting sample of the crystal can prohibit sound transmission in very low frequency range for which the corresponding wave length is over the crystal constant by two order of magnitudes. The underlying mechanism is obviously not due to the Bragg's scattering of traditional phononic crystals, but should be attributed to effective negative mass density (NMD) [10]. To understand the corresponding mechanism, Milton and Willis [11] constructed an illustrative mass-spring model, and showed that the NMD is due to the out-of-phase motion between the observable and hidden parts of microstructure near the resonating frequency. Yao et al. [12] experimentally examined the model and confirmed NMD at the bandgap region of a finite periodic system composed of mass-spring *Corr. Author's Address: 5 South Zhongguancun Street, Haidian District, Beijing, China. liuxn@bit.edu.cn 403 units. A Drude-model of NMD, for which the effective density is all negative below a cut-off frequency, is also found by letting the hidden mass infinitely large (unable to move). Based on this model, Yao et al. [13] designed a subwavelength structure with excellent sound insulating efficiency in very low frequency range. Yang et al. [13] proposed another scheme of NMD by tuning the out-of-plane resonant pattern of a strained membrane attached with small weights, and the NMD is defined by the ratio between the averaged acceleration and restoring force normal to the membrane. Membrane-type EMMs show a great advantage in sound insulation and absorption with lightweight structures [14] and [15]. Milton and Willis [11] also proved that dynamic effective density can also be a tensor by introducing different resonances along perpendicular directions. Huang and Sun [16] modeled the lattice system as an equivalent 2D elastic solid with its effective mass density characterized as a second order tensor. A continuum EMM with anisotropic effective dynamic mass density tensor was also proposed by placing lead cylinders coated with elliptical shaped rubbers in a matrix material [17] and [18]. Anisotropic mass density can also be achieved in non-resonant acoustic metamaterials containing fluid components [19] and [20]. On the other hand, Fang et al. [21] designed an acoustic metamaterial by using a one-dimensional array of subwavelength Helmholtz resonators, and found the abnormal transmission behavior in certain frequency range which can be explained by the effective negative modulus of the acoustic fluid. However, synthesizing the known mechanisms achieving exotic elastic properties to build a bulk EMM with combined NMD, negative bulk modulus (NBM), negative shear modulus (NSM) or desired anisotropy is more complex compared to the EM media, this is mainly due to the inherent coupling between longitudinal and shear wave modes and wave mode conversion. In order to understand the internal mechanism of bulk EMMs, a systematic approach is to study the wave behavior of matrix-inclusion systems with coated sphere or cylinders, for which the analytical scattering solution is available [22] and [23]. It is now clear that the NBN, NMD and NSM are essentially related to the excited local resonant modes of monopolar, dipolar and quadrupolar types, respectively. In this regard, realization of acoustic metamaterial with doubly or triply negative parameters needs the simultaneous activation of more than one types of resonance within a certain frequency band, this either requires mixed fluid and solid components or complicated microstructure. Ding et al. [24] proposed a double-negative acoustic metamaterial by combining an array of alternating bubble-contained water spheres and rubber-coated gold spheres in an epoxy matrix, each composite spheres are carefully designed to have an overlapped frequency of monopole resonance for bubble-contained water spheres and of dipole resonance for rubber-coated gold spheres. Other similar design schemes were proposed [25] and [26] as well. Wu et al. [27] proposed a type EMM by placing rubber coated water beads into a foam matrix, which can possess simultaneous NMD and NSM. Lai et al. [28] designed a hybrid solid made of four types of solid materials, which can selectively trigger the aforementioned resonances with wave direction. However, due to the fluid-solid combination and complex integrity, these EMM designs are theoretically successful but practically remain difficult to be realized for experimental demonstration. By using two-dimensional (2D) network of Helmholtz resonators, Zhang et al. experimentally demonstrated that the acoustic fluid media in the network can possess necessary double negativity or anisotropy for subwavelength imaging [29] or invisible cloaking [30]. Many other EMM based devices have been suggested for promising applications with elastic wave manipulations such as: elastic metamaterials waveguides, elastic wave super lensing, wave cloaking and vibration insulating [31] to [34]. Recently, it is evidenced that a new design principle different from the mentioned resonant scheme can be pursued if microstructure chirality is introduced. In this review, we will summarize the related works. Chirality was first termed by Lord Kelvin [35]: I call any geometrical figure, or group of points, 'chiral', and say that it has chirality if its image in a plane mirror, ideally realized, cannot be brought to coincide with itself. Lots of examples of chirality can be found in nature or manmade objects, such as twisted ropes and chiral nanotubes for three-dimensional (3D) case, triskelion patterns for 2D case, respectively. It should be mentioned that the chirality have been employed earlier in realizing metamaterials to generate negative refraction for EM [36] and for acoustic cases as well [37], which results from the breaking of degeneracy between two circular polarized waves. The mechanism explored in this paper for EMM is however different. As will be demonstrated in the following, an outstanding feature of 2D chiral solids is the coupling between bulk deformation and local rotation. The feature enables the excitation of mixed translational and rotational resonances in the specially designed microstructure, and can be employed to realize double-negative EMMs in a simple and practical manner. The paper is organized in five sections including this introduction. In Section 2, single-negative metamaterials are presented by integrating 2D chiral lattice with local resonators, their applications in broadband vibration isolating are also studied; in Section 3, the mechanism of mixed resonance and simultaneous NMD and NBM is first clarified by discrete models, and the design of realizable continuum versions of double-negative EMMs is given; in Section 4, recent developments on micropolar constitutive model which is potentially more suitable for characterizing chiral elastic materials and chiral EMMs are reported. Finally conclusions are provided in Section 5. 1 SINGLE-NEGATIVE CHIRAL EMM 2D periodic chiral lattice is not only a lightweight material with unique property for static load bearing, but also a highly designable phononic crystal with flexibly topological adaptation [38] and [39]. In particular, tri-chiral lattices, originally proposed by Prall and Lakes [40] to achieve negative Poisson's ratio, are received intensive investigations as ideal candidates for so-called auxetic materials. The geometry of tri-chiral lattice is depicted in Fig. 1 a), where a hexagonal unit cell is highlighted. The lattice is defined by circles of equal radius r linked by straight ligaments of equal length L. The ligaments are required to be tangential to the circles and the angle between adjacent ligaments is equally 60 degrees. The distance between circle centers is denoted as R, while the angle between the line connecting the circle centers and the ligaments is defined as ß. The thickness of the lattice ligaments is denoted as tb. The lattice vectors e can be written in the orthogonal Cartesian vector basis (ib i2) as +1, )R/2, (i, + i2 |R/2. (1) The ratio cosß=L/R is denoted as the topology parameter. By tuning continuously this parameter, a variety of distinct configurations, form traditional triangular lattice to packed circles, can be obtained as shown in Fig. 1b. In order to achieve low-frequency wave attenuation, Liu et al. [41] proposed an EMM model by integrating a tri-chiral lattice with softly coated inclusions, which function as NMD resonators, in the hollow circles. The unit cell of the EMM lattice is shown in Fig. 1c, and the radius of the core cylinder is identified as rc. The first Brillouin zone, the reciprocal lattice vectors (b1, b2) and the irreducible Brillouin zone (IBZ) of the lattice are shown in Fig. 1d. To study the working mechanism of the proposed singlenegative chiral EMM, band structures are calculated by using finite element based Bloch wave analysis. The configuration of the analyzed latticed composite, including geometry and material parameters, are detailed in Table 1. In the figure, results are presented in the normalized frequency Q = rn/rn0 with rn0 being the first order flexural frequency of a simply supported ligament of length L. Fig. 1. a) Geometry of the tri-chiral lattice; b) different lattice configurations induced by variation of the topology parameter; c) unit cell of the EMM lattice; d) reciprocal lattice vectors and Brillouin zone; taken from Liu et al. [41] Table 1. Geometric and material parameters of the single negative chiral EMM lattice topology parameter L/R = 0.9 ligament length L = 26.4 mm Lattice parameters node radius r = 6.4 mm ligament wall thickness tb = 0.5 mm Young's modulus El = 71 GPa Poisson's ratio vl = 0.33 density Pl = 2.7 g/m3 core-node radius ratio rc /r = 0.5 core Young's modulus Ec = 17 GPa Lattice parameters core Poisson's ratio vc = 0.33 core density Pc = 13 g/m3 coating Young's modulus Es = 5 GPa coating Poisson's ratio vs = 0.33 coating density Ps = 0.5g/m3 The band structure of the chiral EMM is shown in Fig. 2a. For reference, the band diagram of a pure lattice structure without coated inclusions is also plotted in Fig. 2b. It is noticed that the band structures in the two lattice materials with and without resonators are almost unchanged at high frequency range Qe[2, 6], while at low-frequency range, a considerably wide band gap Q e [0.81, 1.43] is found for the former. Dispersion curves of the first two lowest modes, which correspond to P (longitudinal) and S (transverse) modes of the pure lattice, are also plotted in Fig. 2a by dashed lines. It is seen that the lowest two branches are split into five branches due to the presence of local resonators. attenuation is the translational resonance of the core in the soft coating layer and the deformation of the lattice is small. In this case, most part of the wave energy is trapped in the coated inclusion due to its local resonance and wave propagation is not allowed. For the locus away from the edge of the band gap, point O in the first branch corresponds to the rigid mode of the structure, while the modes of point A and B in the fourth branch show the propagating wave takes place through the bending of the lattice ligaments. O A B Ж f - + .Щ iW Fig. 3. Mode shapes of three typical branches at high-symmetry points of the Brillouin zone; taken from Liu et al. [41] Fig. 2. Band diagrams of the lattice a) with and b) without coated inclusions; taken from Liu et al. [41] (For interpretation of the colors in figures, the reader is referred to the web version of this article) Wave modes are selectively presented in Fig. 3 to further understand the formation of the low-frequency band gap. In the figure, mode shapes of the first, third and fourth branches located at the high-symmetry points (O, A, B) of the IBZ, which are highlighted in Fig. 2a by dots, are displayed. Un-deformed geometry (red line) is imposed as reference. The first and fourth branches form the boundaries of the first bandgap. The points A and B in the first branch and the point O in the fourth branch are on the edge of the band gap, for which it can be seen that the mechanism of the wave The third branch in Fig. 2a is remarkable. It is found that the slope of this very narrow passing band Q e[1.14, 1.15] along the OA path is negative, which implies that the anti-parallel directions of group and phase velocity. Moreover, this negative band takes place in the low frequency range where the effective mass density of the EMM is negative, thus it is more appropriate to reason that in this band the effective modulus turns negative as well. The wave modes in Fig. 3 indicate that the formation of this branch is closely associated to the rotational motion of the core cylinder. It is also worthwhile to point out that as long as the core possessed a rotational inertia, the third rotation related branch always exists in the band diagram. However, if the chirality is not present, the dispersion curve of this branch will be strictly horizontal, which implies that the group velocity is zero and the rotational mode actually cannot be excited. Analysis of rotation related band is placed in Section 2.1. In order to test vibration alleviation effect of the previous chiral EMM at low frequency, Zhu et al. [42] experimentally investigated the frequency response of a finite sized beam structure. The sandwiched beam consists of a periodic tri-chiral lattice sandwiched into a rectangular frame (470 mm in length, 91 mm in height, 10 mm in width and 0.5 mm in wall thickness), which is fabricated from an aluminum (Al) plate through a water jet cutter. Local resonators made of rubber (Polyteks Poly PT Flex 20 RTV Liquid Rubber) coated metal cylinders are inserted in lattice beam. Steel and tungsten cylinders with the same geometry are used as inclusion cores for the purpose of generating different resonant frequencies. The EMM parameters are specified by Table 2. A specimen containing two sections of different resonators, which is fabricated to test the broadband suppression of vibration, is shown in Fig. 4a. In the experiment, the beam is fixed on one end and excited by a shaker close to the fixed end. White noise excitation signal with bandwidth from 0 to 1000 Hz is generated by the shaker, and the response of the EMM beam is received by an accelerometer attached to the other end. Band structure analysis is conducted as well for comparison. Table 2. Parameters of fabricated chiral EMM lattice with metamaterial resonators topology parameter L/R =0.82 ligament length L = 24.6 mm Lattice parameters node radius r = 8.6 mm ligament wall thickness tb = 0.5 mm Young's modulus El = 71 GPa Poisson's ratio vl = 0.33 density Pl = 2.7 g/m3 diameter of metal cylinder 6.35 mm height of metal cylinder 25.4 mm Lattice density of steel 7.85 g/m3 parameters density of tungsten 15.63 g/m3 rubber Young's modulus 586 Ma loss tangent of rubber < 0.1 Fig. 4b shows the measured FRF of EMM beam with mixed resonators in solid line. For comparison, FRFs from EMM beam with single section of resonators (steel and tungsten cylinders) are also tested and plotted in dashed and dotted lines, respectively. It is found that the frequency region of attenuation of the EMM beam with mixed resonators falls in between 210 Hz and 700 Hz, which is very close to the summation of those regions of the beams with individual resonator section. In the figure, yellow and gray shaded areas indicate the band gap prediction for infinite chiral EMMs with two types of the resonators, respectively. It is seen that the measured FRFs can be well predicted by either individual band gaps or their overlap correspondingly. 11 Ш Ж m Srtliun with Tungsten cylinders I Section with Sltcl «lindert 0 100 200 300 400 500 600 700 800 900 1000 Frequency III /1 Fig. 4. a) Fabricated EMM beam structure with two sections of resonators; b) comparison of measured and predicted FRFs for the broadband EMM beam; taken from Zhu et al. [42] 2 DOUBLE-NEGATIVE CHIRAL EMM It is seen that introducing chirality in EMMs can effectively activate the usually silent rotational mode of resonance in microstructure. This inspires a new mechanism to design EMMs with simultaneously negative NMD and NBM by making use of coupled translational and rotational resonances through appropriate microstructure design. The resulting double-negative EMM will contain only a single type of resonator, and will be considerably simpler than the traditional schemes. 2.1 Discrete Models In order to illustrate how the negative effective modulus can be produced by rotational resonance, consider a 1D chiral mass-spring unit shown in Fig. 5. Three massless springs and a rigid disk with rotational inertia I are pin-connected. The two springs with elastic constant k2 are tangential to a rigid disk with an angle a. During loading process, the pin-joints A, C as well as the disk center are kept in the horizontal axis. Fig. 5. 1D discrete model demonstrating negative stiffness; taken from Liu et al. [43] The rotation of the disk is induced by the equally applied force F. The tension in the spring k2 is given by: f2 = k2(Rd- x cosa). The rotation of the disk is governed by: dt2 2fR = (2) (3) with R being the radius of the rigid disk. Assuming the system is time-harmonically loaded and all quantities have the form as (F,x,Q) = (F,x,6)еш , the rotation of the disk and the end displacement can be related as в = 2Rk2cosa ; 2R2k2 -a21' (4) The balance of the force at the joint A then gives: F = 2kj X - f2cosa. (5) The dynamic effective stiffness of the system, defined by kfff = F / (2x), reads: 2 ( , , k2 cos a kefr = k +—- eff 1 2 1 - w 2 2 2 W„ -w (6) system consists of an infinite host chain of masses ml connected by springs kl, and the resonators, which consist of mass m2 and rotation inertia I and two inclined springs k2 , are inserted in between the adjacent host masses. A dashed box in the figure highlights the nth unit cell in the system. Different from the previous negative stiffness model, here the disks will undergo not only rotation but also translation. Wang [44] considered a similar ID discrete model in which two chiral resonators with different handedness are superimposed in a unit cell to eliminate the overall non-symmetric effect. The dynamic properties and wave behavior are thoroughly examined in their work. Li and Wang [45] later generalized the discrete model to 2D case. Fig. 6. 1D discrete model demonstrating simultaneous NMD and stiffness where ю„ = y]2k2 R2 /1 represents the resonant frequency of the disk. It can be shown that the effective stiffness takes negative value in the following frequency range: J2k1/(2k1 + k2 cos2 a) < — < 1. —n (7) The physics of the negative kf can be explained as follows. Consider the system subjected to compressive deformation, that is, x is positive. Then Eqs. (4) and (7) give that: Assuming the harmonic response as before, the dynamic equation of the system can be characterized by the following equations: -w'mX = kx (x1(n+1) + xf-1) - 2x(n) ) +k2Rcosa((n) -0(n4)) +k2 cos2 a( + x(2-() - 2x(n) ), -w2m2x(n = k2 cos2 a (x((n+1) + x((n) - 2x2n) ), -a2I0(n) = -Rk2cosa(x(n+1) -x(n))-2R2k20(n). (9) (2kj + k2 cos2a) k2R cosa X , (15) where q is the Bloch wave number, L is the length of the unit cell and (x 1, x2, в) are complex amplitudes. By substituting Eq. (15) into Eq. (9) and solving the eigenvalue problem of the coefficients, the dispersion relation can be expressed as: r \ a 1 -a>2 / a2 ( aa 2 2 «j -a k2 cos2 a k - a0 a0 - a sin2 ^, (16) 1.5 1.0 .3 0.5 0.0 a) b) 1 i • • Ì 4-4 0 4o qL Fig. 7. a) Effective mass, b) stiffness and c) dispersion curves of the doubly negative 1D chiral discrete model Thus, the dispersion curves are plotted according to Eq. (16) in Fig. 7c with the parameters as following: L = 0.1, R = 0.08, a = n/6, mx = 0.2, m2 = 0.3, I= 0.0015, kx = 0.2, and k2 = 0.1. It can be found in Fig. 7c that a pass band with negative slope (red solid line) within frequency range (0.88 to 1.00) is generated inside the bandgap region and separates the bandgap region into two small gaps at (0.75 to 0.88) and (1.00 to 1.18). Fig. 7a shows the normalized effective mass (mfmj) as function of frequency predicted by Eq. (14), the negative effective mass is found at frequency range (0.85 to 1.18), highlighted with grey color in the figure. Fig. 7b shows the normalized effective stiffness calculated by Eq. (13), the negative stiffness is found at frequency range (0.90 to 1.00). It is noticed that the frequency range where negative effective mass and stiffness occurs simultaneously is almost overlapped with that of the pass band with negative slope. On the other hand, the separated two small gaps agree well with the regions where only negative mass exists. 2.2 Three-Phase Chiral EMM with Double-Negative Property Enlightened by the discrete model of Fig. 6, Liu et al. [43] proposed a 2D continuum version of EMM with simultaneously NMD and NBM. The analogy between 1D discrete model and continuum EMM is illustrated in Fig. 8a. The host lattice of discrete model is replaced by continuum matrix material with cavities where the resonators is inserted, while the chiral coating mimics the inclined springs of resonators. To make the EMM macroscopically isotropic, the resonators are arranged in a periodic triangular pattern with lattice constant a. The unit cell of the metamaterial is depicted in Fig. 8b. A number of (ns) slots with width ts are cut out from the coating material. The slots are equi-spaced in azimuth and oriented at an angle 6s with respect to the radial direction. The metamaterial lacks any planes of mirror symmetry hence it is said to be Fig. 8. a) Analogy between the 1D discrete model and 2D continuum EMM; b) unit cell configuration of chiral EMM; taken from Liu et al. [43] m 2 2 chiral. The geometrical parameters mentioned above as well as the constituent materials were carefully designed to ensure the translational and rotational resonances occur at overlapped frequency range, and the overlapped frequency was further optimized to be as large as possible. In order to evaluate dynamic effective properties of the EMM with such a complicated microstructure, analytical method is not available. Instead, a numerical-based effective medium method based on the micromechanics approach is adopted. Under long-wavelength approximation, the macroscopic stress, strain, resultant force and acceleration of the unit cell can be determined by averaging local quantities on the unit cell's external boundary as [46]: 25 Yaß= V J ^ ayXßdSY > Eaß= TV J ((ß + U ßds a ) ' dV 2V dV Ь "Ir Fa= V J aaßdSß > Ua = S J Uads> (17) where the Einstein's summation rule on the repeated subscripts is assumed and Greek subscripts range from 1 to 2; o^ , ua and üa are the local stress, displacement and acceleration fields, respectively; dsa=nads with na and ds being the boundary unit normal and line element of the boundary, respectively; xa and V denote the position vector and unit cell's volume, respectively. Considering the macroscopic isotropy, the effective bulk, shear modulus and momentum mass density of the EMM can be defined as: Keff 2 ^aa 1 ' ^eff 2 ^ aß 1 E aß' P ff = Fa 1 Ua (18) where and Есф denote the deviatoric parts of the macroscopic stress and strain, respectively. To evaluate the dynamic properties, the unit cell is time-harmonically loaded on the boundary by the prescribed displacement: ua(x, t ) = u a (x)e"° , Ua(x) = Ua + EaßXß, (19) where the boundary displacement amplitude u a is required to be compatible with a known macro-strain amplitude E aß . The dynamics of the unit cell is then solved by harmonic finite element analysis with frequency being swept over the interested range, from which the effective properties are determined by Eq. (18). 20 N X es is о §10 D" i V (?) -5 0 5 Density 0 Modulus 2 Г Fig. 9. a) Effective density, b) effective moduli and c) band structure of the three-material chiral elastic metamaterial; taken from Liu et al. [43] Fig. 10. Deformation and rotational resonance modes corresponding to three typical values of the effective dynamic bulk modulus: a) quasi-static value, b) negative value and c) positive peak; taken from Liu et al. [43] An example of the final EMM design is specified with the following parameters: epoxy resin (density pm = 1110 kg/m3, bulk modulus Km = 3.14 GPa, shear modulus pm = 0.89 GPa), low-density Polyethylene (pc=920 kg/m3, Kc=0.57 GPa, pc=0.13 GPa) and lead p = 11600 kg/m3, K, = 52.6 GPa, p = 14.9 GPa) are chosen as matrix, coating and core materials, respectively; the core is 5.6 mm in diameter and the coating thickness is 0.7 mm; the triangular lattice constant is a = 10.75 mm; the slot parameters are ns = 12, ts=0.4 mm and 0s = 56°, respectively. Fig. 9a shows normalized effective mass density peff/pm of the EMM as function of wave frequency. It is seen that the negative density appears in range of 9.51 kHz to 21.54 kHz. Fig. 9b shows the normalized effective bulk modulus Keff/pm and shear modulus peff/pm as a function of frequency. It is of interest to note that KeJf becomes negative in frequency range of 14.08 kHz to 14.72 kHz while pf is always positive. The wave dispersion relation along ГК direction is shown in Fig. 9c. The lattice array and the IBZ are also plotted as the inset of Fig. 9c. A stop band is observed for both longitudinal and transverse waves in the frequency range of 9.44 kHz to 21.58 kHz, which matches the frequency range of negative effective density quite well. In addition, a new pass band with negative slope appears in the frequency range of 14.05 kHz to 14.73 kHz, which again matches the overlapped frequency range of the negative effective density and longitudinal modulus (Kef+ pef) well. The origin of the NBM of the proposed EMM can be understood by examining the deformation and traction of a unit cell corresponding to three typical values of Kef as shown in Fig. 10. The un-deformed states are identified by dashed lines. In Fig. 10a (0 kHz), a quasi-static Kef is obtained since the frequency is far from the core rotational resonance frequency. When the frequency approaches the rotational resonance frequency from below (14.5 kHz), a very large clockwise (in-phase) rotation of the core is generated in conjunction with expansion of the unit cell. Such a rotation produces a compressive state in the matrix and compression is also detected on the external boundary of the unit cell. Conversely in Fig. 10c, when the frequency approaches the rotational resonance frequency from above (15.2 kHz), the anticlockwise rotation of the core enhances the tensile state in the matrix and, consequently, a positive peak of Keff occurs. 2.3 Single-Phase Chiral EMM with Double-Negative Property The chirality and coupled translational and rotational resonances can be employed to design doublenegative EMMs through a single type of unit cell. This could simplify significantly the design of EMMs, however the microstructural pattern given in the previous section is still too complicated for experimental validation. To this end, Zhu et al. [47] proposed an innovative chiral EMM design made of only one single solid material, which is easy to be fabricated and tested in plate-based techniques. The unit cell of the proposed EMM is shown in Fig. 11a, which is composed of one matrix material and slotted voids. The chirality of the unit cell is formed by properly slot-cutting in a hexagonal area and leaving three inclined ribs which connect the center piece (functioned as the mass) and the frame. The lattice constant is a, and the widths of the slot, ribs and frames are denoted by 5, r and f, respectively. In the design, the three ribs are very critical since they function as the soft chiral coating introduced in the previous example, and they will support not only the translation resonance but also the rotation resonance of the center piece. Therefore, by carefully choosing the base material and optimizing the geometrical parameters, it is possible to simultaneously achieve translation and rotation resonances and in turn the double negativity at a desired frequency range. A practical example of the proposed singlephase EMM is specified as following: stainless steel with density 7850 kg/m3, Young's modulus 200 GPa, Poisson's ratio 0.3 is chosen as the base material; the size of unit cell and other geometric parameter are a = 12 mm, 5 = 0.5 mm, r = 0.4 mm and f = 1.2 mm, respectively. The effective bulk modulus and the effective mass density of the chiral EMM are evaluated by the method explained in the previous section, and plotted as a function of frequency in Figs. 11c and d, respectively. The effective mass density peff and effective bulk modulus Keff are normalized with the density and Young's modulus of the base material, respectively. The dispersion curves along ГМ direction are also graphed in Fig. 11b, in which the lattice array and its IBZ are also inserted as the inset of Fig. 11b. It is found that the bandgap frequency range (grey area in 37.2 kHz to 53.6 kHz) predicted from the dispersion relation is almost the same as the frequency region (pink area in 37.4 kHz to 54.1 kHz) of the NMD. In Fig. 11c, Keff turns negative in a frequency range of 42.9 kHz to 45.2 kHz (blue strip). Eventually, both the effective longitudinal modulus Ef= Keff +peff and peff become negative in the regime of 43.6 kHz to 45.2 kHz, which implies a pass band with negative Fig. 11. a) The unit cell of the single-phase chiral EMM; b) band structure; c) effective bulk modulus and d) effective density of the proposed chiral elastic metamaterial; taken from Zhu et al. [47] slope for longitudinal wave in this frequency range. The longitudinal wavelength belonging to this band is from 112 mm to 121 mm, and it is much larger than the unit cell size a = 12 mm. A prominent feature of double-negative metamaterial is the negative refraction phenomenon when a propagating wave impinges the interface between the metamaterial and a normal media, originated from the antiparallel direction of energy flow and phase velocity as well as conservation of transverse component of the wave vector. The double negativity of the proposed chiral EMM was numerically and experimentally justified by the negative refraction between an EMM and a normal elastic material. In the simulation, a 30-degree wedged sample composed of 512 EMM unit cells is constructed in a surrounding normal elastic material which is the same as the base material of the EMM, see Fig. 12b. A longitudinal wave beam is generated from the left side edge of the EMM wedge by applying a harmonic displacement excitation. The wave beam propagates inside the metamaterial, and eventually refracts at the inclined edge. The simulation is performed in COMSOL Multi Physics. Fig. 12. a) The EFCs of the EMM superimposed with the EFCs for transverse and longitudinal waves in stainless steel at 43.8 kHz; b) divergence and c) curl of the velocity fields atfc = 43.8 kHz; taken from Zhu et al. [47] To precisely explore refracted wave modes and the way of the incident longitudinal wave propagation through the EMM, the equi-frequency curves (EFCs) of the metamaterial in the frequency range of 39 kHz to 44.6 kHz in the first Brillouin zone are plotted in Fig. 12a. The wave refraction is examined at the frequency fc = 43.8 kHz, in which both the NMD and NBM are present, and the EFC for this frequency is highlighted by red solid curve in the contour. In addition, the EFCs for the P and S wave in the host stainless steel at the same frequency are also plotted in the figure as the blue dotted circle and green dashed circle, respectively. In the figure, Vg and k are the group velocity and wave vector of the incidental wave in the metamaterial, respectively, and Vgt is the group velocity of S wave in the stainless steel. The refracted wave can be determined in the EFCs by using the Snell's law, which states that the components of the wave vector parallel to the refracting interface have to be the same for the incident and refracted waves. Thus the wave vector of the refracted wave can be determined by drawing a line (black dashed line) passing the end point of incident wave vector and being perpendicular to the interface, and then seeking the intersections of the line with the EFCs of the host media. It is seen from the figure that the intersection occurs in the EFC of S waved of the host material, and a refractive wave with angle ф=-37° is predicted, while no P wave is refracted since any intersection with blue circle is impossible. A wave mode conversion from P to S wave is found. It should be noted that the red EFC at this frequency does not possess exactly a circular shape, hence the vectors of the group velocity and phase velocity are not exactly antiparallel. In order to confirm these founding, full wave simulation is conducted, the divergence and curl of the velocity field obtained from finite element solution are shown in Fig. 12b and c to distinguish the P and S wave contents, respectively. It is seen that very weak refracted longitudinal wave field can be observed and much stronger refracted transverse wave field is found, as expected. The refraction angle estimated from the simulation, as shown in Fig. 12c, also agrees well with the prediction based on EFC. The negative refraction was validated by experiment. Due to 2D feature of the problem, the plate-based wave testing technique was adopted. The wedge-shape EMM array (326 mm x 192 mm X 1.5 mm) was fabricated in a thin stainless steel (Grade 304) plate (3048 mm x 1829 mm x 1.5 mm) by using a precision laser cutting technique, which is shown in Fig. 13. The zoomed-in views of the metamaterial at different scales are also shown in the figure. In the experiment, the dimension of the host plate is chosen to be sufficiently large in order to avoid unwanted reflected waves from the boundary. Since the interested frequency is very low (<50 kHz), the lowest symmetric Lamb wave mode (S0) can be a good approximation for the in-plane wave behavior. An incident P wave is launched by symmetrically bonded piezoelectric patches on top and bottom of the plate. In the test, a steady wave excitation with single frequency is avoided since the absorbing boundary for a solid plate is difficult to realize, instead, a tone burst signal with narrow frequency band is chosen. The refracted in-plane waves are measured on the surface of the plate by a 3D laser scanning vibrometer (Polytec PSV-400-3D). The frequency -wavenumber filtering technique is also applied to remove undesired out-of-plane mode generated due to the inaccuracies in the fabrication and experiment setup. Fig. 13. The fabricated triangular array of the chiral EMM; taken from Zhu et al. [47] Fig. 14. a) The measured velocity field at t =1.7 ms with frequency 43.8 kHz; b) Time domain signal envelops of the scanning points on the three selected lines; taken from Zhu et al. [47] The snapshot of the measured amplitude of the in-plane velocity field for t = 1.7 ms at 43.8 kHz is shown in Fig. 14. A local coordinate system x'-y' is defined with x' parallel to the oblique interface. It can be clearly found that the measured refracted wave propagates downward in the negative refraction side with respective to the interface normal. The measured angle of the refraction is about -38°, which is very close to the numerical prediction. Fig. 14b shows the amplitudes of the recorded wave signal envelopes in the time domain at three regularly spaced lines with a distance from each other by 33 mm along y' direction (see the inset in Fig. 14a) lines 1, 2 and 3). From this measurement the group velocity of the outgoing wave in the host material can be estimated as Vg = 3320 m/s, which is exactly the transverse wave speed of the steel plate as predicted previously. 3 MICROPOLAR MODEL FOR CHARACTERIZING PLANAR CHIRAL LATTICES The effective property and wave characteristic of the chiral elastic metamaterial presented in the previous sections are considered in the framework of classical Cauchy elasticity, and this theory succeeds in material design and the interpretation of the observed wave phenomena. However, chiral material cannot be fully described by Cauchy elasticity since it is not able to characterize the handedness of the material [48]. For example, consider the tri-chiral lattice and its mirror reflection shown in Fig. 15, the difference between the two figures should be reflected somehow by constitutive equations. However, classical elastic theory remains the same upon a transformation of mirror reflection, which is equivalent to the inversion of one coordinate axis. It is interesting to note that, for the tri-chiral lattice, if we adopt the sign convention of topology parameter ß according to the relative orientation of the ligament and the link of circle centers, the sign of ß monitors the handedness of the lattice, as indicated in Fig. 15. This operation cannot be achieved by an in-plane rotation due to its chiral nature. To more comprehensively characterize the 2D chiral solids, one should head to the higher order elastic theory, e.g. micropolar theory. For 3D case, it is found early that isotropic chiral elastic material can be well characterized by the non-centrosymmetric micropolar theory [48]. The recently developed 2D chiral micropolar theory [49] to [51], and some new wave phenomena are reviewed in this section. Fig. 15. Chiral lattice with a) ß>0 and its handedness reversed pattern with b) ß<0 Characterization of material chirality is closely related to the concept of pseudo (or axial) tensors which alternate sign under a mirror reflecting transformation or the handedness change of the underlying coordinate system, while ordinary (or polar) tensors are not affected by such actions. In micropolar theory, rotational degree of freedom (DOF) фi is introduced in addition to the displacement ui at a material point [52]. The deformation measures are characterized as strain and curvature: Jm , Kkl = Фк, (20) respectively, and the balances of stress aJi and couple stress Mji are governed by: и, = рд 2ut / dt2 = Jd 2ф / dt2, (21) where ei]k is the Levi-Civita tensor, p and J are the density and micro-inertia, respectively. Subscripts of Latin letters range from 1 to 3, subscripts of Greek letters range from 1 to 2, and a comma in subscript denotes partial differentiation with respect to coordinates. The governing equations is completed by the constitutive equations: CijklSkl -втф ijklik ,l > mj = Bjkiski -D»A ijklrk ,l > (22) where C, D and B are elastic tensors of rank four. Consider 2D problem defined in the xj-x2 plane where u3 = ф1 = ф2 = д / dx3 = 0, Eqs. (20) and (21) reduce to: ^ßa Ф k =Ф a т,a a ßa,ß=Pd Ua 1 dt Uma,a+ eaßaaß = J д Ф1 dt (23) (24) respectively, where ф3 = ф, ка3 = ка, ma3 = ma are defined for brevity, and eaß=e3aß can be considered as the 2D Levi-Civita tensor. Liu et al. [49] proved that 2D isotropic chiral micropolar material should be characterized by the following constitutive equation which is for clarity written in a matrix form: Стц а22 CT12 а21 m m23. 2ß + X X - A A 0 0 X 2ß + X - A A 0 0 - A - A ß + K ß-K 0 0 A 0 0" "u A 0 0 U2,2 ß-K 0 0 "2,1 -ф ß+K 0 0 ul2 +ф 0 Y 0 Фх 0 0 Y. Ф2 .(25) There are four classical micropolar elastic constants (Lame's constants X and p., antisymmetric shear modulus к and higher order modulus y) and a new parameter A characterizing the chiral effect. When the handedness of the material pattern is reversed, the chiral constant A should reverse its sign, and the other constants remain unchanged. From Eqs. (23) to (25), the wave equations of 2D isotropic chiral micropolar material read: d2u p— = (A + 2M)u,xx + Gu + к)иуу + (X + ^-K)v,iy + dt1 +2кф - A(v - v - 2u - 2ф ), T, y V ,xx , yy ,xy Trs' dv 'dt2 p—22 = (A + 2u)u,xx + (U + K)u,yy + (A + U- K)vxy + +2кфy - A(vxx - v yy - 2uxy - 2ф,х ), д 2ф j W = у(ф-хх +Ф,yy ) - 4кФ + 2K(v-x - Uy ) - -2 A(ux + v y ). (26) The form of Eq. (26) looks like those of anisotropic medium, however, they are essentially different and cannot be covered by any anisotropy without chirality, since the parameter A and its sign form a unique pattern in the constitutive matrix in Eq. (25). The positive definiteness of the strain energy density imposes condition А2<(А+м)к on the chiral constant A, which can be either positive or negative, yet its absolute value is bounded. Fig. 16. Tetra-chiral lattice For the tetra-chiral lattice shown in Fig. 16, obviously a 2D orthotropic chiral micropolar constitutive model is needed, which is recently given by Chen et al. [50] and [51] employing the theory of irreducible orthogonal tensor decomposition. The obtained constitutive tensors display a hierarchy structure depending on symmetry of underlying microstructure. Depending on the microstructure symmetry, up to eight material constants, in addition to the five for the isotropic case, are introduced to characterize a 2D orthotropic chiral elastic material. The constitutive Eq. (25) provides a more sophisticated framework to characterize the tri-chiral lattice introduced in Section 2. By assuming the circle of the tri-chiral lattice is rigid, the five effective micropolar constants of the tri-chiral lattice can be analytically obtained as: ■у/ЗЕ X = —n(cos2 ß -Ц)sec3ßcos2ß, л/ЗЕ /л =—n(cos2 ß+ц2 )sec3ß, к = ^Es n(sin2 ß + n2)secß, Y = - 2 EjL_ 4л/з n (З sin2 ß + 2ц2 )secß, A = ~E~ n(n2 - cos2 ß)secßtan ß, (27) where n=tb/R indicates the slenderness ratio of the ligaments. It can be verified that when the chiral lattice is flipped over (sign of ß is reversed), A changes its sign, while all the other parameter remains unchanged. When ß=0 for a traditional triangular lattice (see Fig. 1b), A vanishes as expected. The chiral micropolar theory can reveal unique wave property in 2D chiral solids which cannot be predicted by traditional elasticity theory. The most pronounced difference between the chiral and non-chiral micropolar media is that for the later a non-dispersive longitudinal wave with velocity cp = [(Л+2,м)/р]1/2 can always be decoupled from the other two shear-rotation coupled waves. This is the characteristic of the non-chiral micropolar media, i.e. the microrotation is only coupled with shear but not with dilatation. In the chiral micropolar theory, the rotation is coupled with the dilation deformation due to the non-zero chiral constant A. Hence there would be no longer pure P or pure S waves in such media. The three wave modes are all mixed and dispersive, thus we call P, S or R (rotation) dominated waves, respectively. Moreover, the common feature of circular polarization for 3D isotropic chiral micropolar material [53] is not presented in the current 2D case, i.e. material particles are linearly polarized during the wave motion. However, the loss of mirror symmetry is reflected in another way for the 2D case. Since the medium is in-plane isotropic, the frequency dispersion has to be isotropic and the EFCs of this medium should be concentric circles. On the other hand, the polarization will remain a fixed angle with respect to the wave vector. This feature is schematically shown in Fig. 17a, where the polarization of mixed P/S wave mode forms a chiral pattern without reflective symmetry. P/S mixed polarization accompanied by isotropic dispersion (EFC) is a unique behavior which cannot be predicted by traditional elasticity theory. Fig. 17b shows the dispersion curves (solid lines) of the tri-chiral lattice predicted by using the homogenized wave equation Eq. (26) in which the effective material constants Eq. (27) are taken. The lattice geometry is specified by R = 1.0, ß=0.9, П = 1/20.The ligaments of the lattice are assumed to be massless, while the mass and rotational inertia of the rigid circles are assumed to be unit. The exact solution of the dispersion curves (circles) obtained from Bloch wave analysis of the tri-chiral lattice is also (a) ' V u./k / ч у \ s ' •*y \ / ì 1 >*'* / 4 >•■' г A X / 4 у f 4 Fig. 17. a) Schematic wave behavior of the 2D isotropic chiral material; comparison of b) dispersion curves and c) polarization angle for the chiral and non-chiral homogenization with the discrete model; taken from Liu et al. [49] plotted in the figure as a benchmark. For comparison, dispersion curves (dashed lines) predicted by a non-chiral version of micropolar homogenization [39] are also plotted in the figure. The dispersion curves are grouped in black, red and blue colors for the first, the second and the third branches, respectively. The first and second branches correspond to the displacement dominated modes, while the third one is the rotational dominated wave. The second branch is almost nondispersive. All three models agree well for this branch. However, for the first and third branches, the chiral theory agrees well with those given by the exact solution of the corresponding discrete model and a large discrepancy is found for the non-chiral theory. It is also interesting to notice that the dispersion curves of P and S dominated wave almost coincide at the long wave limit, indicating almost same phase wave velocities. This is the feature of waves in materials with Poisson's ratio v=-1, where the shear modulus is much greater than the bulk modulus. Since the wave is linearly polarized, it is appropriate to examine the wave mode of the displacement dominated wave through the polarization angle Л. The polarization angle of the first two branches predicted by the chiral and non-chiral theories are plotted in Fig. 17c as function of the wave number, where the first and second branches are marked in black and red, respectively. The non-chiral theory always predicts pure P and S waves, thus the polarization angles remain to be 0 and 90 degrees, as expected. For the chiral theory, the S and P dominated wave (for example Л>60° and Л<30°) can be observed when the wave number is small for the first and second branches, respectively. However, for the intermediate wave number, the P/S domination of the two branches become indistinguishable and can even interchange, i.e. the first branch become P-dominated and the second branch become S-dominated. Good agreement between the chiral micropolar homogenization and the exact discrete model (shown by circle and square dots) is found. More recently, based on the chiral micropolar homogenization of the tri-chiral lattice and tetra-chiral lattice, Bacigalupo and Gambarotta [54], reconsidered the single- negative EMM lattice presented in Section 2. By adopting the matrix chiral lattice as a chiral micropolar continuum and considering the resonator as additional degree of freedoms which interact with each material point through appropriate translational and rotational stiffness, a generalized chiral micropolar continuum model with six degree of freedoms is established for the EMM lattice. The dynamic equations of the proposed model are CT11,1 + CT12,2 + k d (1 - U ) = P1U1 ct211 + ct222 + k d (v2 - u2 ) = pU2 , (28) m11 + m22 + CT 21 - ct12 + k e (9 - ф) = ^ф kd (U1 - V1 ) = P2V1 2—( n (n — l) and the root mean squared error (RMSE) is given by, (7) RMS =J — a? +££. (8) The unidirectional planar sinusoid imposed in on the rectangular speckle pattern defines the true displacement as, = A sin I 2n x u = 0 (9) where A and p denote the amplitude and period. Examples of the imposed sinusoid and measured displacement are shown in Fig. 1 for subset-based and Kriging DIC. The sine wave has an amplitude 5 pixels, a period of 75 pixels and zero-mean Gaussian noise is added with a standard deviation 3 grey values. The deformed speckle pattern [8] and [11] is shown in Fig. 2. The red dashed lines shown in Fig. 1 are constructed normal to the measurement at the Kriging control points. This case illustrates one candidate measure of the discrepancy {£,-}= given by the difference between the measurement and the intersection of the normal with the sinusoid. Another possibility would be to use the discrepancy measured only in the direction of the displacement (i.e. the Y-direction discrepancy), but in this case one would expect the discrepancy to be greater than that measured in the direction of the normal. The smaller the measure of discrepancy the better, and therefore normal measure is considered to be superior to the Y-direction measure. Of course, since the sinusoid is by definition a nonlinear function of x, it is necessary to iterate to determine the intersection precisely. The displacement measurement precision and spatial resolution were defined in general terms in the Introduction. In the results presented here the spatial resolution is understood to be one half of the lowest spatial period in pixels at which the DIC measurement is able to reproduce the true displacement to within a RMSB of 5 % or 1.5 % of the amplitude of the sinusoid. The precision is defined as the standard deviation ae (corresponding to the 5 % or 1.5 % discrepancy on RMSs), measured in pixels in either the normal direction or the Y-direction. 3 PERFORMANCE OF KRIGING DIC ALGORITHM In this section, performance of the Kriging-DIC algorithm, in terms of the displacement precision and spatial resolution, is assessed by means of a comparison with the classical subset-based DIC algorithm. A series of images with different sinusoidal deformation fields are used to calculate the displacement precision and spatial resolution of both the subset-based DIC and the Kriging DIC. The sinusoidal displacement field cannot be perfectly reconstructed using polynomial shape functions and therefore provides a good test and a standard procedure for assessing the performance of DIC algorithms. In this study, the original displacement measurement is investigated, rather than strain which is a derived quantity obtained by post-processing. The reference images are defined by an experimental speckle pattern while the parameters of the deformed images with sinusoidal displacement fields are shown in Table 1. A region of interest with uniformly distributed sample points (centres of subsets) is selected which contains several periods of the sinusoidal displacement. Both the subset-based DIC and Kriging DIC methods are applied to measure the displacements at the sample points and then calculate the Y-direction discrepancy and the normal discrepancy respectively. In both cases the displacement precision is quantified in terms of standard deviation ae. The DIC algorithm parameters are given in Table 2 and the normalized sum of n а. ю d) - i кч \ : J : ; , i ; f1 H- / I v 1 \ . Ili, /Iii \ у \ ; i i^^n i V 'i i i — jy^ff \ \ i i i i ! 1 1 \ 1 ^ 1/ 1 / T T j \ i 4 N i ! /, / 1 T T M N M H ! 1 о ф о го Q. -10 .2 200 О 220 240 260 280 300 320 340 X-direction locations on the speckle image (pixel) 360 380 400 260 280 300 320 340 X-direction locations on the speckle image (pixel) Fig. 1. The normal discrepancy: a) subset DIC size 31*31, b) subset DIC size 41*41, c) subset DIC size 51*51, d) subset DIC size 61*61, e) 1st order Kriging DIC Fig. 2. Deformed speckle pattern with imposed sinusoidal displacement field having a period of 75 pixels Fig. 3. The RMSE of the normal discrepancy In terms of the percentage (5 % and 1.5 %) of the sinusoidal amplitude vs the period of sinusoidal deformation Fig. 4. 80 100 120 140 160 Period of the deformation sine wave (pixel) The displacement precision (STD of normal discrepancy) vs the period of sinusoidal deformation squared differences (NSSD) criterion was instead of the classical SSD for reasons of increased robustness [34]. Table 1. Parameters of the sinusoidal deformed images Table 2. Parameters of DIC algorithms Parameter Value Amplitude 5 pixels Period 25 25 200 pixels Gaussian noise (standard deviation) 3 grey values Subset-based DIC Kriging DIC Criterion NSSD NSSD Sample (control) points 31x10 31x10 Subset size 31 —— 61 (pixels) Regression order 0, 1st and 2nd Correlation function Exponential Shape function 2nd-order Intensity interpolation 6x6 bi-cubic 6x6 bi-cubic 40 60 80 100 120 140 160 180 200 Period of the deformation sine wave (pixel) Fig. 5. The RMSE of the Y-direction discrepancy in terms of the percentage (5 % and 1.5 %) of the sinusoidal amplitude vs the period of sinusoidal deformation 40 60 80 100 120 140 160 180 200 Period of the deformation sine wave (pixel) Fig. 6. The displacement precision (STD of Y-direction discrepancy) vs the period of sinusoidal deformation Figs. 3 and 4 and Figs. 5 and 6 show the root mean squared error (RMSE) and the standard deviation of the normal discrepancy and the Y-direction discrepancy respectively. As expected, the RMSE and the standard deviation both decrease as the period of the sinusoid increases. In Figs. 3 and 5 the 5 % and 1.5 % errors are shown and it is seen that the Kriging models with 0th, 1st and 2nd order regression polynomials are able to accurately represent sinusoidal displacements with low periods (high spatial frequency). The spatial resolution of a particular DIC code is given by half the period corresponding to the intersection of the RSME curve with the horizontal 5 % or 1.5 % line. Figs. 7 and 8 show the normal-direction measurement precision plotted against the spatial resolution for 5 % and 1.5 % errors. This result is particularly pleasing as it shows the Kriging DIC method to produce significantly better precision for the same spatial resolution as the 31*31 pixel subset DIC code. The equivalent 7-direction comparison is presented in Figs. 9 and 10, where again the Kriging DIC approach produces improved results over the subset DIC code, but not such a great improvements as is apparent from Figs. 7 and 8. The 7-direction measurement tends to overestimate the error, which is better estimated using the normal direction measure, and therefore the result shown in Figs. 7 and 8 is considered to represent a better performance comparison than in Figs. 9 and 10. It is clear also that increasing the order of regression in Kriging DIC enables a small improvement in spatial resolution to be achieved at considerable cost to measurement precision. 20 25 30 35 40 45 50 55 Spatial resolution based on the normal-distance RMSE (one half of the sinusoidal period in pixel) Fig. 7. Displacement precision vs spatial resolution (one half of the periods) based on the normal-discrepancy RMSE under the criterion of 5 % sinusoidal amplitude, for subset-based DIC using different subset sizes and Kriging DIC with different order regression functions 30 35 40 45 50 55 60 65 70 75 Spatial resolution based on the normal-distance RMSE (one half of the sinusoidal period in pixel) Fig. 8. Displacement precision vs spatial resolution (one half of the periods) based on the normal-discrepancy RMSE under the criterion of 1.5 % sinusoidal amplitude, for subset-based DIC using different subset sizes and Kriging DIC with different order regression functions Fig. 9. Displacement precision vs spatial resolution based on the Y-direction RMSE under the criterion of 5 % sinusoidal amplitude, for subset-based DIC and Kriging DIC 40 45 50 55 60 65 70 75 80 85 90 Spatial resolution based on the y-distance RMSE (one half of the sinusoidal period in pixel) Fig. 10. Displacement precision vs spatial resolution based on the Y-direction RMSE under the criterion of 1.5 % sinusoidal amplitude, for subset-based DIC and Kriging DIC 4 CONCLUSIONS The performance of the Kriging DIC algorithm in terms of measurement precision and spatial resolution is considered in this study. A comparison is made between Kriging DIC and the classical subset-based DIC on the basis of a series of speckle images with imposed sinusoidal deformation of various spatial frequencies. Two new measures of displacement discrepancy are used to redefine the measurement precision and spatial resolution for both local and global DIC. The difference between the true sinusoid and the measurement is determined most accurately when measured along the normal at the Kriging control points. This enables a better assessment to be made of the relative merits of different DIC algorithms. Results obtained using the RMSE criterion of 5 % and 1.5 % of the sinusoidal amplitude show Kriging DIC to be superior to classical sub-set based DIC in both displacement measurement precision and spatial resolution. 5 ACKNOWLEDGEMENTS The authors wish to acknowledge the contribution of Lucas Wittevrongel in providing the sinusoidal deformation images used in this study. The first author acknowledges the support of the China Scholarship Council and the University of Liverpool. 6 REFERENCES [1] Aster, R.C., Borchers, B., Thurber, C.H. (2005). Parameter Estimation and Inverse Problems. Academic Press, Amsterdam. [2] Passieux, J.C., Bugarin, F., David, C., Perie, J.N., Robert, L. (2015). Multiscale displacement field measurement using digital image correlation: Application to the identification of elastic properties. Experimental Mechanics, vol. 55, no. 1, p.121-137, DOI:10.1007/s11340-014-9872-4. [3] Sutton, M.A., Orteu, J.J., Schreier, H.W. (2009). Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications. Springer, New York. [4] Rethore, J., Hild, F., Roux, S. (2007). Shear-band capturing using a multiscale extended digital image correlation technique. Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 49-52, p. 5016-5030, D0I:10.1016/j.cma.2007.06.019. [5] Wagne, B., Roux, S., Hild, F. (2002). Spectral approach to displacement evaluation from image analysis. European Physical Journal - Applied Physics, vol. 17, no. 3, p. 247-252, DOI:10.1051/epjap:2002019. [6] Cheng, P., Sutton, M.A., Schreier, H.W., McNeill, S.R. (2002). Full-field speckle pattern image correlation with B-spline deformation function. Experimental Mechanics, vol. 42, no. 3, p. 344-352, DOI:10.1007/BF02410992. [7] Mortazavi, F. (2013). Development of a Global Digital Image Correlation Approach for Fast High-Resolution Displacement Measurements. PhD thesis, École Polytechnique de Montréal, Montreal. [8] Wittevrongel, L., Lava, P., Lomov, S.V., Debruyne, D. (2015). A self adaptive global digital image correlation algorithm. Experimental Mechanics, vol. 55, no. 2, p. 361-378, DOI:10.1007/s11340-014-9946-3. [9] Hild, F., Roux, S. (2012). Comparison of local and global approaches to digital image correlation. Experimental Mechanics, vol. 52, no. 9, p. 1503-1519, DOI:10.1007/ s11340-012-9603-7. [10] Hild, F., Roux, S. (2006). Digital image correlation: from displacement measurement to identification of elastic properties - a review. Strain, vol. 42, no. 2, p. 69-80, DOI:10.1111/j.1475-1305.2006.00258.x. [11] Bornert, M., Brémand, F., Doumalin, P., Dupre, J.-C., Fazzini, M., Grédiac, M., Hild, F., Mistou, S., Molimard, J., Orteu, J.-J.. Robert, L., Surrel, Y., Vacher, P., Wattrisse, B. (2009). Assessment of digital image correlation measurement errors: Methodology and results. Experimental Mechanics, vol. 49, no. 3, p. 353-370, DOI:10.1007/s11340-008-9204-7. [12] Triconnet, K., Derrien, K., Hild, F., Baptiste, D. (2009). Parameter choice for optimized digital image correlation. Optics and Lasers in Engineering, vol. 47, no. 6, p. 728-737, DOI:10.1016/j.optlaseng.2008.10.015. [13] Besnard, G., Hild, F., Roux, S. (2006). "Finite-Element" displacement fields analysis from digital images: Application to portevin-le chatelier bands. Experimental Mechanics, vol. 46, no. 6, p. 789-803, DOI:10.1007/s11340-006-9824-8. [14] Zhou, P., Goodson, K.E. (2001). Subpixel displacement and deformation gradient measurement using digital image/ speckle correlation (DISC). Optical Engineering, vol. 40, no. 8 p. 1613-1620, DOI:10.1117/1.1387992. [15] Perlin, K. (1985). An image synthesizer. Proceedings of the SIGGRAPH Conference, San Francisco, p. 287-296, DOI:10.1145/325165.325247. [16] Orteu, J.-J., Garcia, D., Robert, L., Bugarin, F. (2006). A speckle texture image generator. Proceedings of the Speckle International Conference, NTmes, DOI:10.1117/12.695280. [17] Schreier, H.W., Braasch, J.R., Sutton, M.A. (2000). Systematic errors in digital image correlation caused by intensity interpolation. Optical Engineering, vol. 39, no. 11, p. 29152921, DOI:10.1117/1.1314593. [18] Pan, B., Xie, H.-M., Xu, B.-Q., Dai, F.-L. (2006). Performance of sub-pixel registration algorithms in digital image correlation. Measurement Science and Technology, vol. 17, no. 6, p. 16151621, DOI:10.1088/0957-0233/17/6/045. [19] Bergonnier, S., Hild, F., Roux, S. (2005). Digital image correlation used for mechanical tests on crimped glass wool samples. Journal of Strain Analysis for Engineering Design, vol. 40, no. 2, p. 185-198, DOI:10.1243/030932405X7773. [20] Schreier, H.W., Sutton, M.A. (2002). Systematic errors in digital image correlation due to undermatched subset shape functions. Experimental Mechanics, vol. 42, no. 3,p. 303-310, DOI:10.1007/BF02410987. [21] Réthoré, J., Hild, F., Roux, S. (2008). Extended digital image correlation with crack shape optimization. International Journal for Numerical Methods in Engineering, vol. 73, no. 2, p. 248-272, DOI:10.1002/nme.2070. [22] De Bievre, P. (2009). The 2007 International Vocabulary of Metrology (VIM), JCGM 200:2008 [ISO/IEC Guide 99]: Meeting the need for intercontinentally understood concepts and their associated intercontinentally agreed terms. Clinical Biochemistry, vol. 42, no. 4-5p. 246-248, DOI:10.1016/j. clinbiochem.2008.09.007. [23] Sutton, M.A. (2010). ASTM E2208 Standard Guide for evaluating Non-Contacting Optical Strain Measurement Systems. ASTM International, West Conshohocken, DOI:10.1520/E2208-02R10E01. [24] Lava, P., Cooreman, S., Coppieters, S., De Strycker, M., Debruyne, D. (2009). Assessment of measuring errors in DIC using deformation fields generated by plastic FEA. Optics and Lasers in Engineering, vol. 47, no. 7-8, p. 747-7 53, D0I:10.1016/j.optlaseng.2009.03.007. [25] Canal, L.P., Gonzalez, C., Molina-Aldareguia, J.M., Segurado, J., LLorca, J. (2012). Application of digital image correlation at the microscale in fiber-reinforced composites. Composites Part A:Applied Science and Manufacturing, vol. 43, no. 10, p. 1630-1638, D0I:10.1016/j.compositesa.2011.07.014. [26] Wang, D.Z., DiazDelaO, F.A., Wang, W.Z., Mottershead, J.E. (2015). Full-field digital image correlation with Kriging regression. Optics and Lasers in Engineering, vol. 67, p. 105115, D0I:10.1016/j.optlaseng.2014.11.004. [27] Chen, D.J., Chiang, F.P., Tan, Y.S., Don, H.S. (1993). Digital speckle-displacement measurement using a complex spectrum method. Applied Optics, vol. 32, no. 11, p. 18391849, D0I:10.1364/A0.32.001839. [28] Hild, F., Roux, S. (2012). Comparison of local and global approaches to digital image correlation. Experimental Mechanics, vol. 52, no. 9, p. 1503-1519, DOI:10.1007/ s11340-012-9603-7. [29] Bruck, H.A., McNeill, S.R., Sutton, M.A., Peters, W.H. (1989). Digital image correlation using Newton-Raphson method of partial differential correction. Experimental Mechanics, vol. 29, no. 3, p. 261-267, D0I:10.1007/BF02321405. [30] Vendroux, G., Knauss, W.G. (1998). Submicron deformation field measurements: Part 2. Improved digital image correlation. Experimental Mechanics, vol. 38, no. 2 p. 86-92, D0I:10.1007/BF02321649. [31] Matheron, G. (1993). Principles of geostatistics. Economic Geology, vol. 58, no. 8, p. 1246-1266, D0I:10.2113/ gsecongeo.58.8.1246. [32] Sacks, J., Welch, W.J., Mitchell, T.J., Wynn, H.P. (1989). Design and analysis of computer experiments. Statistical Science, vol. 4, no. 4, p. 409-435, D0I:10.1214/ss/1177012413. [33] Wackernagel, H. (1998). Kriging Weights in Multivariate Geostatistics: An Introduction with Applications. Springer Berlin, Heidelberg, p. 93-99, D0I:10.1007/978-3-662-03550-4. [34] Pan, B., Qian, K.M., Xie, H.M., Asundi, A. (2009). Two-dimensional digital image correlation for in-plane displacement and strain measurement: a review. Measurement Science & Technology, vol. 20, no 6, D0I:10.1088/0957-0233/20/6/062001. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 430-439 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3754 Original Scientific Paper Optimization of the Vibro-Impact Capsule System Yang Liu1* - Sheikh Islam1 - Ekaterina Pavlovskaia2 - Marian Wiercigroch2 1 Robert Gordon University, School of Engineering, UK 2 University of Aberdeen, Centre for Applied Dynamics Research, UK Optimization of the vibro-impact capsule system for the best progression is considered in this paper focusing on the choice of the excitation parameters and the shape of the capsule. Firstly the fastest and the most efficient progressions are obtained through experimental investigations on a novel test bed. Control parameters, the amplitude and the frequency of harmonic excitation, and one of the system parameter, namely the stiffness ratio, are optimized. The experimental results confirm that the control parameters for the fastest progression are not the same as those for the most efficient progression from the energy consumption point of view. Therefore, the capsule system can be controlled either in a speedy mode or in an energy-saving mode depending on the operational requirements. In the second part of the paper, optimization of the capsule shape is studied using computational fluid dynamics (CFD) simulations. Here the aim of achieving the best progression is addressed through minimizing the drag and the lift forces acting on a stationary capsule positioned in the pipe within a fluid flow. The CFD results indicate that both drag and lift forces are dependent on capsule and arc lengths, and finally, an optimized shape of the capsule is obtained. Keywords: capsule system, vibro-impact, experiment, optimization, CFD simulation Highlights • Optimization of the vibro-impact capsule system for the best progression is considered in this paper. • Optimization focuses on the choice of the excitation parameters and the shape of the capsule. • The experimental results confirm that the control parameters for the fastest progression are not the same as those for the most efficient progression from the energy consumption point of view. • Optimization of the capsule shape is studied using computational fluid dynamics (CFD) simulations. • The drag and the lift forces acting on a stationary capsule positioned in the pipe within a fluid flow are minimized. • The CFD results indicate that both drag and lift forces are dependent on capsule and arc lengths. 0 INTRODUCTION The inspection and maintenance of pipelines which are usually called "pigging" pose an ongoing challenge that many industries face, including the oil and gas industry [1] and [2]. The use of in-line intelligent pipeline inspection mechanism (known as "pig") is a common method of inspection to determine the condition of oil pipelines [3]. However, developing such a self-propelled mechanism with the accessibility of pipelines in different diameter sizes [4] and [5] is not an easy task, particularly for the design of its driving mechanisms (e.g. legs). Furthermore, a large proportion of pipelines are installed in remote and hostile environments which make traditional inspection mechanisms inaccessible. The vibro-impact capsule system driven by autogenous internal excitation is a promising solution in which there is a growing interest in recent years, e.g. [6] to [9]. Such capsule mechanisms have been demonstrated as an effective way to move in complicated environments, but much less attention has been paid to the optimization of capsule systems. Hence, this paper aims to address this issue focusing on two aspects, i.e. optimization of the capsule dynamics through the appropriate choice of the excitation parameters, and minimizing the resistant forces acting on the capsule in fluid flow via the determination of capsule's optimum geometric parameters. Optimizations for capsule systems in terms of the average speed of capsule progression have been considered by several researchers. For example, Chernous'ko [10] studied the optimum rectilinear motion of a two-mass system to obtain its maximum mean velocity. Li et al. [11] proposed a minimum energy solution for an internal mass-driven capsule system. In [12], the dynamics of a mobile capsule system containing a movable internal mass was optimized with respect to its maximal average steady-state velocity. However, there is less published work which studied optimizations for such capsule systems experimentally. In [9], using a mathematical model of the capsule system, Liu et al. found that the parameters of the vibro-impact capsule system for the best progression and for the minimum energy consumption were different. To consider the capsule motion in different environmental conditions, vibro-impact responses of the capsule system with various friction models were investigated in [13]. In [14], a position feedback control method suitable for 430 *Corr. Author's Address: Robert Gordon University, School of Engineering, Garthdee Road, AB10 7GJ, Aberdeen, UK, y.liu8@rgu.ac.uk dealing with chaos control and coexisting attractors was studied for enhancing the desirable forward and backward capsule motion. Recently, the dynamical response of the capsule system was studied by means of path-following techniques in [15] focusing on two practical problems which were maximizing the rate of progression and directional control of the system by following a typical period-1 trajectory. In this paper, the experimental verification of the theoretical finding from [9] will be presented, where the recently designed test bed proposed in [16] is used for this experimental study. In [16], the comparisons between the experiments and numerical simulations were presented showing a good agreement, and the conducted bifurcation analysis indicated that the behaviour of the system was mainly periodic and that a fine tuning of the control parameters can significantly improve the performance of the system. Differently, optimal control parameters with respect to capsule average speed and power efficiency will be identified experimentally in this paper through varying the frequency and amplitude of external excitation and stiffness ratio. As the capsule geometry plays a significant role in the dynamic behaviour of the system in the fluid, computational fluid dynamics (CFD) has been adopted by many researchers to investigate capsule movement in pipeline. In [17], Feng et al. studied two-dimensional simulations of the motion of elliptic capsules carried by a Poiseuille flow in a channel. Khalil et al. [18] studied the pressure distribution around a capsule by using a three-dimensional steady-state turbulent flow CFD model. Borregales et al. [19] investigated the pigging process inside pipeline through two-dimensional transient CFD simulations, and found that sliding at constant velocity allows pig to move easily inside a straight pipe. Li et al. performed unsteady turbulent flow CFD simulations of a single coal log by using a dynamic mesh method in [20]. In this paper, CFD modelling will be carried out for optimizing capsule shape (e.g. capsule and arc lengths) in order to obtain the minimum drag and lift forces in a dynamic flow condition. CFD simulations were carried out considering that the capsule moves in the same direction as the fluids flow. So the findings in this paper could be used for prototype design and fabrication. The rest of the paper is organized as follows. In Section 1, the experimental apparatus and the mathematical model of the vibro-impact capsule system are briefly introduced, and then experimental optimization of the capsule dynamics is presented. In Section 2, the CFD analysis is described and the results are discussed. Finally, some conclusions are drawn in Section 3. 1 DYNAMICS AND ENERGY OPTIMIZATIONS 1.1 Experimental Apparatus To verify the theoretical finding obtained in [9], the test bed [16] shown in Fig. 1a is used in this study. The experimental rig consists of a linear DC servomotor mounted on a base frame connected with a standing frame which holds a support spring with an adjustable stiffness k2. The motor has a movable rod with the mass m1 harmonically excited with a desired frequency m and amplitude pd through the electro-magnetic fields generated by the coils within the motor. Although a nonlinear resistance force keeps the rod in place when the motor is switched on, we assume that this force could be linearized around the working point and characterised by constant coefficients kj and c multiplied by the displacement and velocity, respectively. A gap S exists between the rod and the support spring, and the rod contacts with the support spring when their relative displacement is Fig. 1. a) Photograph of the test bed; and b) schematic of the experimental setup [16] larger or equals to the gap. A detailed description of its physical parameter identification can be found in [16]. A schematic of the experimental setup is shown in Fig. 1b, where the absolute displacement of the rod is xb and the absolute displacement of the base frame is x2 which is measured by a linear variable differential transformer (LVDT) displacement transducer. The relative displacement of the rod and the base frame X] X2 is measured via three hall sensors within the motor. The acceleration of the rod is obtained using an accelerometer mounted directly on the rod. 1.2 Non-Dimensional Equations of Motion Similar to the system described in [9], [13] and [14], in this study the vibro-impact capsule is modelled as a two degrees-of-freedom dynamical system depicted in Fig. 2, where a movable internal mass m1 is driven by a harmonic force with amplitude pd and frequency m generated by a linear actuator. The actuator contains a movable part connected to the internal mass and a fixed part mounted on the rigid capsule m2. We simplify the model of the actuator here and represent the interaction between the internal mass and the capsule by using a linear spring with stiffness k and a viscous damper with damping coefficient c. xx and X2 represent the absolute displacements of the internal mass and the capsule, respectively. The internal mass contacts a weightless plate connected to the capsule by a secondary linear spring with stiffness k2 when the relative displacement xj-x2 is larger or equals to the gap S. When the force acting on the capsule exceeds the threshold of the dry friction force f between the capsule and the supporting environmental surface, bidirectional motion of the capsule will occur, and the dynamic friction force fd will be applied to the capsule. As the considered system operates in bidirectional stick-slip phases, a detailed consideration of these phases can be found in [9]. Fig. 2. Physical model of the vibro-impact capsule system In order to re-scale the capsule dynamics and extend the experimental findings to a general vibro-impact driven system, we introduce the following non-dimensional variables: _ k, _ dx■ k, t = m0t, x = -f x, yt = — = —X;, fs dT ®0fs y==4^x, fd=4, dT Wofs fs and parameters: К J mx _ m , m = —, mo а = Pd fs ' z = 8 = ^8, ß fs и 7 m2 2mlm0 where i = 1, 2. Then the equations of motion can be written in a dimensionless form as (see [9] for the detailed derivations): xi = У1, У1 =acos(rnz) + (x2 -xj) + 2Z(y2 -y ) -h1ß (x1 - x2 -8 ), x2 = y2 [h2 (1 - h1 ) + h3 fy ], y2 =[h2 (1 - h1 ) + h3 h ][-fd (X2 - X1 ) -2Z (y2 - yi ) + h1 ß (X1 - X2 - 8)] / 7, (1) , 3) and the auxiliary Heaviside functions hi (i = 1, are given by h = h (xi - x2 -S h2 = h (|(x2 - X1 ) + 2С(У2 - yi h3 = h (|(X2 - Xi ) + 2C(y2 - yi ) -ß(( - X2-S)-i). Based on the analysis conducted in [13], where four various friction models were used to describe the motion of the capsule in different environments, the friction of the experimental test bed in [16] was found to be well approximated by the Coulomb Stribeck model given by: ( Ш fd = 1 + e sign (x2 ), (2) where vs is the non-dimensional Stribeck velocity. 1.3 Experimental Optimization In order to optimise the performance of the system, we introduce the average progression of the capsule system over one period of the external excitation Ta = 2л/ю, Pvg = T [ X2 (Ta )-x2 (O)] , (3) and the power efficiency which is the ratio of the capsule progression to the work done by the external excitation over the period Ta, Eavg - X (Ta )- X (0) j Jacos (ют)у1 T)dT (4) After non-dimensionalizing the experimental results, both performance indices, Pavg and Eavg were calculated in order to determine the optimum control parameters. Optimizations of the frequency of excitation ä for ß = 8.92 and ß= 13.35 are presented in Fig. 3 and 4, respectively. The optimal frequency of the fastest forward progression for ß = 8.92 is ä = 0.29, and the optimal frequency for ß = 13.35 is ä = 0.3. It is worth noting that the optimal progression for ß = 8.92 is chaotic including both forward and backward motions in every period of excitation, while only periodic forward motion exists for the optimal progression for ß = 13.35 as shown in the blow-up window of Fig. 4. Optimization of the amplitude of excitation pd for ß = 13.35 is ä = 0.61, and 8 = 0.34 is presented in Fig. 5 where both forward and backward progressions were recorded. It can be seen from the figure that the best forward progression is achieved at a = 0.27 and the best backward progression is obtained at a = 0.41. This experimental result reveals that the largest amplitude of excitation will not generate the fastest Fig. 3. Time histories of the capsule displacement, x2 obtained experimentally by varying the frequency of excitation ä at у=4.82, Z= 0.16, ß=8.92, a=0.97, and 8 = 0.52 Fig. 4. Time histories of the capsule displacement, x2 obtained experimentally by varying the frequency of excitation m at у=4.82, Z= 0.16, ß=13.35, a=0.97, and 8 = 0.52 Fig. 5. Time histories of the capsule displacement, x2 obtained experimentally by varying the amplitude of excitation а at y=4.82, Z= 0.16, ß=13.35, m = 0.61, and S = 0.34; phase trajectories are presented in blow-up windows Fig. 6. Average progress, Pavg and energy efficiency, Eavg obtained experimentally by varying the frequency of excitation, m at a)ß=8.92, and b) ß=13.35, у=4.82, а=0.97, and S = 0.514 capsule progression, thereby a proper selection of optimal parameters should be considered for different configurations of the system. Fig. 6 presents the performance indices Pavg and Eavg recalculated for the experimental results shown in Figs. 3 and 4. As can be seen from Fig. 6a, the fastest progression is m = 0.29 and the most efficient progression is m = 0.23 for ß = 8.92. As the stiffness ratio ß increases to 13.35, the fastest and the most efficient progressions are recorded at m = 0.31 and m = 0.35, respectively. All the capsule displacements in time histories are displayed in the embedded windows where the fastest progression is indicated by black line, the most efficient progression is shown by grey line and the rest of displacements are denoted by light grey lines. It can be seen that more energy is spent per period during the fastest motion due to spell of backward motion. This causes the fastest progression shown by black line less efficient comparing with the one indicated by grey line. Fig. 7 shows a series of experiments by varying the amplitude of excitation. In order to compare the capsule speed in one direction, we chose a set of control parameters which has forward progression only. It can be seen from the figure that both the fastest and the most efficient progressions were achieved at а = 0.29, so only one capsule progression marked by black line is presented in the embedded window. The experimental results as a function of stiffness ratio are shown in Fig. 8. As k was kept as a constant throughout the experiments, the variation of the stiffness ratio were implemented by changing the length of the support spring, i.e. the stiffness of the support spring k2. As can be observed from the figure, the fastest and the most efficient progressions were achieved at ß = 16.1 and ß = 8.92, respectively. It is also noted that when the stiffness ratio is small, both performance indices are much better than the ones for large stiffness ratio. Fig. 7. Average progress, Pavg and energy efficiency, Eavg obtained experimentally by varying the amplitude of excitation, a at ß=13.35, y = 4.82, ä = 0.277, and 8 = 0.514 Fig. 8. Average progress, Pavg and energy efficiency, Eavg obtained experimentally by varying the stiffness of the support spring, ß aty=4.82, a=0.97 ä = 0.304, and 8 = 0.514 The results of the presented experimental study allow one to obtain the optimal excitation forcing in order to drive the capsule system either in fastest or the most efficient ways. Although the experiments were conducted in air, our mathematical modelling [13] indicates that similar results could be obtained if the capsule is placed on the heavily lubricated surface. However, if the capsule is fully submerged and moving in the fluid, there will be lift and drag forces creating the resistance to the capsule motion. Therefore, in the next section, we will consider optimization of the stationary capsule shape using CFD approach. 2 CFD OPTIMIZATION 2.1 CFD Model A two-dimensional steady-state model of the vibro-impact capsule system inside a pipe was developed in this section in order to investigate the effect of shape of the capsule on the hydrodynamic coefficients. The purpose of the study was to optimize capsule shape for reducing resistance force from the fluid. Here one of the key assumptions was that the test bed used in Section 1 is fully enclosed in an outer capsule. Although in real applications, the capsule will be moving along the pipeline within the fluid flow, in the first approximation, the stationary capsule was considered in this study. Table 1. Model parameters for 2D simulation Model Parameters Symbol Value Capsule length Lc 118 mm Capsule diameter Dc 80 mm Arc length rrhs 40 mm Pipe internal diameter Dp 140 mm Pipe length Lp 2000 mm Pipe/Capsule clearance ^ clear 1 mm Inlet velocity V 0.5 m/s Outlet pressure Po 0 MPa Fig. 9. 2D capsule simulation layout and boundary conditions Design of Experiment (DOE) technique is implemented to explore the effects of multi variable variation on the reduction of drag forces, and to improve the efficiency of the vibro-impact mechanism. Two-dimensional modelling approach was selected to reduce the complexity of the modelling process and costly simulation time. Fig. 9 shows the two-dimensional numerical model with boundary conditions, and the physical dimensions of the computational domain given in Table 1 are based on the dimension of our original prototype. The assumptions made in solving the 2D computational model are given as follows. • The optimization study is performed under steady state condition of the capsule, i.e. the operating conditions and the capsule position inside the pipe do not change with time. • The capsule is operated under isothermal conditions. • A pressure outlet condition is assumed and no slip condition is applied to the external walls. • The flow is turbulent based on the fluid properties and boundary conditions. • A small gap is assumed between the capsule bottom and the inner wall of the pipe in order to reduce the high frictional resistance due to contact. As the flow condition was calculated to be turbulent, the k-e turbulent model [21] was adopted for the study. The model includes the turbulent kinetic energy represented by the turbulent viscosity: Я =PCnк2 /S and the turbulent kinetic energy, dk — + Uj dt ' dx = T, dUi dx, --S- dx, , ч dk (v + ) dx, (5) (6) and the turbulent dissipation expressed as: ds TT ds — + Uj- dt ' dx, = Cs, - Ce2 s+A sl k j dx k dx lv KJAs j Iv KJAsj a ds k dx, (7) The flow domain was created as a rectangular cross-section of the pipeline with the capsule position at approximately halfway between the inlet and the outlet. The left face and the right face of the pipeline were defined as being a velocity inlet and pressure outlet, respectively. The upper face of pipeline, lower pipeline and the capsule were defined as standard walls. error is considered less than 10~3 for convergence. The multi-factor DOE optimization technique was utilized in simulations in order to establish the most suitable design parameters for the capsule. Two parameters, capsule length Lc (between 118 and 150 mm) and arc length Rrhs (between 40 and 80 mm) were considered as the factors for each simulation, and the best design for reduced drag force allowing the capsule to move faster was obtained. These measurement regions were selected based on the design constraints of the vibro-impact mechanism of the capsule. The goal of the optimization is to improve the existing capsule design such that the hydrodynamic coefficients are made as low as possible. Since two design parameters were weighed against each other in each study, three dimensional response output plots of the results were created to show which parameters allow the most optimal operation. The goal-driven optimization has created 50 sample simulations, and 3 potential candidates were chosen by using the objective function. The computational domain has been meshed with 106818 elements including close wall biasing. A grid sensitivity test using over 200000 cells has proven that the grid size is sufficient to provide grid independency. In addition, each simulation took approximately 1800 iterations with an average convergence time of 20 minutes, and velocity contours, streamlines, pressure contours and response surface were achieved through ANSYS Post processing. 2.2 Computational Procedure and Modelling Parameters 2.3 Results and Discussion The governing equations were solved using the finite volume method for CFD solver (ANSYS FLUENT) which is based on the SIMPLE (Semi Implicit Method for Pressure Linked Equation) algorithm. The relative This section outlines and discusses the results for the DOE study which investigates the effects of the capsule length Lc and the arc length RRhs on the drag and the lift coefficients. Fig. 10 presents the response V Fig. 10. Response surfaces for a) drag and b) lift coefficients 436 Liu, Y. - Islam, S. - Pavlovskaia, E. - Wiercigroch, M. Fig. 11. Effect of capsule and arc lengths on a) drag and b) lift coefficients; coefficients are marked by diamonds and linear fits are shown by lines surfaces for both coefficients. It can be seen from the figure that the minimum drag coefficient is achieved when both arc and capsule lengths increase. Fig. 11 shows the parametric correlation between the drag and the lift coefficients with the capsule and the arc lengths. The plots were created using the results obtained from 10 design points out of 50 analysed. The results confirm again that both coefficients are dependent on the arc and the capsule lengths. Fig. 12 presents the pressure contours, the velocity contours, and the streamlines for the original and the optimum capsule designs. Comparing both designs, the downstream nose cone has been enlarged and the capsule length has been increased for the optimum design. Pressure and velocity contours show a slight reduction in passing fluid pressure and velocity at downstream. This small variation reduces the drag coefficient, and hence will improve the motion of the capsule in pipe. The goal-driven optimization studies identified three best candidates for the minimum drag and lift forces, and the optimized design variables for the best candidate of capsule design are given in Table 2. A graphic comparison between the original and the optimum designs is presented in Fig. 13. Fig. 12. Plots of a) pressure contours, b) velocity contours, and c) streamlines for the original (upper) and the optimum (lower) designs Table 2. Optimized capsule measurements Model parameter Symbol Original value Optimum value Capsule length Lc 118 mm 138.16 mm Arc length rrhs 40 mm 79.15 mm Drag Lift Fig. 13. Comparison of the original (darker) and the optimum (lighter) capsule designs 3 CONCLUSIONS Optimization of the vibro-impact capsule system was studied in this paper experimentally using a novel test bed and numerically through CFD simulations. The paper is focused on addressing two optimization issues, namely the choice of the excitation force parameters for driving the capsule and the choice of the capsule shape to minimize the resistance force. The aim in both cases is to obtain the best progression of the capsule system. Optimization for the progression of the capsule system was carried out experimentally by varying the frequency and the amplitude of excitation. The experimental results revealed that the best progression of the system for stiffness ratio ß = 8.92 was achieved at the excitation frequency ä = 0.29 when the motion of the capsule was chaotic including both forward and backward motions in every period of excitation, while only periodic forward motion exists for ß = 13.35 and the optimal frequency is ä = 0.3. Our comparisons on the average progression Pavg and the power efficiency Eavg have confirmed that the control parameters for the fastest progression are not the most efficient ones, so the capsule system can be controlled either in a speedy mode or an energy-saving mode depending on the requirement of operation. The optimum capsule geometric parameters for the best progression were identified through minimizing the drag and the lift forces using CFD simulation by assuming that a stationary "ideal" capsule was positioned in the pipe with fluid flow. Here we assumed that the experimental test bed was encapsulated and only an outer capsule was considered for CFD simulation. An optimum design could be achieved by minimizing the drag and lift resistance forces acting on the capsule. The CFD results indicate that both drag and lift forces are dependent on the capsule and the arc lengths. Finally, an optimum set of capsule and arc lengths was obtained and the improvement was demonstrated using a histogram representation. Parametric studies for this optimum design through three-dimensional modelling and time dependent simulations would be the future work which could provide a comprehensive understanding of how the capsule interacts with the fluids and accurately predicts the performance of the capsule in real scenarios. Consideration of the moving capsule in the fluid flow would be the next step in our modelling for the shape optimization. Also, creating a model where the capsule is fully submerged and moving in the fluid with proper lift and drag forces acting on the capsule would be an interesting study which would allow refining the driving force optimization strategy. The uncertainty of the input parameters, environmental friction and fluid velocity will be investigated in the dynamical model and the CFD analysis in the future work, respectively, and the findings will be reported in due course. 4 ACKNOWLEDGMENTS Dr. Yang Liu would like to acknowledge the financial support for the Small Research Grant (31841) by the Carnegie Trust for the Universities of Scotland. 5 REFERENCES [1] Shukla, A., Karki, H. (2016). Application of robotics in onshore oil and gas industry - a review Part I. Robotics and Autonomous Systems, vol. 75, p. 490-507, D0l:10.1016/j. robot.2015.09.012. [2] Shukla, Д., Karki, H. (2016). Application of robotics in onshore oil and gas industry - a review Part II. Robotics and Autonomous Systems, vol. 75, p. 508-524, D0I:10.1016/j. robot.2015.09.013. [3] Russell, D., Errington, N. (2012). Assessing the effectiveness of pipeline cleaning programs. Pipeline & Gas Journal, vol. 239, no. 6. [4] Zhang, Y., Jiang, S., Zhang, X., Ruan, X., Guo, D. (2011). A variable-diameter capsule robot based on multiple wedge effects. IEEE-ASMETransactions of Mechatronics, vol. 16, no. 2, p. 241-254, DOI:10.1109/TMECH.2009.2039942. [5] Park, J., Hyun, D., Cho, W., Kim, T., Yang, H. (2011). Normalforce control for an in-pipe robot according to the inclination of pipelines. IEEE Transactions on Industrial Electronics, vol. 58, no. 12, p. 5304-5310, D0l:10.1109/TIE.2010.2095392. [6] Nagy, Z., Leine, R.I., Frutiger, D.R., Glocker, C., Nelson, B.J. (2012). Modeling the motion of microrobots on surfaces using nonsmooth multibody dynamics. IEEE Transaction on Robotics, vol. 28, no. 5, p. 1058-1068, D0I:10.1109/ TR0.2012.2199010. [7] Carta, R., Sfakiotakis, M., Pateromichelakis, N., Thoné, J., Tsakiris, D.P., Puers, R. (2011). A multi-coil inductive powering system for an endoscopic capsule with vibratory actuation. Sensors and Actuators A: Physical, vol. 172, no. 1, p. 253258, D0I:10.1016/j.sna.2011.03.036. [8] Sfakiotakis, M., Pateromichelakis, N., Tsakiris, D. (2014). Vibration-induced frictional reduction in miniature intracorporeal robots. IEEE Transactions on Robotics, vol. 30, p. 1210-1221, DOI:10.1109/TRO.2014.2334931. [9] Liu, Y., Wiercigroch, M., Pavlovskaia, E., Yu, H. (2013). Modelling of a vibro-impact capsule system. International Journal of Mechanical Sciences, vol. 66, p. 2-11, D0I:10.1016/j.ijmecsci.2012.09.012. [10] Chernous'ko, F. (2002). The optimum rectilinear motion of a two-mass system. Journal of Applied Mathematics and Mechanics, vol. 66, no. 1, p. 1-7, D0I:10.1016/S0021-8928(02)00002-3. [11] Li, H., Furuta, K., Chernousko, F.L. (2006). Motion generation of the Capsubot using internal force and static friction. Proceedings of the 45th IEEE Conference on Decision and Control, p. 6575-6580, D0I:10.1109/CDC.2006.377472. [12] Fang, H., Xu, J. (2011). Dynamics of a mobile system with an internal acceleration-controlled mass in a resistive medium. Journal of Sound and Vibration, vol. 330, no. 16, p. 40024018, D0I:10.1016/j.jsv.2011.03.010. [13] Liu, Y., Pavlovskaia, E., Hendry, D., Wiercigroch, M. (2013). Vibro-impact responses of capsule system with various friction models. International Journal of Mechanical Sciences, vol. 72, p. 39-54, D0I:10.1016/j.ijmecsci.2013.03.009. [14] Liu, Y., Pavlovskaia, E., Wiercigroch, M, Peng, Z. (2015). Forward and backward motion control of a vibro-impact capsule system. International Journal of Non-Linear Mechanics, vol. 70, p. 30-46, D0I:10.1016/j.ijnonlinmec.2014.10.009. [15] Paez Chavez, J., Liu, Y., Pavlovskaia, E., Wiercigroch, M. (2016). Path-following analysis of the dynamical response of a piecewise-linear capsule system. Communications in Nonlinear Science and Numerical Simulation, vol. 37, p. 102114, D0I:10.1016/j.cnsns.2016.01.009. [16] Liu, Y., Pavlovskaia, E., Wiercigroch, M. (2016). Experimental verification of the vibro-impact capsule model. Nonlinear Dynamics, vol. 83, no. 1, p. 1029-1041, D0I:10.1007/s11071-015-2385-6. [17] Feng, J., Huang, P.Y., Joseph, D.D. (1995). Dynamic simulation of the motion of capsules in pipelines. Journal of Fluid Mechanics, vol. 286, p. 201-227, D0I:10.1017/ S002211209500070X. [18] Khalil, M.F., Kassab, S.Z., Adam, I.G., Samaha, M.A. (2009). Prediction of lift and drag coefficients on stationary capsule in pipeline. The 13th International Water Technology Conference, p. 435-461. [19] Borregales, M.A., Ensalzado, R., Asuaje, M. (2014). Prediction of lift and drag coefficients on stationary capsule in pipeline. ASME 2014 International Mechanical Engineering Congress and Exposition, p. V007T09A082, D0I:10.1115/ IMECE2014-37452. [20] Li, W., Lu, S., Liu, Y., Wang, R., Huang, Q., Yan, J. (2015). CFD Simulation of the unsteady flow of a single coal log in a pipe. The Canadian Journal of Chemical Engineering, vol. 93, no. 11, p. 2084-2093, D0I:10.1002/cjce.22312. [21] Botros, K.K., Golshan, H. (2010). Field validation of a dynamic model for an MFL ILI tool in gas pipelines. The 8th International Pipeline Conference, p. 325-336, D0I:10.1115/IPC2010-31018. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 440-451 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3752 Original Scientific Paper Nonlinear MDOF system Survival Probability Determination Subject to Evolutionary Stochastic Excitation Ioannis P. Mitseas14* - Ioannis A. Kougioumtzoglou2 - Pol D. Spanos3 - Michael Beer145 1 Leibniz University Hannover, Faculty of Civil Engineering and Geodetic Science, Germany 2 Columbia University, Department of Civil Engineering and Engineering Mechanics, USA 3 Rice University, Department of Mechanical Engineering, USA 4 University of Liverpool, Institute for Risk and Uncertainty and School of Engineering, United Kingdom 5 Tongji University, School of Civil Engineering and Shanghai Institute of Disaster Prevention and Relief, China An approximate technique for assessing the reliability of nonlinear multi-degree-of-freedom (MDOF) systems subject to a non-stationary stochastic excitation vector is developed. The proposed technique can be construed as a two-stage approach. First, relying on statistical linearization and utilizing a dimension reduction approach the nonlinear n-degree-of-freedom system is decoupled and cast into (n) effective single-degree-of-freedom (SDOF) linear time-variant (LTV) oscillators. Second, utilizing the effective SDOF LTV oscillator time-varying stiffness and damping elements in conjunction with a stochastic averaging treatment of the problem, the MDOF system survival probability and firstpassage PDF are determined. Overall, the developed technique appears to be efficient and versatile since it can handle readily, at a low computational cost, a wide range of nonlinear/hysteretic behaviors as well as various stochastic excitation forms, even of the fully non-stationary in time and frequency kind. A 3-DOF system exhibiting hysteresis following the Bouc-Wen model is included in the numerical examples section. Comparisons with pertinent Monte Carlo simulations demonstrate the accuracy of the technique. Keywords: first-passage problem, nonlinear stochastic dynamics, evolutionary stochastic processes, nonlinear/hysteretic systems, survival probability Highlights • A stochastic dynamics methodology for determining the survival probability of nonlinear MDOF systems. • Approximate analytical expressions provided for estimating the time-varying survival probability. • Survival probabilities and the associated first-passage PDFs are determined at a low computational cost. • The developed technique is characterized by enhanced versatility as it can handle readily a wide range of nonlinear behaviors as well as various stochastic excitations with arbitrary EPS forms. 0 INTRODUCTION Excitations acting upon dynamical systems such as wind, wave, and seismic loads commonly exhibit evolutionary features. In this setting, not only the intensity of the excitation but also its frequency content exhibit strong variability. This fact necessitates the representation of this class of structural loads by non-stationary stochastic processes. Further, structural systems under severe excitations can exhibit significant nonlinear behavior of the hysteretic kind. Thus, of particular interest to the structural dynamics community is the development of techniques for determining the response and assessing the reliability of nonlinear/hysteretic systems subject to evolutionary stochastic excitations (e.g., [1] to [3]). Further, in engineering dynamics, the evaluation of the probability that the system response stays within prescribed limits for a specified time interval is advantageous for reliability based system design applications. In this regard, the first-passage problem, that is, the determination of the above time-variant probability known as survival probability, has been a persistent challenge in the field of stochastic dynamics for many decades. Monte Carlo simulation techniques are among the most potent tools for assessing the reliability of a system (e.g. [4]). Nevertheless, there are cases where the computational cost of these techniques can be prohibitive, especially when large-scale complex systems are considered; thus, rendering the development of alternative efficient approximate analytical/numerical techniques for addressing the first-passage problem necessary. Indicatively, one of the early approaches, restricted to linear systems, relies on the knowledge of the mean up-crossing rates and on Poisson distribution based approximations (e.g., [5] to [7]). Further attempts to address the firstpassage problem range from analytical ones (e.g., [8]) to numerical ones (e.g., [9]). Furthermore, techniques based on the concepts of the numerical path integral (e.g., [10] to [13]), of the probability density evolution (e.g., [3]), or of stochastic averaging/linearization (e.g., [14]) constitute some of the more recent approaches. 440 *Corr. Author's Address: Institute for Risk and Reliability, Faculty of Civil Eng. and Geodetic Science, Callinstr. 34, 30167, Hannover, Germany, mitseas@irz.uni-hannover.de In this paper, an approximate analytical technique for determining the survival probability and firstpassage probability density function (PDF) of nonlinear multi-degree-of-freedom (MDOF) systems subject to an evolutionary stochastic excitation vector is developed. Specifically, first relying on a statistical linearization based dimension reduction approach the original MDOF system is decoupled and cast into (n) effective single-degree-of-freedom (SDOF) linear time-variant (LTV) systems corresponding to each and every degree of freedom of the original MDOF system. Second, a stochastic averaging based approximate technique is utilized to derive the nonlinear MDOF system survival probability and first-passage PDF at a low computational cost. The remainder of this paper is organized as follows: In section 1.1 the statistical linearization technique for nonlinear MDOF systems is presented. Next, in section 1.2 a stochastic averaging/statistical linearization treatment of the problem, through a system dimension reduction approach is briefly delineated. In section 1.3, it is shown that the nonlinear MDOF system non-stationary marginal, transition and the joint response amplitude probability density functions (PDFs) can be approximated by closed-form expressions. Further, section 2 provides analytical closed-form expressions for the time-dependent survival probability of the nonlinear MDOF structural system as well as for the corresponding first-passage PDF. In section 3, illustrative examples comprising a 3-DOF system exhibiting Bouc-Wen hysteresis and subject to evolutionary stochastic excitations are considered. Pertinent MCS data demonstrate the reliability of the proposed technique. Finally, section 4 provides with concluding remarks. 1 MDOF SYSTEM DIMENSION REDUCTION In this section, the basic elements of an approximate dimension reduction/decoupling technique developed by some of the authors for determining the non-stationary response amplitude PDF of nonlinear MDOF systems subject to evolutionary stochastic excitation are reviewed for completeness; see [15] and [16] for a more detailed presentation. 1.1 Statistical Linearization Treatment Consider an n-degree-of-freedom nonlinear system governed by the equation: where y denotes the response acceleration vector, y is the response velocity vector, y is the response displacement vector, defined in relative coordinates; M, C and K denote the (n * n) mass, damping and stiffness matrices, respectively; g(y, y) is an arbitrary nonlinear (n * 1) vector function of the variables y and y . F(t)T=(/1(t), f2(t), .ft)) is a (n * 1) zero mean, non-stationary stochastic excitation vector process defined as F(t)=yä(t) where YT = (Yi, Y2, •••, Yn) is an arbitrary (n * 1) vector of constant weighting coefficients, and ä (t ) is a non-stationary process with an evolutionary power spectrum (EPS) Sa (cu, t). In this regard, F(t) possesses the EPS matrix: S f (a, t ) Yi2Sa (a, t) 0 0 гХ (a, t ) 0 0 YSa (a, t) .(2) Further, the non-stationary stochastic excitation process is regarded to be a filtered stationary stochastic process according to the concept proposed by Priestley [17]; see also [18]. Thus, the excitation EPS matrix of Eq. (2) takes the form: SF (cu, t ) = A (cu, t (m) A (cu, t ) (3) where the superscripts (T) and (*) denote matrix transposition and complex conjugation, respectively; A(m, t) is the modulating matrix which serves as a time-variant filter; and SF (ет) is the power spectrum matrix corresponding to the stationary stochastic vector process F (t ). Note that both separable and non-separable EPS can be defined considering Eq. (3). In this manner, excitations exhibiting variability in both the intensity and the frequency content can be considered. Focusing next on the frequency domain, the response determination problem is defined as seeking the corresponding system response EPS matrix of the form: S (a>,t) S.1.1 К t ) S.1.2 S.2.1 (®> t ) S.,.1 К t ) (К t ) (К t) S.1., Кt ) >>t ) (К t) (co, t) .(4) My + Cy + Ky + g(y, у) = F (t). (1) According to the statistical linearization method (e.g., [1] to [3]), a linearized version of Eq. (1) takes the form: My + (C + Ceq )y + (k + Keq )y = F(t). (5) where H(m) is the frequency response function (FRF) matrix defined as: Next, adopting the standard assumption that the response processes are Gaussian, the time-dependent elements of the equivalent linear matrices Ceq and Keq are given by the expressions: ceq=e I dy j and keq- = E № - к (6) (7) Further, for a linear MDOF system subject to evolutionary stochastic excitation a matrix input-output spectral relationship of the form: Sy Ю t ) = H gen (ю, t )St (to) Hen (co, t ), (8) can be derived (e.g., [1] and [3]), where H gen (c t ) = \h (t-T)A(co,x)e^dT. (9) 0 In Eq. (9) denotes the impulse response function matrix. Furthermore, the time dependent cross-variance of the response can be evaluated by the expression: E[У'У}] = jSy,y, i0,()dm- (10) It can be readily seen that Eqs. (6) to (10) constitute a coupled nonlinear system of algebraic equations to be solved numerically for the system response covariance matrix. Note in passing that instead of the frequency domain Wiener-Khinchin relationship of Eq. (8), a state-variable formulation can be adopted yielding a system of differential equations of the Lyapunov kind (e.g., [1] and [19]) for the system response covariance matrix. Nevertheless, although a pre-filtering treatment can be applied for considering non-stationary stochastic excitation processes of the separable kind (e.g., [1]), excitations possessing a non-separable EPS (e.g. realistic cases of earthquake excitations) cannot be accounted for, at least in a straightforward manner. Next, omitting the convolution of the impulse response function matrix with the modulating matrix can lead to substantial reduction of computational effort, especially for the case of MDOF systems (e.g., [16] and [20]). In this manner, Eq. (9) takes the form: H (ffl) = (-ffl2 M + m(c + Ceq ) + (( + Ka (12) Consequently, taking into account Eqs. (3) and (8), Eq. (11) becomes: (со,t) = H(сО,t)H(w). (13) Note that the Eq. (13) can be regarded as a quasi-stationary approximate relationship which, in general, yields satisfactory accuracy in cases of relatively stiff systems (e.g., [20] to [22]). Note in passing that the spectral input-output relationship of Eq. (13) is exact for the case of stationary processes (e.g., [1] to [3]). Further, adopting the aforementioned quasi-stationary approach, it can be readily seen that for the ith degree of freedom, using Eqs. (2), (10) and (13) yields: E [ y2 (t )]=J ( H и H2 +... —ю + |Hin (ет)|2 y2n)Sa (co,t)dm, (14) Ю t E [ y? (t )]=_p (j H и (с)2 7l2 +... —к> + |Hn (ет)|2 yi)Sa (со,t)da>, (15) and Eqs. (14) and (15) hold true in the approximate quasi-stationary sense delineated earlier. Clearly, Eq. (13) constitutes an approximate formula for determining the MDOF system response EPS matrix at a low computational cost; thus, circumventing computationally intensive Monte Carlo simulations. 1.2 Effective SDOF Linear Time-Variant System Following next the system dimension reduction approach developed in [16], an auxiliary effective SDOF LTV system corresponding to the ith degree of freedom can be defined as: У +ßeq,i (t ) У ( ) yi =äi ( )' (16) Hgen (ffl, t ) = H (co) A (co, t ), (11) where the time-varying equivalent stiffness and damping elements of the effective LTV system can be determined by equating the variances of the response displacement and velocity expressed utilizing the quasi-stationary FRF of Eq. (16) with the corresponding ones determined via Eqs. (14) and (15); this yields: * [ yi (t )] = = 1 С (t)-®2)2 + ßi (t)C xy.2Sä (со,t)da>,(17) and .K; (t )-®2)2+ß (t H2 xy.2Sä (a>, t )dm, (18) Clearly, Eqs. (17) and (18) in conjunction with Eqs. (14) and (15) constitute a nonlinear system of two algebraic equations to be solved for the evaluation of the LTV system time-varying equivalent stiffness m2qi (t) and damping ßeqA(t) coefficients. Note that determining the time-varying natural frequency meq,i(t) is especially important for a number of reasons such as tracking and avoiding moving resonance phenomena (e.g., [23]), determining peak system response estimates based on design spectrum compatible excitation power spectra (e.g., [24]), or developing efficient approximate techniques for determining nonlinear system survival probability and first-passage PDF (e.g., [14]). Next, a stochastic averaging technique (e.g., [15] and [16]) is applied for casting the second-order stochastic differential equation (SDE) of Eq. (1) into a first-order SDE governing the evolution in time of the response amplitude ai(t). In this regard, and based primarily on the assumption of light damping, it can be argued that the response yi (t) of the effective LTV system of Eq. (16) exhibits a pseudo-harmonic behavior described by the equations: y, (t) = ai (t)cos( (t)t + Vi (t)), (19) and У (t) = -meq,i (t()x sm(,i (()t + Vi (()). (20) In Eq. (19) the response amplitude ai (t) is a slowly varying function with respect to time defined as: ( ■ 1л V a2 (t) = y] (t)- .Ж (t ). ю \ eq, (21) whereas ф (t) stands for the phase of the response yi (t). Further, relying on a combination of deterministic and stochastic averaging (e.g., [16]) a first-order SDE governing each and every degree-of-freedom response amplitude process ai (t) takes the form: *(t ) = --2 (t )a (t)+ 2a ( K,i (t ) + ^SF (,i (t), t) i (t) n(t ). (22) In Eq. (22), n(t) stands for a stationary, zero mean and delta correlated Gaussian white noise process of unit intensity, i.e., E(n(t)) = 0 ; and E(q(t)q(t+T)) = S(t), with S(t) being the Dirac delta function. Associated with the above SDE (Eq. (22)) is the Fokker-Planck (F-P) partial differential equation governing the response amplitude transition PDF of the Markovian process a; that is, д — P (ai,2' h Kl> t1 ) = д да. 1 д2 +--; 2 да,1 (t )■ а, + nSF (,i (t)'t) 2а,<, (t) p (a 2 > t2\ an, t1 J nSF ( (t)'t) <, (t ) p(ai2,t2 \ ai1,t1 ) (23) Further, considering the case p(ai2, t2 | ai,1 = 0, tj = 0) = p(ai, t), the marginal system response amplitude PDF has been shown to follow a time-dependent Rayleigh distribution of the form (e.g., [16], [25] and [26]): t \ a. p (a., t ) = —j— exp c (t) ч 2c (t) (24) where ci(t) accounts for the non-stationary variance of the LTV system of Eq. (16). As it was shown in [15] and [16] substituting Eq. (24) into Eq. (23) and manipulating yields the following nonlinear ordinary differential equation (ODE): ci (t) = ß (t)ci (t)- nSF ((Peg,i (t) t) i (t) (25) to be solved for the non-stationary LTV system response variance ci(t) via standard numerical schemes such as the fourth order Runge-Kutta. 1.3 Transition and Joint Nonlinear System Response PDFs Taking into account that no change of state can occur if the transition time is zero i.e., p(ai2,t1 | au, t1) = S(ai2- ai,1) and following a similar analysis as the one in [25], the transition response amplitude PDF p(ai,2, t2 | ai,1, t1) for the ith degree-of- 1 1 2 freedom of the original MDOF system is assumed to be of the form: p (a, 2, 1 a,i, ti ) = f „2 i ( t2 ) exp -h? (tl, t2 ) 2С (ti,t2) I0 x 2hi (t1, t2 ) 'i (t1, t2 ) ,(26) where cj(/1, t2) and hi(t1, t2) are functions to be determined and I0 represents the modified Bessel function of the first kind and of zero order. Next, substituting Eq. (26) into the F-P Eq. (23) and manipulating (see also [9],[14] and [25]) yields the linear first-order ODEs: dci (tl, t2 ) dt„ ßeqi (c ( t2 ))c ( t2 )- nSF (i (j (tl, t2 )), t2 ) Ki (Ci (t1' t2 )) = 0 (27) and ^^ + 2ßv ( t2)hi ) = 0. (28) Relying on the assumption that the equivalent damping and stiffness coefficients follow a slowly varying with respect to time behavior, the following approximations over a small time interval [tiJ_1, tij\ are introduced; i.e., ßeq,l(th])=ßeqß—i) and Veq,i(tij) = ®eq,Äj-i) for t e [ty-b j Next, based on the slowly varying with time behavior of the EPS, Sf,, (m,t) is also treated as a constant over the interval [ti,j-1, ti,j\. Further, based on the above assumptions, introducing the variable тi,j-=t,- - ti;-1, and applying a first-order Taylor expansion around the point т,-=0, Eqs. (27) and (28) become (see [14] for a detailed derivation): , ((-1, h j ) = nSF ( (-1 ) 'i,j-1 )) T t ) ' meq,i ^'i,j-1 ) and hi ( (j-l, t,j ) = a.j-^1 -ßeq,i (j )j . (29) (30) Furthermore, considering Eqs. (25) and (29) and applying a first-order Taylor expansion for the response variance ci(t) around the point t=ti j-1 yields: ^ (j h, j ) = C. fo )-C. (j )(!-ßeqi {'U-! )) .(3D Relying next on the Markovian assumption for the process ai, the joint-response amplitude PDF p(ai,j-1, tj-1 ; ai,j, t-) is given by: P (j-^ 0-1 ; ai,j > tj ) = P ( ai j-i> h-i) P(ai,j ' 0 1 ai,j-i' ji ).(32) Utilizing Eqs. (24) and (26), Eq. (32) becomes: P («( O-i ; a, j, t j ) = - a, j-a ,j exp c j h (-v t j ) ( ~atjci ((j-i ) - ah-ici ((j-i, t j ) _ 2c, ((-i h ((-l,(j ) h.2((jh((j-i)1r (ah((j-l, j 2c,- ((j-i h (j ) Further, setting , ((j-l, (j ) (33) 2 ri.j = g ( ( (Uj Eq. (31) yields: С (,.ti,j ) = ci (,j )(l - r,$ ). (34) (35) Next, considering Eqs. (29) and (30) and Eqs. (34) and (35), the joint response amplitude PDF p(ai, j-1, ti, ;-1 ; ai -, ti, j) of Eq. (33) is given in the form: p (a-j-i. ti,j-i;ai, j> 0 j ) = a, ja j x exp f X c ((,-i ) (, j )(i-1 j )' c (, j-i )+«5-c (i, j )л 2C (j ) ((ij )(i-rlj ) f 2 a,- I0 X a, a j-ir ,j (i - r j ) (36) 2 NONLINEAR MDOF SYSTEM RELIABILITY ASSESSMENT In this section the approximate analytical technique developed by some of the authors in [14] for nonlinear SDOF survival probability determination is generalized herein to account for MDOF systems by utilizing the dimension reduction/decoupling technique outlined in section 1. In this regard, the survival probability Pf is defined as the probability that the system response a amplitude a, stays below a prescribed barrier B over the time interval [0, T], given that a(t = 0) < B. Further, the first-passage PDF and the survival probability pi (T) are related according to the expression: ^ (T ) = . ' K ' dT (37) Next, adopting the discretization scheme employed in [9] yields intervals of the form: \_ti,j-i, j ], j = 1,2,..., m, tifl = 0, Teq,i (-1 ) ti,m = T, ti, j - ti,j-1 = 2 (38) where the response amplitude ai is assumed to be constant over [tiJ-b t j due to its slowly varying in time behavior. In Eq. (38) TeqA represents the LTV system equivalent natural period given by: (t ) = 2n ю ,, (t )■ (39) Note in passing that a smaller time interval can be chosen if higher accuracy is required. In this regard, the survival probability Pf is assumed to have a constant value over the same time interval as well. Obviously, the survival probability is given by if (T) = П [1 - FB ], (40) j=i where FB is defined as the probability that the response amplitude a, will exceed the prescribed barrier B over the time interval [tiJ_1, t^], given that no crossings have occurred prior to time tiJ_1. Next, invoking the Markovian property of the response amplitude a, one gets: FB = P™bla {tj )> в n a, (j )< B]_ Hjj Prob[ai [tiJ_l )< B] H (41) j-i where П denotes the intersection symbol. Utilizing Eq. (24) H1Bj_1 can be determined analytically in a straightforward manner; that is, HU-1 = \P (j-i , ti, j -1 ) dai, j -1 =1 - exp B2 A 2ci (j-i ) ,(42) whereas HB_ j is defined as a double integral of the form: rc B Hu-1, j = \dai, jx \p ( a л ' li,j-1 ; ai, j ' j ) dai,j-i. (43) Further, taking into account Eq. (36) and expanding the Bessel function I0(x) in the form (e.g., [27]): " (x / 2)) 'o (x) = Z k !Г(к +1) (44) analytical treatment of the involved integrals is possible yielding: Hj-1, j = A,0 + Z4> (45) where 4,o = exP ( ( -B l - exp 2c (,7 )( - Г* ) ' -B2 2c,- ()(l - Г2 ) x(l- r,2,), (46) L and A =_ (о.mo --йr n:,(2»y with L = 4"c. (t. . , )"+1 c (t. . )(1 - r2. )"+2 x 'и ' v j-1 ! j !\ uj ) b2 -, (47) Г Л 1 + и (j )( - r2j ) -Г. [1 + и,0] Л )(1 - r2))) Г, [1 + и] + B2 Л- ( X v v C. (ti,J )(1 r,2 ^ Г, [1 + )]-Г, 1 + ), B2 2c, (t,j )(1 - Г2. ) X B 2 ,(48) In Eq. (48) rt[y,z] represents the incomplete Gamma function defined as Г [y, z] = jf^e^'dt. Concisely, the developed technique comprises the following steps: i. Determination of the MDOF system non-stationary response covariance matrix (Eqs. (10) and (13)) via a statistical linearization treatment of the problem. ii. Determination of the equivalent linear time-varying elements ßeqj(t) and meqi(f) by solving the system of algebraic equations (Eqs. (17) and (18)). k=0 =1 X + X 0 IV. Determination of c,(t) via numerically Integrating the first-order ODE Eq. (25). Determination of the equivalent natural period Teq i(t) (Eq. (39)) and discretization of the time domain via Eq. (38). Further, v. Determination of the parameters H H i-1, i via Eqs. (42) and (43). j-i and vi. ",J B ( \ Determination of the survival probability P (T) via Eq. (40) and of the corresponding firstpassage PDF pi (T) via Eq. (37). 3 NUMERICAL APPLICATIONS In this section, a nonlinear three-degree-of-freedom system following the Bouc-Wen hysteretic model (e.g., [28] and [29]) subject to evolutionary stochastic excitation is considered to demonstrate the reliability of the technique. The survival probabilities and the first-passage PDFs obtained via the developed approximate technique are compared with survival probability and first-passage PDF estimates obtained via pertinent Monte Carlo simulations (10,000 realizations). The Monte Carlo simulations were conducted by utilizing a spectral representation methodology; additional details can be found in [30]. Further, a standard fourth-order Runge-Kutta numerical integration scheme is employed for solving the nonlinear system differential equation of motion (Eq. (1)), whereas the barrier level B is expressed as a fraction X of the maximum over time and over DOF value of the non-stationary response displacement standard deviation, i.e. B = Xmax(ai (t )) with __i and t (t) = yjci (t). Considering displacements defined in relative coordinates, the 3-DOF nonlinear system is governed by Eq. (1) where / = ( У2 Уз Z1 z2 z3 ) M = Mn M12 M21 M22. (49) (50) where and ml 0 0 " Mu = m2 m2 0 , (51) m3 m3 m3 Ml2 = M и = M 22 = 0з,3 . (52) K = Kii K i: K„ K,. where K„ = akx -ak2 0 0 ak2 акз 0 ak (54) K12 = (1 - a)kj -(1 - a)k, 0 0 0 (1 - a)k2 -(1 - a)k3 0 (1 - a)k3 and K21 K 22 03,3' (55) (56) In Eqs. (54) and (55) stands for the rigidity ratio which can be viewed as a form of post-yield to pre-yield stiffness ratio (a = 1 corresponds to the linear system). Further, the damping matrix of the structural system C is assumed to be proportional to the stiffness matrix; that is, C = C C '"'Il 12 C C 21 22 where C11 = cK11, C = C = 0 12 '-'21 03,3, and C = ^22 1 0 0 0 1 0 0 0 1 (57) (58) (59) (60) In Eq. (58) c is taken equal to 0.2 x 10 2. For the specific example y, = m,, and the loading vector becomes F (t) = (( (t) Л (t ) f (t ) 0 0 0). (61) Further, g (jy ) = = (000-gì (УРzi)-g2 (Уг,z2)-g3 (KZ). (62) In the Bouc-Wen model the additional state zt is associated with the displacement y, via the equation: z,- = gi (,zi ), (63) where g i (, Z ) = -y\ y\zi\zi I"-1 - ßy'i \zi \" + Ay i. (64) The parameters у, ß, A and n are capable of representing a wide range of hysteresis loops (e.g., [28] and [29]). In this example the values a = 0.15, ß=y=0.5, n = 1 and A = 1 are considered. The equivalent linear matrices take the form (e.g., [1] to [3]): = C C eqll eq12 C C eq 21 eq 22 where and C = C = C = 0 ^eq11 ^eq12 eq 22 "3,3' (65) (66) Ceq1 0 0 C = eq 21 0 Ceq 2 0 . (67) 0 0 c , eq3 Further, Keq = Keq11 Keq12 Keq 21 Keq22 where and Keqll Keq 12 K eq 21 03,3 Keq 22 Ki о о 0 к92 о 0 0 k eq3 (68) (69) (70) The elements c and k in Eqs. (67) and (70) are given by the expressions: E(y,z, ) л ш and ^ E (y,2 ) +ß E & - A, (71) (72) respectively. 3.1 A 3-DOF Hysteretic System under Evolutionary Stochastic Excitation of the Separable Form In this example, the excitation EPS Sa (m,г) takes the form Sä (^1 ) = |w()|2 SCP where SCP (a>) represents the widely used in engineering applications Clough-Penzien power spectrum (e.g., [31]) and w (t) denotes a time-modulating envelope function given by: w(t) = k ( - e-2 ), (74) where b = 0.1 and b2 = 0.3; and к is a normalization constant so that w (t)max = 1. The Clough-Penzien spectrum is given by: SCP (®) = S0 mg 4 + 4tfg )X У 2 X (2 -m2 ) + 4V x_(m/ mf)4_, (l -(m/ mf )2 ) + 4£f 2 (m/ mf )2 (75) where S0 is the amplitude of the excitation spectrum, modeled as a white noise process. The parameters values used are ^g=0.7, mg=2 rad s-1, l(t) and damping ßeq>l(t) elements corresponding to each DOF are plotted, respectively. Underlying the analytical approximate approach is the attempt to capture the time evolution as well as the essential characteristics of the frequency content of the nonlinear system response. Note that the ability of the technique to provide with time-varying natural frequencies rneqi(t) can be of particular importance if seen in conjunction with recent theoretical developments regarding the concept of the mean instantaneous frequency (MIF) (e.g., [33] to [35]). In this regard, rneq>(t) together with the MIF of the excitation can be potentially employed for evaluating the effects of temporal non-stationarity in the frequency content of the excitation on the system response as well as for tracking moving resonance phenomena (e.g., [23] and [36]). Further, in Figs. 11, 12 and 13 the survival probabilities Ff (T) for every DOF of the hysteretic MDOF system are plotted for various barrier levels, respectively; comparisons with MCS (10,000 realizations) demonstrate a satisfactory degree of accuracy. Fig. 11. Survival probability for various values of the parameter X for the first DOF; comparisons with MCS (10,000 realizations) Fig. 12. Survival probability for various values of the parameter X for the second DOF; comparisons with MCS (10,000 realizations) Fig. 13. Survival probability for various values of the parameter X for the third DOF; comparisons with MCS (10,000 realizations) 4 CONCLUDING REMARKS An approximate analytical technique for determining the time-varying survival probability and associated first-passage PDF of nonlinear/hysteretic MDOF systems subject to evolutionary stochastic excitation has been developed. Specifically, based on an efficient dimension reduction approach and relying on the concepts of stochastic averaging and statistical linearization, the original nonlinear w-degree-of-freedom system has been decoupled and cast into (w) effective single-degree-of-freedom (SDOF) linear time-variant (LTV) oscillators corresponding to each and every DOF. In this regard, time-varying effective stiffness a>lq t (t) and damping ßeqi (t) elements corresponding to each and every Dof have been defined and computed, while the non-stationary marginal, transition and joint response amplitude PDFs have been efficiently determined in closed-form expressions. Finally, the MDOF system survival probability and first-passage PDF have been determined approximately in a computationally efficient manner. Overall, the developed technique exhibits enhanced versatility since it can handle readily a wide range of nonlinear behaviors as well as various stochastic excitations with arbitrary non-separable EPS forms that exhibit strong variability in both the intensity and the frequency content. A 3-DOF system exhibiting hysteresis following the Bouc-Wen model has been included in the numerical examples section. Comparisons with pertinent Monte Carlo simulations have demonstrated the reliability of the technique. 5 ACKNOWLEDGEMENTS The first author acknowledges the financial support of the State Scholarships Foundation - IKY in Greece. 8 REFERENCES [1] Roberts J.B., Spanos P.D. (2003). Random Vibration and Statistical Linearization. Dover Publications, New York. [2] Soong, T.T., Grigoriu, M. (1993). Random Vibrations of Mechanical and Structural Systems. Prentice-Hall, New Jersey. [3] Li, J., Chen, J. (2009). Stochastic Dynamics of Structures, Dover Publications, New York, D0l:10.1002/9780470824269. [4] Schueller, G.I., Pradlwarter, H.J., Koutsourelakis, P.S. (2004). A critical appraisal of reliability estimation procedures for high dimensions. Probabilistic Engineering Mechanics, vol. 19, no. 4, p. 463-474, D0I:10.1016/j.probengmech.2004.05.004. [5] Corotis, R., Vanmarcke, E.H., Cornell, C.A. (1972). First passage of nonstationary random processes. Journal of the Engineering Mechanics Division, vol. 98, no. 2, p. 401-414. [6] Vanmarcke, E.H. (1975). On the distribution of the firstpassage time for normal stationary random processes. ASME Journal of Applied Mechanics, vol. 42, no. 1, p. 215-220, DOI:10.1115/1.3423521. [7] Barbato, M., Conte, J.P. (2001). Structural reliability applications of non-stationary spectral characteristics. ASCE Journal of Engineering Mechanics, vol. 137, no. 5, p. 371-382, D0I:10.1061/(ASCE)EM.1943-7889.0000238. [8] Kovaleva, A. (2009). An exact solution of the first-exit time problem for a class of structural systems. Probabilistic Engineering Mechanics, vol. 24, no. 3, p. 463-466, D0I:10.1016/j.probengmech.2009.01.002. [9] Solomos, G.P., Spanos, P.D. (1983). Structural reliability under evolutionary seismic excitation. International Journal of Soil Dynamics and Earthquake Engineering, vol. 2, no. 2, p. 110116, D0I:10.1016/0261-7277(83)90007-4. [10] Iourtchenko, D.V., Mo, E., Naess, A. (2006). Response probability density functions of strongly non-linear systems by the path integration method. International Journal of NonLinear Mechanics, vol. 41, no. 5, p. 693-705, D0I:10.1016/j. ijnonlinmec.2006.04.002. [11] Naess, A., Iourtchenko, D., Batsevych, O. (2011). Reliability of systems with randomly varying parameters by the path integration method. Probabilistic Engineering Mechanics, vol. 26, no. 1, p. 5-9, D0I:10.1016/j.probengmech.2010.05.005. [12] Kougioumtzoglou, I.A., Spanos, P.D. (2013). Response and first-passage statistics of nonlinear oscillators via a numerical path integral approach. ASCE Journal of Engineering Mechanics, vol. 139, no. 9, p. 1207-1217, D0I:10.1061/ (ASCE)EM.1943-7889.0000564. [13] Kougioumtzoglou, I.A., Spanos, P.D. (2014). Stochastic response analysis of the softening Duffing oscillator and ship capsizing probability determination via a numerical path integral approach. Probabilistic Engineering Mechanics, vol. 35, p. 67 - 74, D0I:10.1016/j.probengmech.2013.06.001. [14] Spanos, P.D, Kougioumtzoglou, I.A. (2014). Survival probability determination of nonlinear oscillators subject to evolutionary stochastic excitation. Journal of Applied Mechanics, vol. 81, no. 5, 051016, D0I:10.1115/1.4026182. [15] Spanos, P.D., Lutes, L.D. (1980). Probability of response to evolutionary process. Journal of Engineering Mechanics, vol. 106, no. 2, p. 213-224. [16] Kougioumtzoglou, I.A., Spanos, P.D. (2013). Nonlinear MDOF system stochastic response determination via a dimension reduction approach. Computers and Structures, vol. 126, p. 135-148, D0I:10.1016/j.compstruc.2012.11.020. [17] Priestley, M.B. (1965). Evolutionary spectra and non-stationary processes. Journal of the Royal Statistical Society, vol. 27, no. 2, p. 204-237. [18] Dahlhaus, R. (1997). Fitting time series models to non-stationary processes. The Annals of Statistics, vol. 25, no. 1, p. 1-37, D0I:10.1214/aos/1034276620. [19] Gajic, Z., Qureshi, M. (1995). Lyapunov Matrix Equation in System Stability and Control. Academic Press, New York. [20] Mitseas, I.P., Kougioumtzoglou, I.A., Beer, M. (2016). An approximate stochastic dynamics approach for nonlinear structural system performance-based multi-objective optimum design. Structural Safety, vol. 60, p. 67-76, D0l:10.1016/j. strusafe.2016.01.003. [21] Hammond, J.K. (1973). Evolutionary spectra in random vibration. Journal of the Royal Statistical Society, vol. 35, no. 2, p. 167-188. [22] Jangid, R.S., Datta, T.K. (1999). Evaluation of the methods for response analysis under non-stationary excitation. Shock and Vibration, vol. 6, no. 5-6, p. 285-297, DOI:10.1155/1999/312010. [23] Tubaldi, E., Kougioumtzoglou, I.A. (2015). Nonstationary stochastic response of structural systems equipped with nonlinear viscous dampers under seismic excitation. Earthquake Engineering & Structural Dynamics, vol. 44, no. 1, p. 121-138, D0l:10.1002/eqe.2462. [24] Giaralis, A., Spanos, P.D. (2010). Effective linear damping and stiffness coefficients of nonlinear systems for design spectrum based analysis. Soil Dynamics and Earthquake Engineering, vol. 30, no. 9, p. 798-810, D0l:10.1016/j.soildyn.2010.01.012. [25] Spanos, P.D., Solomos, G.P. (1983). Markov approximation to transient vibration. Journal of Engineering Mechanics, vol. 109, no. 4, p. 1134-1150, D0I:10.1061/(ASCE)0733-9399(1983)109:4(1134). [26] Spanos, P.D. (1978). Non-stationary random vibration of a linear structure. International Journal of Solids and Structures, vol 14, no. 10, p. 861-867, D0l:10.1016/0020-7683(78)90076-8. [27] Abramowitz, M., Stegun, I.A. (1970). Handbook of Mathematical Functions. Dover Publications, New York. [28] Wen, Y.K. (1980). Equivalent linearization for hysteretic systems under random excitation. Journal of Applied Mechanics, vol. 47, no. 1, p. 150-154, D0l:10.1115/1.3153594. [29] Ikhouane, F., Rodellar, J. (2007). Systems with Hysteresis: Analysis, Identification and Control using the Bouc-Wen Model. John Wiley and Sons, Chichester, D0l:10.1002/9780470513200. [30] Shinozuka, M., Deodatis, G. (1991). Simulation of stochastic processes by spectral representation. Applied Mechanics Reviews, vol. 44, no. 4, p. 191-204, D0l:10.1115/1.3119501. [31] Clough, R.W., Penzien, J. (1993). Dynamics of Structures, McGraw-Hill, New York. [32] Liu, S.C. (1970). Evolutionary power spectral density of strong-motion earthquakes. Bulletin of the Seismological Society of America, vol. 60, no. 3, p. 891-900. [33] Qian, S. (2002). Introduction to Time-Frequency and Wavelet Transforms. Prentice Hall, New Jersey. [34] Kijewski-Correa, T., Kareem, A. (2006). Efficacy of Hilbert and wavelet transforms for time-frequency analysis. Journal of Engineering Mechanics, vol. 132, no. 10, p. 1037-1049, D0l:10.1061/(ASCE)0733-9399(2006)132:10(1037). [35] Spanos, P.D., Giaralis, A., Politis, N.P., Roesset, J.M. (2007). Numerical treatment of seismic accelerograms and of inelasticseismic structural responses using harmonic wavelets. Computer-Aided Civil and Infrastructure Engineering, vol. 22, no. 4, p. 254-264, D0l:10.1111/j.1467-8667.2007.00483.x. [36] Beck, J.L., Papadimitriou, C. (1993). Moving resonance in nonlinear response to fully non-stationary stochastic ground motion. Probabilistic Engineering Mechanics, vol. 8, no. 3-4, p. 157-167, D0l:10.1016/0266-8920(93)90011-J. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 452-462 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3735 Original Scientific Paper Evaluation of Substructure Reduction Techniques with Fixed and Free Interfaces Fabian M. Gruber* - Daniel J. Rixen Technical University of Munich, Institute of Applied Mechanics, Germany Substructure reduction techniques are efficient methods to reduce the size of large models used to analyze the dynamical behavior of complex structures. The most popular approach is a fixed interface method, the Craig-Bampton method (1968), which is based on fixed interface vibration modes and interface constraint modes. In contrast, free interface methods employing free interface vibration modes together with attachment modes are also used, e.g. MacNeal's method (1971) and Rubin's method (1975). The methods mentioned so far assemble the substructures using interface displacements (primal assembly). The dual Craig-Bampton method (2004) uses the same ingredients as the two aforementioned free interface methods, but assembles the substructures using interface forces (dual assembly). This method enforces only weak interface compatibility between the substructures, thereby avoiding interface locking problems as sometimes experienced in the primal assembly approaches using free interface modes. The dual Craig-Bampton method leads to simpler reduced matrices compared to other free interface methods and the reduced matrices are sparse, similar to the classical Craig-Bampton matrices. In this contribution we evaluate the primal (classical) formulation of the Craig-Bampton method, the MacNeal method, the Rubin method and the dual formulation of the Craig-Bampton method. The presented theory and the comparison between the four substructuring methods will be illustrated on the Benfield truss, on a three-dimensional beam frame and on a two-dimensional solid plane stress problem. Keywords: dynamic substructuring, component mode synthesis, model order reduction, dual assembly, Craig-Bampton method, free interface method Highlights • Derivation of the Craig-Bampton, MacNeal, Rubin and dual Craig-Bampton method in a common framework. • Highlighting the differences and emphasizing the similarities between the four methods. • Assessing the accuracy of the methods. • Evaluation of the sparsity pattern of the reduced matrices. • Discussion of the negative eigenvalues inherent to the dual Craig-Bampton method. 0 INTRODUCTION Nowadays the increasing performance of modern computers makes it possible to solve very large linear systems of millions of degrees of freedom (DOF). Nevertheless, since dynamic analysis requires solving a lot of linear systems and since the refinement of finite element models is increasing faster than the computing capabilities, dynamic substructuring still remains an essential tool for analyzing dynamical systems in an efficient manner. Building reduced models of subparts of a structure enables sharing models between design groups. Moreover the reduction of the DOF of substructures is also important for building reduced order models for optimization and control. If a single component of a system is changed, only that component needs to be reanalyzed and the system can be analyzed at low additional cost. Thus dynamic substructuring offers a flexible and efficient approach to dynamic analysis. Dynamic substructuring techniques can be classified in two categories depending on the underlying modes which are used [1]. The term mode can refer to all kind of structural shape vectors. The first class consists of methods using fixed interface vibration modes and interface constraint modes to represent the substructure dynamics. The method commonly used is the Craig-Bampton method (CBM) [2] which assembles the substructures in a primal way using interface displacements in order to enforce interface compatibility. The second class consists of methods using free interface vibration modes and attachment modes. Common representatives of that class are MacNeal's method (MNM) [3] and Rubin's method (RM) [4] using a primal assembly process as well. Herting generalizes in [5] the concept of component mode synthesis to include any kind of interface boundary condition for the modes. In contrast to the aforementioned methods, the dual Craig-Bampton method (DCBM) [6] uses the same ingredients as MacNeal's and Rubin's method, but assembles the substructures in a dual way using interface forces. As a consequence, the DCBM enforces only weak interface compatibility between the substructures, thereby avoiding interface locking problems as sometimes experienced in the primal assembly approaches. Furthermore, the dual Craig-Bampton method leads to simpler reduced matrices 452 *Corr. Author's Address: Technical University of Munich, Institute of Applied Mechanics, Boltzmannstr. 15, 85748 Garching, Germany, fabian.gruber@tum.de compared to other free interface methods and the reduced matrices are sparse, similar to the classical Craig-Bampton matrices. In this contribution we evaluate the primal (classical) formulation of the Craig-Bampton method, the MacNeal method, the Rubin method and the dual formulation of the Craig-Bampton method. The presented theory and the comparison between the four substructuring methods will be illustrated on three different examples. In section 1, the differences between primal and dual assembly are stated. Following this, the formulation of the CBM, the MNM, the RM and the DCBM will be outlined in section 2 explaining general properties of the fixed interface method (subsection 2.1) and of the free interface methods (subsection 2.2). These properties will be illustrated subsequently in detail in section 3 using the Benfield Truss (subsection 3.1), a three-dimensional beam frame (subsection 3.2) and a two-dimensional solid plane stress problem (subsection 3.3). Finally a brief summary and conclusions are given in section 4. the equation of motion, Eq. (1), of one substructure partitioned in the same manner writes M (s) lylbb M (s) Mb s) M() r«bs)" + Г K (s) bb K «1 U(s)_ K (s) _ ib K (s) bb _ ' f(sy + _ f(s). 0 (3) Defining the local Boolean localization matrix A^ which is selecting the boundary DOF of substructure s gives the relation (4) which will be used later. 1.1 Primal Assembly The Eqs. (1) and (2) of all substructures N can be assembled in a primal way as: M u + K u = f , aa a a J a ' (5) where 1 PRIMAL AND DUAL ASSEMBLY OF SUBSTRUCTURES Consider a finite element model of a global domain. This domain is divided into N non-overlapping substructures such that every node belongs to exactly one substructure except for the nodes on the interface boundaries. The linear/linearized equation of motion of one substructure is written as M(s) U(s - K(s) u(s)= fw- g ( ) (1) where the superscript (s) is the label of the particular substructure ( = . Mw , Kw and u(s} are the mass matrix, the stiffness matrix and the displacement vector of the substructure, respectively. f ^ is the external force vector and g^ is the vector of reaction forces on the substructure due to its connection to adjacent substructures at its boundary DOF. The local displacements of each substructure can be divided in local internal DOF u and boundary DOF : (s) им = H 0 Ub 0 I (2) of the The boundary DOF ^ ^ ^ substructure are a subset of the assembled boundary DOF ub of the global domain. is a Boolean localization matrix (connectivity matrix) relating the assembled boundary DOF ub of the global domain to the substructure boundary DOF Ub}. Consequently, M = K = Mbb M(i) Mib MN ) Mib IK bb K(1) K ib ~ (N ) K ib M1) Mbi (1) M ? (1) K bi ? (1) with M( N ) M bi (N ) MM ? (N )' Kbi K 0 (N ) ub " f b ' «я f (1) «a = , fa = J i M bb = ^ Mib = M^D, Mbi = LfMb), ib b N K bb=YLb} kbSL?, Kib = Kbi = Lf Kb), ib b N K b = fb (s) f(S). (6) (7) (8) 0 s=1 s=1 s=1 The reaction forces g^ on the interfaces of the substructures cancel out during assembly. This assembly is called primal assembly since the compatibility between the substructures is enforced using the same boundary displacements for adjacent substructures. 1.2 Dual Assembly Another way to enforce the interface compatibility between the substructures is to consider the interface connecting forces as unknowns. These forces must be determined to satisfy the interface compatibility condition (displacement equality) and the local equation of motion of the substructures: У " B(su}= o, M(s)Us) + Ks)u(s) + B(s)T X = fs\ (9) (10) В(s) is a signed Boolean matrix (constraint matrix) acting on the substructure interface DOF. X is representing the interconnecting forces between substructures which is corresponding to the negative interface reaction force vector g(s) in Eq. (1) meaning ) = BT X, g g^=-Bp X, (11) and X is the vector of all Lagrange multipliers acting on the interfaces which are the additional unknowns. The displacement vector is partitioned according to Eq. (2). With the block-diagonal matrices M = M (i) 0 0 M (N ) (12) K = K (1) K (N) (13) and ' u«" - f u = , f= UN \ _ f(N В = Г B(1) B (N Г the set of Eqs. (9) and (10) can be written as: " M 0" U ' K BT u ' f " + = 0 0 _x_ B 0 x_ 0 (14) (15) (16) In this hybrid formulation the Lagrange multipliers X enforce the interface compatibility constraints and can be identified as interface forces [6]. The dual assembled system in Eq. (16) is equivalent to the primal assembled system in Eq. (5) since both systems express the same local equilibrium for each substructure and enforce the same interface compatibility. 2 COMPONENT REDUCTION METHODS 2.1 Craig-Bampton Method (CBM) Considering the partitioned equation of motion in Eq. (3), the internal DOF ufs of every substructure can be seen as being excited by its boundary DOF uf}, namely m^ uU(s + Kpu? = - K$ubs - m««« (17) This indicates that uf^ of each substructure can be approximated by a superposition of a static response and of eigenmodes associated to Mi? K(s). The static response is given by UfL=-K w- aib "b *ib "b , (18) where the columns of matrix are the static response modes also called constraint modes [2]. The fixed interface normal modes ф^ are obtained as eigensolutions of the generalized eigenproblem К^ф^ =ю()гМ^фР . The columns of the n( ) v n( ) m affi v Ф^ rr\ntain fbp first fixed ' X пф matrix Фф interface normal modes ф^ which can also be considered as the free vibration modes of the substructure s clamped on its boundary DOF ubs. The approximation of u\"s therefore writes u^* usa and the displacements approximated by Г J (s) 1bb 0(s)l Ж шМ _ ib 1ф _ n(s)_ (19) of the substructure are (20) with the vector of modal parameters П of dimension corresponding to the amplitudes of the fixed Ф interface normal modes Фф. The CBM reduction matrix TCB for reducing the primal assembled system can be defined as: 0 0 ш(1)г(1) H \n(N ) j(N ) ib b Фф} Фф) JN ) (21) matrix g(S) can be used instead of K(s) , which is defined by: ю (s)2 = K(s)+" 1 1 4-4=1 ю (s)2 (27) and the CBM reduced matrices Kred,CB = TCB KaTCB , (22) MredcCB = TBMJCB , (23) are found in [2]. 2.2 Free Interface Methods Considering the equation of motion (Eq. (10)) of substructure s, every substructure can be seen as being excited by the interface connection forces and the external forces. This indicates that the displacements of each substructure u(s) can be expressed in terms of local static solutions u^l and in terms of eigenmodes associated to the entire substructure matrices Kand M(s) : (s) (s) u ' = ul,a, ■ 1И(')-г(') Jj=1 with the static solution ) =-K (s )+BY я. и, Zf(') r R j=1 - (a j u j ■ (24) (25) The static response u)lt is obtained by solving Eq. (10) under the assumption of no inertia forces and no external forces acting on the substructure. K(s) is equal to the inverse of K when there are enough boundary conditions to prevent the substructure from floating when its interface with neighboring substructures is free [6]. If a substructure is floating, Kw is a generalized inverse of K(s) and R(s) is the matrix containing the rigid body modes as columns. The vector a(s) contains the amplitudes aj of the rigid body modes and the vector П contains the amplitudes of the local eigenmodes в being eigensolutions of the generalized eigenproblem (s) K WgW = > m(s>ejs>. An approximation is obtained by retaining only the first free interface r(s )a(s. normal modes в(s). Calling ©(s) the matrix containing only these eigenmodes, the approximation of the displacements of the substructure is given by: u{s) « -K(s)+B{s)TЯ + R{s)a{s) + &%{s). (26) Since a part of the subspace spanned by ©(s) is already included in K(s) the residual flexibility ,(s) Note that, by construction G^ = G^ , which is computed using the second equality in Eq. (27). For further properties of g(s) see [6]. As a result the approximation of one substructure writes: J' ) {r(-) в(-) -gs)b(-)t ] a -r R(' ' в'-' G(')A(' ' a Я (')' 1(-) J') (28) A^ is the matrix containing the residual flexibility attachment modes of substructure s, since the Boolean localization matrix as defined in Eq. (4) simply picks the columns of G^ associated to the boundary DOF [7]. The approximation in Eq. (28) can now be used to reduce the substructure DOF. Using the orthogonality properties of the modes in Eq. (28) the equation of motion of one substructure in Eq. (1) becomes M M iyl free n(s) + K (s) ^ "-free П g (s) S b g (s) S b = If(f(s) g (29) with the matrices fW = T(sf ' free A1 *W = Ж (s W = I 0 0 0 I 0 0 0 M « r ,bb _ 0 0 0 " 0 Q(s)2 0 0 0 G (s) r,bb ») A t s) л S)T . (30) (31) (32) (33) u a CB T b Gis the residual flexibility and My(lb is the interface inertia associated to the residual flexibility related to the boundary DOF, respectively, and being a diagonal matrix containing the remaining n^ (s) eigenvalues m (s) the reduction matrix tRs consistently to the mass and stiffness matrix resulting in a true Rayleigh-Ritz method as was observed in [8]. 2.2.2 MacNeal Method (MNM) 2.2.1 Rubin Method (RM) In order to assemble the substructure equation of motion in Eq. (29) in the global system a second transformation is applied by the RM [4]. The force DOF gbs'' are transformed back to the boundary displacements и^' using Eq. (28) [7]: u(s ) = A(s)u(s) = Rb as + П + g(s). (34) Rs) and ©is) are the subparts of R(s) and ©(s) related to the boundary DOF, respectively. From this equation the interface force DOF can be solved as: gw = k^ RV^-eiV'), (35) with = G(lb . The transformation matrix T2(s) from force dOf back to the boundary displacements Ubs leaving aand П unchanged writes then: The MNM [3] is nearly identical to the RM except for a small change. First we will derive the preliminary MNM reduced matrices K/^ш and following the derivation of the RM to show the similarities between these two methods. The reduced stiffness matrix of both the RM and the MNM are identical (given in Eq. (40)) >(•) = K (») KT -1 = к red ,MN "~red,R> (42) but the MNM reduced mass matrix M^ MN is obtained differently. The residual mass term M(s)b of the matrix M^ in Eq. (30) is neglected resulting in a modified matrix labeled M(s) = free,MN I 0 0 0 I 0 0 0 0 (43) instead of M^ for the MNM [7]. The preliminary MNM reduced mass matrix writes now: T(s) _ ±2 - K I 0 0 0 I 0 ■(s) R(s) r,bbb -K(s) ©(s) r,bb b K (s) r,bb (36) Application of this transformation to the matrices of Eqs. (30) and (31) gives the RM reduced matrices of one substructure K (s) _ T (sf K (s) T (s) red,R J2 "-free* 2 , M) = T)T M(s) T(s) lrl red,R 1 2 lrl free1 2 • (37) (38) The RM reduction matrix for one substructure writes therefore: T(s) = T (s)T ±R 2 ' and the RM reduced matrices As) - T(sf ,(s)T (s) Ks> = T[s> Ks>Ts red R 1V M M = TM red,R R *R > (39) (40) (41) are found [7]. These matrices can be directly assembled using primal assembly to get the RM reduced matrices KredR and MredR of the global system. This process was outlined in section 1.1 and applied in section 2.1 for the CBM. The RM applies M(s) = T red ,MN 2 ('f M (') T (') = M W lr± free,M^± 2 free,MN ' (44) This gives in fact inconsistent equations of motion since the mass and stiffness matrices are not reduced with the same basis. The assembly of the MNM reduced matrices KK(tsIMN and M(^dMN in the global system proceeds in the same manner as for the RM. Observing that the boundary DOF ub have no associated inertia in Eq. (44), those DOF can be condensed out of the equation of motion of the assembled problem and the final MNM reduced matrices KredMN and Mred гШ are obtained [3]. Thus the size of the assembled MNM system is reduced further by the number of DOF of ub. 2.2.3 Dual Craig-Bampton Method (DCBM) Replacing gbi according to Eq. (11) by the Lagrange multipliers X in the equation of motion of one substructure in Eq. (29), the reduced substructure matrices can be directly coupled using the dual assembly procedure [6] as outlined in section 1.2. Assembling all substructures N in a dual fashion by keeping the interface forces X as unknowns, the entire structure can consequently be approximated by a a (N ) ,(N ) (45) with the DCBM reduction matrix TD T = DCB rw @w r(n ) 0 ©(N ) 0 -G{1B(1)T I . (46) The approximation of the dynamic equations of the dual assembled system in Eq. (16) is M n(1) aN > + Kred,DCB > nN N X X = Tcb f, (47) with the DCBM reduced mass and stiffness matrix M = T1 red,DCB DCB ' M 0" " I 0" T= 0 0 DCB 0 Mr _ red,DCB JDCB K B B 0 with M_ Mr =Yb()G()M «g«B (s)T (48) (49) (50) and KredDCB are diagonal for the parts related to the different substructures. The coupling between the substructures is only achieved by the rows and columns related to X . The DCBM applies the reduction matrix TDCB consistently to the mass and stiffness matrix resulting in a true Rayleigh-Ritz method. The DCBM enforces only a weak compatibility between the substructures and does not enforce a strong displacement compatibility between the interfaces compared to many other common reduction methods [6]. Considering the system of Eqs. (9) and (10) multiplied by the reduction matrix TDcb , the last row of Eq. (47) results from M (1)u (1) + K (1U(1) + ß(1'Ä = f W 3 = f (1) M(N u(N ) + K(N )+ ß( )J 3 = f(N ) УN ßMU*=0 _ ć—ls=1 multiplied from left by the last row of TTCB which is [_ B (1)g« ... - GN ] I ]. (52) Replacing the strong interface compatibility condition of Eq. (9) by the weak form according to the multiplication of Eq. (51) by Eq. (52) can be interpreted as follows. Denote Afthe residual forces of substructure 5 resulting from the weak satisfaction of the local equilibrium of the substructure approximating the dynamics by a small number of free interface normal modes. Name A«(s) = G^ Afs the displacements these residual force Afwould create locally. Then the weak interface compatibility condition (Eqs. (51) and (52)) states that a compatibility error (i.e. an interface displacement jump) equal to the incompatibility of Au^s} is permitted [6]. Compared to MacNeal's and Rubin's method [3] and [4], the weak interface compatibility of the DCBM avoids locking problems occurring during the application of the aforementioned methods. Therefore, the approximation accuracy is improved [6]. But the fact that a weak interface compatibility is allowed in the DCBM implies that the infinite eigenvalues related to the Lagrange multipliers X in the non-reduced problem in Eq. (16) are now becoming finite and negative [9]. In practice those negative eigensolutions will appear only in the higher eigenvalue spectrum if the reduction space is rich enough [9]. Nevertheless, the reduction basis has to be selected with care avoiding potential non-physical effects of the possibly occurring negative eigenvalues. If Mr in Eq. (48) is neglected strong interface compatibility is enforced again and the DCBM reduced system with Mr = 0 is equivalent to the MNM [6]. Then static condensation can be applied again to remove X (as it was done for ub at the end of the derivation of the MNM in section 2.2.2) from the assembled system since no mass is associated. Thus the size of the assembled system is reduced again by the number of DOF of X . 3 EXAMPLES AND DISCUSSION 3.1 Benfield Truss The Benfield truss [10] of Fig. 1 is used to compare the results obtainable by the CBM, the MNM, the RM and the DCBM. The planar truss consists of two substructures having uniform bay section whereas all members have constant area and uniform stiffness and mass properties. The left component consists of five equal bays and has a total of 18 joints and the right component consists of four equal bays and has a total of 15 joints [10]. The lowest eigenfrequencies a of the entire structure shall be approximated by the different methods. The relative error interface methods for this example, the RM performs always better than the DCBM and the DCBM performs again always better as the MNM. The CBM and the DCBM result in similar frequency errors. The sparsity pattern of the reduced stiffness matrix Kred and reduced mass matrix MredMN of the CBM (Fig. 3), the MNM (Fig. 5), the RM (Fig. 6) and the DCBM (Fig. 4), respectively, illustrate the differences of the assembled reduced structures. Both the reduced stiffness matrix Kred and the reduced mass matrix Mred applying the CBM and the DCBM, respectively, have only diagonal entries for the subparts of each substructure. On the one hand the coupling between the substructures using the CBM is entirely achieved by the last rows and last columns in -"rel,j у mfuu,j of the /h eigenfrequency the mass matrix M , CB (Fig. 3b) and the remaining is used as a criterion to assess the accuracy of the different methods. Thereby full,j is the jt eigenfrequency of the full (non-reduced) system and ared j represents the jth eigenfrequency of the reduced system obtained by each method. Using 5 elastic (fixed or free interface normal modes) per substructure the relative errors srel depicted in the semi-log graph in Fig. 2 are resulting. Fig. 1. Benfield truss [10] Fig. 2. Relative error of eigenfrequency using 5 normal modes per substructure for the approximation of the lowest eigenfrequencies of the Benfield truss Since all methods give the correct rigid body modes only the relative errors of the elastic modes are plotted. All methods give a relative error less than 1 % for the first six eigenfrequencies. Comparing the free part is diagonal [2]. On the other hand the coupling applying the DCBM is entirely achieved by the last rows and last columns in the stiffness matrix KredDCB (Fig. 4a) and again the remaining part is diagonal [6]. The corresponding degrees of freedoms are either the interface displacements ub or the interface forces X but no direct coupling between the modal parameters of adjacent substructures occurs which ensures the sparse structure. Fig. 3. Sparsity pattern of the reduced matrices applying the CBM using 5 normal modes per substructure Fig. 4. Sparsity pattern of the reduced matrices applying the DCBM using 5 normal modes per substructure In contrast the sparsity pattern of the reduced stiffness matrix K. and the reduced mass matrix Mred obtained by the MNM and the RM, respectively, show full matrices. The MNM gives indeed an entirely diagonal reduced mass matrix MredMN (Fig. 5b) but causes always a full coupling between all DOF of all substructures via the reduced stiffness matrix Kredm (Fig. 5a). This makes the reusability of reduced models obtained by the MNM very inefficient and therefore nearly impossible from a practical point of view. The RM also causes a coupling between the substructures via interface displacements ub in the reduced stiffness matrix KredR (Fig. 6a) as well as in the reduced mass MredR (Fig. 6b). 5 10 а) К red,M N b) Mred,MN Fig. 5. Sparsity pattern of the reduced matrices applying the MNM using 5 normal modes per substructure 10 20 10 a) KredtR b) MredtR Fig. 6. Sparsity pattern of the reduced matrices applying the RM using 5 normal modes per substructure Moreover all DOF belonging to one reduced substructure are coupled with all other DOF of the same substructure which is why the reduced matrices of the RM are full for the substructure blocks and not diagonal. This result concerning the sparsity of the reduced matrices is outlined in Table 1 which shows the number n of non-zero elements in the reduced matrices Kred and Mred and the sum ntotal of both obtained by the different methods for this example. The reduced matrices of the CBM, the MNM and the DCBM contain a similar number of entries while the RM causes even for such a simple example a remarkable high number of entries. The number of entries of the MNM are comparable to the CBM and the DCBM but will increase dramatically if the number of substructures is increased since Kred will always be completely full. Table 1. Number n of non-zero elements in the reduced matrices obtained by the different methods for the Benfield truss (5 normal modes per substructure) CBM MNM RM DCBM n in Kred 40 216 314 196 n in Mred 118 16 354 50 ntotal 158 232 668 246 3.2 Beam Frame In [6] a three-dimensional frame made of steel beams (Young's modulus 210 GPa, Poisson's ratio 0.3, and density 7500 kgm3) schematically shown in Fig. 7 is considered. Each cell in the frame has a height of 0.35 m and a width and depth of 0.5 m. All outer beams have a hollow circular cross-section with the outside and inside diameters 0.02 m and 0.018 m. The diagonal members inside the cells have a solid circular cross section with diameter 0.008 m. The frame is divided into 5 substructures and again the objective is the approximation of the lowest eigenfrequencies w of the frame. Approximation of the eigenfrequencies of this system using the four presented methods with 4 normal modes per substructure (rigid body modes are not counted as free interface normal modes) is carried out. The results presented in [6] for the DCBM were obtained based erroneously on an incomplete set of free interface modes in the substructures using a simple Lanczos eigensolver: modes related to multiple frequencies were not determined by the simple Lanczos algorithm applied in that work. Fig. 7. Frame made of steel beams divided into 5 substructures [6] Considering the geometry of the beam frame structure of the arms in Fig. 7, the symmetry of the substructures is identifiable resulting in multiple eigenfrequencies of the substructures. The simple Lanczos eigensolver (no blocks, no restarts) used in [6] was not capable of determining multiple eigenvalues reliably [11]. Hence, when the simple Lanczos algorithm is used to determine the free modes of the substructure, modes that are multiple are missing and the representation of the lower spectrum is incomplete. Normal modes associated to the higher eigenfrequencies computed by the simple Lanczos eigensolver are used leading to poor approximation accuracy. Using the block Lanczos method and taking the lowest eigenfrequencies with associated subspaces of eigenvectors, the accuracy of the eigenfrequencies of the reduced global system obtained by the two reduction methods increases in the low frequency range [11]. In this work we are using an eigensolver capable of determining multiple eigenvalues for the comparison of the results obtained by the different methods. Again the eigenfrequencies ared of the reduced system applying the CBM, the MNM and the DCBM as well as the eigenfrequencies afull of the full (non-reduced) eigenfrequency s system are computed. The relative errors Srei,j - rredj - ^full, j \/ ®m,,j comparing the eigenfrequencies of the non-reduced model with the eigenfrequencies ared of the reduced system obtained by the reduction methods are plotted in Fig. 8. Nevertheless, as already observed in [6], this confirms that the accuracy of the DCBM in the low frequency range is two orders of magnitude better than the CBM. Comparing only the free interface methods, the RM performs slightly better than the DCBM. Both the RM and DCBM give a much better approximation accuracy than the MNM. 3.3 Two-dimensional solid In order to emphasize the weak coupling between the substructure interfaces applying the DCBM, as described in section 2.2.3, the problem of a two dimensional rectangle (density p, Young's modulus E, Poisson's ratio v=0.3) decomposed in 12 substructures as illustrated in Fig. 9 is considered. Each substructure is discretized by 16x9 bilinear four-noded elements (plane stress) and the structure is clamped on the left side in both directions. Again the objective is to approximate the lowest eigenfrequencies a of the entire structure with the four different methods. Using 8 normal modes per substructure (not including potential rigid body modes) the relative errors Srel,j -\Vred,j -mfull,j\/ ®full j depicted in Fig. 10 are resulting. Making use of such a large reduction basis the DCBM and the RM give excellent results in the low frequency range in comparison to the CBM and the MNM. Approximating the 20 lowest eigenfrequencies of the full structure, neither compatibility problems nor trouble with negative eigenfrequencies occur during the reduction process with the DCBM. Similar graphs depicting the relative errors erel as in Fig. 10 are obtained using different numbers of normal modes per substructure for the four methods. When considering a very small reduction basis for the substructures (using only a small number of normal modes), non-physical negative eigenvalues of the reduced problem show up in the low frequency Fig. 8. Relative error erelj of eigenfrequencyj using 4 normal modes per substructure for the approximation of the lowest eigenfrequencies of the entire beam frame Fig. 9. Two dimensional solid problem (dimensionless width lx = 16, height ly = 9 decomposed in 12 substructures; each substructure is discretized by 16x9 bilinear four-noded elements (plane stress) and the structure is clamped on the left side in both directions range. Using only 2 normal modes per substructure as reduction basis for this example negative values emerge among the 20 smallest absolute values of eigensolutions a2 applying the DCBM. In this case the first 14 of the lowest absolute values of a2 are positive and the associated eigenmodes are approximating the true eigenvalues and eigenmodes of the full system accurately. But, as shown in Table 2, the 15th eigenvalue is negative and the associated eigenmode depicted in Fig. 11 shows the non-physical behavior of this eigensolution. Fig. 10. Relative error srelj of eigenfrequencyj using 8 normal modes per substructure for the approximation of the lowest eigenfrequencies of the entire structure Table 2. Number j and associated eigenvalue a2 of the two dimensional solid reduced with the DCBM „2 P ' E 0.000874 0.009729 14 0.245011 15 -0.250785 16 0.261676 The eigenmodes corresponding to negative eigensolutions with higher absolute values show similar behavior. All these negative eigenvalues are related to the weak compatibility on the interfaces and not meaningful from a physical point of view. Consequently detecting and filtering out those negative eigenvalues is an additional step in the DCBM compared to the other methods based on primal assembly. This extra step is cheap and therefore does not increase the effort of the reduction process using the DCBM compared to the other four methods keeping the DCBM very efficient. Fig. 11. Eigenmode number 15 approximated by the DCBM using 2 free interface normal modes per substructure e ( m2l5 = -0.250785 — ) P 4 CONCLUSIONS In this paper the general concepts of the Craig-Bampton method (CBM), the MacNeal method (MNM), the Rubin method (RM) and the dual Craig-Bampton method (DCBM) were briefly presented, compared and discussed using three examples. The DCBM is outperforming the CBM using the same number of normal modes per substructure as reduction basis with comparable computational effort and having similar sparsity pattern of the reduced matrices. Comparing the free interface methods, the RM performs slightly better than the DCBM but results in full matrices. Both the RM and DCBM give a much better approximation accuracy than the MNM while the MNM generated always full coupled reduced matrices. Properties of the DCBM were outlined and an additional necessary step, namely filtering out the negative eigensolutions, during the reduction process was illustrated. Non-physical negative eigenvalues of the reduced dual assembled problem are intrinsic in the reduction process using the DCBM caused by the weak compatibility on the interfaces between the substructures. Filtering out these negative eigenvalues is the decisive factor for the excellent approximation quality of the DCBM. The numerical effort adding this additional step is negligible keeping the efficiency of this method. 5 REFERENCES [1] Craig, R.R. (2000). Coupling of substructures for dynamic analyses: An overview. Proceedings of AIAA/ASME/ASCE/ AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, Atlanta, p. 1573-1584. 2 [2] Craig, R.R., Bampton, M.C. (1968). Coupling of substructures for dynamic analyses. AIAA Journal, vol. 6, no. 7, p. 13131319, D0l:10.2514/3.4741. [3] MacNeal, R.H. (1971). A hybrid method of component mode synthesis. Computers & Structures, vol. 1, no. 4, p. 581-601, D0I:10.1016/0045-7949(71)90031-9. [4] Rubin, S. (1975). Improved component-mode representation for structural dynamic analysis. AIAA Journal, vol. 13, no. 8, p. 995-1006, D0I:10.2514/3.60497. [5] Herting, D.N. (1985). A general purpose, multi-stage, component modal synthesis method. Finite Elements In Analysis and Design, vol. 1, no. 2, p. 153-164, D0I:10.1016/0168-874X(85)90025-3. [6] Rixen, D.J. (2004). A dual Craig-Bampton method for dynamic substructuring. Journal of Computational and Applied Mathematics, vol. 168, no. 1-2, p. 383-391, DOI:10.1016/j. cam.2003.12.014. [7] Voormeeren, S.N., Van Der Valk, P.L.C., Rixen, D.J. (2011). Generalized methodology for assembly and reduction of component models for dynamic substructuring. AIAA Journal, vol. 49, no. 5, p. 1010-1020, DOI:10.2514/1.J050724. [8] Craig, R.R., Chang, C.J. (1977). On the use of attachment modes in substructure coupling for dynamic analysis. 18th Structural Dynamics and Materials Conference, San Diego, p. 89-99, D0I:10.2514/6.1977-405. [9] Rixen, D.J. (2011). Interface Reduction in the Dual Craig-Bampton method based on dual interface modes. Linking Models and Experiments, vol. 2, Conference Proceedings of the Society for Experimental Mechanics Series, p. 311-328, D0I:10.1007/978-1-4419-9305-2_22. [10] Benfield, W.A., Hruda, R.F. (1971). Vibration analysis of structures by component mode substitution. AIAA Journal, vol. 9, no. 7, p. 1255-1261, D0I:10.2514/3.49936. [11] Gruber, F.M., Rutzmoser, J.B., Rixen, D.J. (2015). Comparison between primal and dual Craig-Bampton substructure reduction techniques. Proceedings of the 11th International Conference on Engineering Vibration, Ljubljana, p. 12451254. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 463-470 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3770 Original Scientific Paper Modelling the Belt - Envelope Interactions during the Postal Mail Conveying by a Sorting Machine Lionel Manin1* - Stefan Braun1 - Damien Hugues2 1 Univ. Lyon, INSA-Lyon, France 2 Solystic SAS, France The purpose of this study is to predict the trajectory of an envelope when conveyed between two flat belts in a mail sorting machine. Each of the two sides of the envelope is in contact with one belt. The friction coefficient can vary between the two sides, e.g. for the cases of a postcard or a plastic windowed envelope. The two belts are assumed to travel horizontally with the same velocities, but in real world small velocity differences can occur. Also, some deviation from the horizontality of the belts can be observed. These phenomena can lead to unexpected situations where the relative position of the envelope to the belt changes during transportation. The work focuses on the belt-envelope mechanical interactions in order to predict the envelope position during its conveying. A 2D model is developed which considers the dynamic equilibrium of the envelope and simulates the transient behavior at different locations in the sorting machine. In the dynamic friction model, the contact surfaces between the envelope sides and the belts have to be determined at each time step. The belt/envelope friction forces are applied at the centers of these contact surfaces One side of the envelope is considered stuck (in adhesion) to the belt while the other is considered slipping. The position of the envelope center of gravity and the angle with the horizontal axis are calculated. For evaluation of the proposed dynamic model, several cases of operation are simulated. Keywords: mail sorting machine, belts, envelope, conveying, simulation Highlights • Modelling the belt-envelope interactions in mail sorting machine. • Calculation of belt-envelope contact surfaces based on polygon clipping theory. • 2D modelling, simulation of the envelope trajectory. 0 INTRODUCTION The research presented in this paper is based on the mechanics of mail conveying in sorting machines as used in postal sorting centres. The general concept of these machines is illustrated in Fig. 1. A lot of different technologies are used in such systems: automatic address scanning and recognition, coding printer, etc. For robust machine operation, the correct conveying of the letters has to be ensured, which is equivalent to sustaining a stable position of the envelopes with respect to the belts that transport them. To give an idea, the machine can handle up to 53,000 mail pieces/ hour depending on the length of the pieces processed. The dynamics of belt transmission has been given a lot of interest and the related litterature is dense. Some of the earliest published works on the topic [1] to [3] were focused on the belt mechanics quasi-static behavior, and the belt-pulley contact friction. A reference paper about the first research on belt vibrations can be found in [4]. Progressively, an following the industry needs, the field of investigation moved to dynamics, durability, vibrations, nonlinear dynamics and stability of belt transmissions as exposed in [5] to [7]. The nonlinear phenomena often encountered in flat or poly-v belt transmissions is caused by parametric excitations due to pulley speed and torque fluctuations that introduce a pulsating tension in the belt [7] to [8]. In most of the studies, an analytical formulation is realized. Several other studies were developed using a finite element modelling [9], or also multi-body dynamics theory [10]. The belt dynamics has also been investigated through the topic 'dynamics of axially moving continua' which can be applied to string-like or beam-like or plate-like moving systems. In these works, most of the time only a part of the continua is considered with some boundary conditions [11]. In [12] a thorough recent review on the dynamics of axially moving continua was done, it starts from the earliest works and it details the main contributions of works published on that topic until 2014. It is suprising to see that among this very extensive litterature, no work considers the transportation of mail and the interaction with moving continua. In the scope of this research, a simulation tool is developed, it predicts the trajectory of a letter being conveyed between two flat belts and crossing several direction switches (Fig. 2) where the envelope/belt contact surface is reduced (Figs. 4a and 5). *Corr. Author's Address: Univ Lyon, INSA LyonCNRS UMR5259, LaMCoS, F-69621, France, lionel.manin@insa-lyon.fr 463 aj b) Fig. 2. aj Envelope conveying, bj envelope passing through a direction switch 1 MODEL DESCRIPTION The proposed dynamic model is two-dimensional. The envelope, the front and back side belts are considered as vertical planes as illustrated in Fig. 3. The belts have constant and independent velocities and can be inclined with respect to the horizontal axis. The aim is to determine the envelope position as a function of time. The dynamic equilibrium of the envelope is written at the mass centre G. The forces considered are the weight of the envelope and the friction forces of the belts onto the envelope sides. The contact surfaces have to be determined at each time step of the simulation process. The proposed method is presented in the next paragraphs. 1.1 Parameter Definition The envelope is defined by: its height hp and width lp, its mass mp and inertia Ip, the posi-tion of the plastic window if any and its size (hf * lf), the mass centre coordinates (xG, yG), the left bottom corner coordinates (xOp, y0p), the initial angle with respect to the horizontal axis вр. The belts are defined each by the coordinates of the two lower corner points and the width (Fig. 3). From these geometry input parameters, the angles between the belts and the horizontal axis are deduced. When the envelope is travelling through a direction switch (Fig. 2) the contact surface with the front belt(s) are more complicated to determine than a) Fig. 3. Geometry parameters for the belts and the conveyed mail envelope, general case including a sorting switch b) Fig. 4. a) Contact cases between belts and envelope when travelling through a direction switch, b) conveying with several direction switches the contact surface with the backside belt. This is due to the front belt being divided into two neighboring belts, introducing a dis-continuity. As the envelope passes this discontinuity, five contact cases can occur as shown in Fig. 4a. These cases are implicitly taken into account with our parameter definition. A standard mail sorting machine contains several successive direction switches which can be rep-resented by the scheme in Fig. 4b. In that case, multiple front belts are considered as on the machine. 2.2 Forces Acting on the Envelope Theoretically the front and rear belts have the same velocities. However, a small difference can exist which is estimated at around 5 %. We have made several assumptions when setting up the model: the envelope is adhering to the slowest belt (rear or front) and is sliding with respect to the fastest belt(s). To give an example: if the rear belt is slower than the two front belts (Fig. 5) then the envelope is Fig. 5. Polygonal contact surfaces and friction forces between belts and envelope subjected only to the friction forces applied by the front belts. As a consequence of these assumptions, the envelope minimum traveling velocity is equal to that of the slowest belt. The minimum contact pressure that prevents the envelope from falling freely can be estimated from Eq. (1), where Sp total is the total contact surface with the envelope paper, Sf the contact surface with the plastic window, and the friction coefficients fp and f (belt/envelope, plastic window/ envelope). If Eq. (1) is satisfied then the envelope does not fall. P\SpJo,al ■ fp + Sf ■ ff ) mp ■ g. (1) The friction forces are calculated for each contact surface S using the belt/envelope contact pressure p (supposed constant and uniform), the angles of the belts ecj, the friction coefficients f and the velocities of the belts (Eqs. (2) and (3)). In Eq. (4) the terms inside the parenthesis are considered only if the Eq. (1) is not satisfied. Fxi = sign (Vfmntj - V_ )• f • p • 5,. • cosecj, (2) Fyi = sign( - V_)• ] • p • 5,. • sinecj. (3) Z^ (+)( ■ g+ff ■ p ■ Sf+fp ■ p ■ SP_.««i) V i M« (G) = X(( • (( - ) + Fyi ■ ( - xG)) + + E((j '(o - ySj ) + Fy,j Sj - Xg )) (5) (4) • the sum of the moment of the external forces calculated at point G equalling the envelope inertia times the angular acceleration. The vector of the external forces applied to the envelope is given in Eq. (4), where the subscript p denotes the friction forces on the envelope paper, the subscript f denotes the friction forces on the envelope plastic window. The moment of the friction forces expressed at point G is given in Eq. (5). The coordinates of each contact surface centre Gt are (xSi, ySi). Finally, the equations of motion of the envelope are given in Eq. (6). An Euler iterative numerical scheme is used to solve Eq. (6), as detailed in Eqs. (7) to (9). One should note that at each time step the contact surfaces and their centre coordinates change and have to be re-calculated. Some initial conditions are considered depending on the running conditions simulated (Eq. (10)). \xn+1 =xn + dt • Fex,„ I X„+1 = Xn + dt •K dt dya_ dt dt Уо \xn+1 = xn + dt • FeK I Xn+1 = Xn + dt rXn (7) (8) (9) n p dt2 d 2в„ dt F., Mext (G )_ where I p = m (( + hp). (6) As already mentioned, a key point is the determination of the different contact surfaces St and their mass centres Gt where the friction forces are applied. This is developed in the next paragraph. The envelope is subjected to its proper weight and to the friction forces applied by the front and rear belts. The dynamic equilibrium of the envelope in the vertical plane can be written with: • the sum of external forces equalling the mass of the belt times its acceleration, XG_ 0 \V rear yG_ 0 , X0 = 0 ®p_ 0 _ 0 (10) 0 0 0 , Fext_ 0 = 0 0 0 2.4 Calculation of the Contact Surfaces The contact surfaces S are represented by polygons Pi (Fig. 5). They are the intersections between the belts and the two sides of the envelope with one side including a plastic window. All these parts are polygons that are defined by the coordinates of their four corners. The proposed method relies on getting the intersection polygons and the coordinates of their x G в X0 ~ ext ,y X0 ~ Fig. 6. Contact surfaces between the front belts and the envelope paper (green/daker) and plastic window (yellow/lighter) 2.5 3 X [m] Fig. 7. Representation of the mail sorting machine with the parameters of Table 1 points in order to determine the contact surfaces and their mass centres (Fig. 6). In literature, this problem is known as polygon clipping, often found in computer graphics [13]. In order to determine the intersection polygons, we use a polygon clipping technique from computer graphics. Many algorithms have been proposed in literature to solve this problem. For our implementation, we use the Vatti clipping algorithm [14]. This algorithm enables us to determine the intersection of arbitrarily shaped polygons in a two dimensional space. Once the intersection polygons are determined, the surface areas and mass centres are calculated using the 'polyarea' and 'centroid' functions of Matlab. one without and one with a plastic window. The envelope is adhering to the rear belt and the rear belt is slower than the front belt, no direction switch is Table 1. Geometry parameters of the belts Front belts Yi X4 Y4 lb lb [m] [m] [m] [m] [m] [m/s] 1 0.500 0.044 1.500 0.044 0.030 4.00 2 1.540 0.044 3.000 0.044 0.030 4.00 3 3.040 0.044 4.500 0.044 0.030 4.00 Rear belts Xi Yi X4 Y4 lb lb [m] [m] [m] [m] [m] [m/s] 1 0.500 0.044 5.000 0.044 0.030 3.99 2.5 Input Data The data defining the system are gathered in an Excel file which is read by the simulation program under the Matlab environment. An example is given in Table 1. From these input data, the system can be represented graphically as shown in Fig. 7. Table 2. Envelope and simulation parameters 3 RESULTS 3.1 Test Cases In order to validate the model developed, some simulation test cases have been considered. They are based on real observations made on a mail sorting machine, and also on some intuition on how the envelope position will evolve under different running conditions. The test cases are presented in Fig. 8. All tests were done for two different types of envelope: Name Value Description mp 0.02 envelope mail mass [kg] lP 0.219 mail width [m] hp 0.108 mail heigth [m] X лор 1 x coordinate left bottom corner [m] Уор 0.005 y coordinate left bottom corner [m] во 0 initial mail angle [deg] xf 0.1 window x position (m] yf 0.02 window y position [m] If 0.1 window width [m] hf 0.035 window heigth [m] g 9.81 gravity [m/s2] fp 0.4 paper/belt friction coefficient ff 0.41 plastic/belt friction coefficient dt 0.001 simulation time step [s] P 60 contact pressure [Pa] Fig. 8. Simulation test cases, (1) inertia centre centred with the belts, (2) inertia centre above belt centre, (3) inertia centre below belt centre, adhesion to the rear belt Vfront > Vrear considered. When the inertia centre of the envelope is centred with the belts, no rotation is observed for the envelope without window (Fig. 8 (1-b)), but a slight rotation is observed for the envelope with a window (Fig. 8 (1-d)). This is due to the difference of friction coefficients between the paper (fp) and the plastic (f) with the belt. When ff> fp , the consequence is the generation of an anticlockwise moment that initiates the envelope rotation. The influence of the envelope inertia centre position with respect to the belt is highlighted with the simulation cases 2 and 3 in Fig. 8. A rotation is generated by the unbalance of the friction force moments. This phenomena is amplified when the envelope has a plastic window. The final positions obtained for the several test cases correspond to those that were expected and confirm the viability of the model. The vertical position of the inertia centre yG and the rotation angle 6p are plotted versus the conveying distance in Fig. 9. When the envelope is passing through a direction switch, there are three phases where the contact surfaces are reduced (Fig. 4a) and there is a risk of falling for the envelope if the vertical friction forces are smaller than the envelope weight, as illustrated in Fig. 9. Several running conditions can be simulated, we only present here some results to show the capabilities of the model. The friction coefficient values were given by the sorting machine manufacturer that had made some measurements. The important point for the simulations is that the paper/ belt friction coefficient should differ from the plastic/ belt friction coefficient, so that the effect of the plastic window on the envelop trajectory can be identified. 3.2 Simulation of a Conveying Several generic cases of mail conveying are considered (Fig. 10). The graphical interface enables the visualization of the successive positions taken by the envelope. Some alerts have been implemented in the code so that when the envelope reaches a position which is no more realistic, the simulation is stopped. For example, if the lower point of the envelope has a negative vertical coordinate, or the rotation angle is larger than 30° then the simulation is stopped. 4 CONCLUSIONS The work described in this paper deals with the simulation of the mail conveying in a postal sorting machine. It is based on an industrial demand and therefore corresponds to real issues encountered in such machines. The intention was to develop a simple simulation model in order to predict the envelope trajectory in the mail sorting machine under different Fig. 9. Envelope vertical position and angle qp while passing through a direction switch Fig. 10. Snap shots of the envelope while being conveyed operating conditions. The model developed considers simple or windowed envelopes. Some assumptions were made in order to simplify the analysis and the dynamic equilibrium. The forces considered are the envelope's weight and the friction forces applied by the front and rear belts between which the envelope is maintained. The data definition permits the simulation of machines with several direction switches. These switches happen to be the most critical situations in the conveying process. The proposed model is very flexible, as the dimensions of the envelope can be changed as well as all the other parameters describing the geometry of the conveying path. The core of the model relies on the calculation of the contact surfaces between the envelope sides and the belts for each simulation time step. This is based on the theory of polygon clipping that permits the computation of intersection surfaces. Evaluation of several test cases showed that the results obtained by simulation are consistent with real world observations and permit to explain several phenomena encountered. 5 ACKNOWLEDGEMENTS The authors are very grateful with the Solystic Company that has partly financed the work presented. 6 REFERENCES [1] Firbank, T.C. (1970). Mechanics of the belt drive. International Journal of Mechanical Sciences, vol. 12, no. 12, p. 10531063, DOI:10.1016/0020-7403(70)90032-9. [2] Childs, T.H.C. (1980). The contact and friction between flat belts and pulleys. International Journal of Mechanical Sciences, vol. 22, no. 2, p. 117-126, DOI:10.1016/0020-7403(80)90048-X. [3] Gerbert, G. (1999). Traction Belt Mechanics, MVD, Chalmers, Göteborg. [4] Abrate, S. (1992). Vibrations of belts and belt drives. Mechanism and Machine Theory, vol. 27, no. 6, p. 645-659, DOI:10.1016/0094-114X(92)90064-O. [5] Parker, R.G., Lin, Y. (2004). Parametric instability of axially moving media subjected to multifrequency tension and speed fluctuation. Journal of Applied Mechanics, vol. 68, no. 1, p. 4957, DOI:10.1115/1.1343914. [6] Pellicano, F., Catellani, G., Fregolent, A. (2004). Parametric instability of belts: theory and experiments. Computers & Structures, vol. 82, no. 1, p. 81-91, DOI:10.1016/j. compstruc.2003.07.004. [7] Michon, G., Manin, L., Remond, D., Dufour, R., Parker, R. (2008). Parametric instability of an axially moving belt subjected to multifrequency excitations: experiments and analytical validation. Journal of Applied Mechanics, vol. 75, no. 4, DOI:10.1115/1.2910891. [8] Michon, G., Manin, L., Parker, R.G., Dufour, R. (2008). Duffing oscillator with parametric excitation: analytical and experimental investigation on a belt-pulley system. Journal of Computational and Nonlinear Dynamics, vol. 3, no. 3, p. 1-6, D0l:10.1115/1.2908160. [9] Leamy, M.J., Wasfy, T.M. (2002). Transient and steady-state dynamic finite element modeling of belt-drives. Journal of Dynamic Systems, Measurement, and Control, vol. 124, no. 4, p. 575-581, D0I:10.1115/1.1513793. [10] Čepon, G., Manin, L., Boltežar, M. (2009). Introduction of damping into the flexible multibody belt-drive model: A numerical and experimental investigation. Journal of Sound and Vibration, vol. 324, no. 1-2, p. 283-296, D0I:10.1016/j. jsv.2009.02.001. [11] Čepon, G., Boltežar, M. (2007). Computing the dynamic response of an axially moving continuum. Journal of Sound and Vibration, vol. 300, no. 1-2, p. 316-329, D0I:10.1016/j. jsv.2006.08.014. [12] Marynowski, K., Kapitaniak, T. (2014). Dynamics of axially moving continua. International Journal of Mechanical Sciences, vol. 81, p. 26-41, DOI:10.1016/j. ijmecsci.2014.01.017. [13] Agoston, M.K. (2005). Computer Graphics and Geometric Modeling - Implementations and Algorithms, Springer, London, D0I:10.1007/b138805. [14] Vatti, B.R., (1992) A generic solution to polygon clipping. Communications of the ACM, vol. 35, no. 7, p. 56-63, D0I:10.1145/129902.129906. Strojniški vestnik - Journal of Mechanical Engineering 62(2016)7-8, 471-482 © 2016 Journal of Mechanical Engineering. All rights reserved. D0l:10.5545/sv-jme.2016.3774 Original Scientific Paper Assessment of the Fatigue Parameters from Random Vibration Testing: Application to a Rivet Joint Martin Česnik - Janko Slavič - Miha Boltežar* University of Ljubljana, Faculty of Mechanical Engineering, Slovenia In order to estimate a structure's fatigue life when excited with an acceleration profile the fatigue parameters must be known. However, material's fatigue exponent and fatigue strength are not always readily available, especially for complex structures that include riveted or welded joints for which additional fatigue tests are needed. This study introduces a new fatigue-parameter assessment method based on random vibration loading and its application to a blind-hole rivet joint that diminishes the need for additional fatigue tests. The presented procedure requires a simple experimental setup; however, a more extensive analysis of the experimental results is necessary. The method of fatigue parameter assessment is presented and applied on real, experimentally obtained data from vibration tests of rivet-joint specimens, excited with a random base-vibration load in the frequency range of a single natural frequency. Special attention was given to the modelling of the rivet joint and the uncertainties arising from the riveting process were considered. With the presented procedure it is possible to obtain the fatigue parameters solely from the results of random-vibration testing with different acceleration profiles and therefore diminishing the need for additional classic fatigue tests. The obtained fatigue parameters indirectly include the stress concentration factor and the damping-loss-factor increase during the damage accumulation. Additionally, by applying the random-vibration load the influences of the natural-frequency shift and the small nonlinearities of the structure are reduced, which can present a major issue in classic harmonic-vibration fatigue testing. Keywords: vibration fatigue testing, fatigue parameters, random base-excitation, frequency-domain counting method Highlights • A method for fatigue parameter assessment from vibration tests with electro-dynamic shaker is presented. • The tested specimens were excited with acceleration-controlled random signal. • Frequency based Tovo-Benasciutti method was applied to the parametrized numerical model to determine expected fatigue life of a particular specimen. • The method was applied to actual testing of blind-hole rivet joint with nonlinear response. 0 INTRODUCTION In a product's design and development stages standard accelerated vibration tests [1] and [2] are performed in order to prove that no fatigue failure due to the vibration load will occur within the expected service life. Clearly, for the case of vibration loading the structure's own dynamic response directly imposes a certain stress load to the structure that is greatly amplified in the vicinity of natural frequencies [3] and [4]. The introduced stress state together with the material's fatigue parameters directly leads to the damage accumulation rate [5], which can be obtained in time- or in frequency domain [6]. Since the structure's dynamic response to the random vibration excitation is commonly described in terms of a response power spectrum the damage-counting methods in the frequency domain [7] and [8] are the most convenient for a vibration fatigue analysis. In this way the analyst is able to predict whether the structure will withstand the given vibration test beforehand at the design stage. During a service life the fatigue failure often occurs at the construction elements used for joining different parts (e. g. welds, solder joints, rivets), where high stress concentrations appear. The fatigue assessment for the quasi-static harmonic or random loading of welded [9] and riveted joints is well researched field for safety assurance of airplanes [10] and [11] and offshore structures [12]. However, when the stress load leading to the fatigue failure depends not only on the external loading but also on the structure's own dynamic response the phenomenon exceeds classical fatigue scope and should be regarded as vibration fatigue problem [13]. Accordingly, much initial research in the field of vibration fatigue has been applied to the case of microelectronic solder joints (e. g. ball grid array (BGA), plastic ball grid array (PBGA)) in order to increase the reliability of the electronic components. Liu et al. [14] observed the propagation of a crack in a solder joint for the case of a sinusoidal vibration load. Additionally, Chen et al. [15] combined harmonic vibration testing with a numerically obtained stress response to assess the S-N curve of solder material. Using the harmonic-excitation approach for assessing material's fatigue parameters Yu et al. [16] estimated the solder joint fatigue life under random vibration loading. By further work of George et al. [17] and Česnik et al. [4] harmonic vibration testing has been shown to present an adequate alternative to classic fatigue testing for fatigue parameter assessment. On the other hand, several studies dealing with vibration fatigue for the case of random vibration have been presented. Pagnacco et al. [18] optimised the thickness distribution of a plate structure to prolong the structure's fatigue life when exposed to random vibration. Furthermore, Paulus and Dasgupta [19] proposed a semi-empirical model for a fatigue-life estimation in the case of random vibration. Han et al. [20] performed a vibration fatigue analysis of a spot-welded structure with an iterative approach. Recently, significant advances with the applicative note to the random vibration fatigue analysis of real components from automotive industry were made by Braccesi et al. [21]. In studies [4] and [15] to [20] the fatigue parameters were obtained either with standard fatigue tests or with harmonic vibration excitation of the specimen. In any case an additional test with a different experimental setup was needed to obtain the fatigue parameters. In this study a new methodology is presented that enables the assessment of fatigue parameters (e. g. fatigue exponent b and fatigue strength C) only from a series of random-signal base-excitation vibration tests in the frequency range of specimen's dynamic response. Using the developed methodology the need for additional fatigue tests is diminished. The presented method introduces several enhancements to the application of the vibration fatigue analysis to real structures. Firstly, by identifying the fatigue parameters it provides an analyst with a physically consistent evaluation of the random-vibration tests that resulted in the fatigue failure. Secondly, it allows the fatigue testing with electro-dynamic shaker, even when the tested structure shows large changes of modal parameters and high nonlinearities [22], which in the case of harmonic base-excitation [4] and [17] present a main source of unstable response. Thirdly, the presented methodology improves Allegri and Zhang's approach [23] with inverse power law by applying physically more consistent Tovo-Benasciutti counting method to the real experiment of the dynamic structure on the electro-dynamic shaker. To show the applicative contribution of this methodology the vibration fatigue of blind-hole rivet is studied, which fastens a metal plate to the rigid base structure. A blind-rivet joint presents an economically favorable alternative to the bolted joint when a shell structure is required to be fastened to the solid base structure, for example a thermal shield in the vehicle's exhaust system. Here, the blind-rivet joint is the part of the structure where the lowest local stiffness occurs; hence a large strain and stress occur in the rivet itself when the structure is exposed to a vibration load, so leading to the vibration fatigue failure of the rivet. Opposed to the several research that study the influence of rivets to the fatigue life riveted plate [24] to [28], this research deals with the fatigue damage accumulation on the rivet itself as a result of the dynamic response of the riveted plate structure. The nonlinearity aspect of this phenomenon has already been studied by Proso et al. [22]. What is more, significant uncertainties already occur at the stage of riveting process, inherently altering the dynamic responses of individual blind-hole rivet specimens. Therefore, in this research a parametric model of the rivet is proposed that is for the process of fatigue parameter assessment individually adjusted to each tested specimen. This manuscript is organised as follows. Section 1 deals with a linear structure's response to a random base-excitation, which is then applied to the vibration fatigue analysis. Additionally, the frequency based cycle counting method is used to obtain the rivet-joint fatigue parameters. In Section 2, first the experiments of the rivet-joint specimens are presented. Afterwards a general numerical model of the specimen is introduced with special attention given to the individualization of each specimen's numerical model. Lastly, according to the experimental results the fatigue parameters of the rivet joint are assessed. Section 3 draws the conclusions. 1 THEORETICAL BACKGROUND A novel method of fatigue-parameter assessment from a random-signal base-excitation testing can be divided into three main steps. Firstly, the power spectral density (PSD) of the structure's stress response to the random-signal base-excitation is obtained. Secondly, the structure's stress response is introduced into the proposed fatigue-life estimation method (Fig. 1). Lastly, by associating the experimentally obtained specimens's fatigue lives with the numerical fatigue-life estimation the representative values of fatigue strength and fatigue exponent are obtained for the set of tested specimens. This section provides the underlying theoretical background of the new fatigue-parameter assessment method. 1.1 Stress Response to Random Base-Excitation In this study linear and viscously damped model is proposed to obtain structure's dynamic response. The excitation of the structure is given in terms of the kinematic profile at fixed degrees-of-freedom (DOF) with the acceleration PSD, denoted by Saa(a). This type of excitation is commonly used for performing standard environmental testing [1], therefore giving the developed method an applicative value. In general, the equilibrium equation for the multi-degree-of-freedom (MDOF) system in the case of the base-excitation can be written as [29]: MX + Cx + Kx = -x0 (t )Mg. (1) Here, M, C and K denote the mass, damping and stiffness matrix, respectively. Vector x represents the displacements for the system's degrees of freedom, vector g is the direction vector and a scalar function x0 (t ) represents the acceleration for the constrained DOFs. x0(t) is assumed to be a Gaussian random stationary process. The aim of the following deduction is to obtain the stress-response PSD to the applied acceleration profile with the structural-modification-using-response-functions (SMURF) method [30], which represents a faster alternative to the harmonic analysis of the numerical model. In order to calculate the stress response using SMURF method one must initially perform a numerical modal analysis of the unconstrained structure. The results of the modal analysis are a vector of natural frequencies ar and modal matrices of the mass-normalised displacements u Фх and stresses u Фа. Here superscripted u denotes the unconstrained structure. Using the results of the modal analysis and by employing a mode-superposition method [31] one can write a displacement uX(a) and stress u ст(ет) response to the force excitation F(a) of the unconstrained structure as: UX(®) = uHxf (®) F(®), 11 o(ffl) = 11Haf (ю) F(ffl), (2) (3) where uHrf (a) and uH af (ю) represent the displacement and stress transmissibility matrices to the force excitation. The frequency is denoted with m throughout the manuscript. The jkth element of the transmissibility matrix [30]: uHCTf (ю) can be written as 'H.f * O) = 1 £ + i(c-c^^Jì-^ ) U i * U i * та, jr та ,kr £ + i (c-WrJì-g ) c (4) where denotes the damping-loss factor of the rth natural frequency, "фст>jr is the jrth element of the modal matrix u Фп and * is a complex conjugation. A similar definition also follows for the jkth element of the transmissibility matrix uHxf (a). With known transmissibility matrices uH at (ю) and uHxf(a) it is possible to obtain the transmissibility matrix Haa (a) of the constrained structure that relates the structure's stress response to the base acceleration excitation as: = H И • ao (». (5) Fig. 1. Outline of the fatigue life T calculation The step of calculating Haa (a) is defined with the leading equation of the SMURF method [30]: Ha(©) = ®2 •[UH(©) • 'Hf-'c©)]. (6) By applying the subscript c only those elements of the matrix uHaf (ю) that define the transmissibility from the constrained DOFs to all the stress DOFs are included in the equation. The subscript cc denotes that only the transmissibility functions between the constrained DOFs are used. Once the transmissibility Haa (a) is retrieved the structure's random-signal response to the random signal base-excitation can be obtained. Let the Saa (a) be a known acceleration PSD matrix, that is usually defined in the shaker-control software. By combining U r=1 + Eqs. (5) and (6) the stress-response PSD matrix can be obtained as [32] and [33]: Saa (ю) = H;a (ю) • Saa (ю) • H Ga (ю)т (7) Once the stress-response PSD matrix to the random base-excitation is obtained for an arbitrary point on the structure with Eq. (7), the fatigue analysis passes over from the whole dynamic structure to an individual point on the structure, whose stress state is further analysed. For the case of the planar stress state superscript (p) the stress-response PSD matrix at a single point on the structure has the form: S*p>) = * c) S„=* (c) S„=* c S*w * c) S*w * (c) sby> * c) s*. * c) s^ * (c) s* (c) . (8) Accordingly, for a three-dimensional stress state the stress-response PSD matrix is of size [6*6] [34]. To perform a fatigue analysis of the multi-axial stress in the frequency domain it is essential to obtain the PSD matrix of an equivalent Von Mises stress [23] and [34] of the whole structure with: Sa_ (ю) = Trace[QSna (ю)]. (9) In Eq. (9) Q is a constant matrix, which is for a planar stress state given as: Q 1 -1/2 0| -1/2 0 0 (10) Using this definition, the equivalent Von Mises stress is a stationary zero-mean Gaussian process [23] and [34]; therefore, the existing frequency methods [7] for fatigue-life calculations can be adopted. Before this, some general properties of the equivalent stress-response PSD for an individual point on the dynamic structure need to be defined. 1.1.1 Characteristics of a Stress-Response PSD Since the damage intensity is calculated with frequency-domain counting methods some characteristics of the stress-response spectrum are given here. The shape of the equivalent stress PSD ( Sa (w), Eq. (9)) can be characterized with a set of spectral moments; for a stationary and a zero mean-valued random process oeq (t) the mth moment is defined as [7]: К = I mmS„ (m) dm. Jo "'I From Eq. (11) it is obvious that the presented calculation is fully performed in the frequency domain as frequency m represents the independent variable. Certain spectral moments of the equivalent stress PSD can be used to describe the key properties of a stress loading in the time domain [7]. The number of zero crossings with a positive slope per unit time is defined with: _L h 2n V h (12) Additionally, the combination of spectral moments gives a number of bandwidth parameters, generally defined as: а = A„ (13) 1.2 Stress Response to Random Base-Excitation When the stress loading within a structure is random, the fatigue life can be estimated in the time or frequency domain [5]. For the case of random base-excitation, the frequency-domain methods are faster and based only on the characteristics of the equivalent stress PSD and on the fatigue properties of the material. Whether the time- or frequency-domain method is used for the fatigue-life estimation, the general rules for fatigue-damage accumulation hold. The most fundamental principle is Palmgren-Miner's linear damage-accumulation rule [7]: D = У ^. ,N, (14) In Eq. (14) D denotes the total fatigue damage; the fatigue failure is recognized when D equals unity. nt denotes the number of cycles at a certain load amplitude and N is the total number of cycles of a certain load amplitude that would cause fatigue failure. In order to obtain N for an arbitrary load amplitude the material fatigue parameters should be considered, as stated in Basquin's equation [35]: o = C N~llb. (15) Here, a denotes the stress amplitude and N the number of total cycles to the fatigue failure. Eq. (15) introduces two material parameters; the fatigue = strength C and the fatigue exponent b. Introducing Eq. (15) into Eq. (14) the damage intensity can be written as [7]: D = 2-C-\p j-y Pa (a) da, (16) where pa (a) denotes the amplitude distribution of the load history. From Eq. (16) it is evident that the damage intensity relates directly to the bth moment of the stress-amplitude distribution. According to Pitoiset and Preumont [34] the equivalent stress can be proposed in Eq. (16). In the time domain the stress distribution pa (a) is obtained from the rainflow matrix of the stress history [5]; however, in the frequency-domain method the bth moment of pa(a) is calculated with the use of statistical characteristics of the equivalent stress PSD. In this study the Tovo-Benasciutti method [36] is used, as it has been proved to give the most accurate results for a number of different random-load profiles [8]. With the Tovo-Benasciutti method the damage intensity is obtained as: D = ( + (1 - B,b )aB2*-1 ) Dnb . (17) Here, Dnb is the damage intensity for the narrowband stress profile and Btb is a factor obtained from the bandwidth parameters of the stress PSD: Dnb =a2vpC-b ((2X0 ) + 2 j, (18) B _ (a a2 ) tb t 1\2 ' a2 - 1) (1.112(1 + a1a2 - (a + a2))e211a + (a1 -a2)), (19) where ax and a2 are the bandwidth parameters, defined in Eq. (13). 1.3 Fatigue-Parameter Assessment The most often used approach for vibration fatigue analysis is firstly to obtain the fatigue parameters on specially prepared samples with classic fatigue-testing machines and later use these fatigue parameters on a numerical model and in equations for the damage intensity. However, when designing a new or modified product, often only the information about the product's fatigue life during accelerated random-vibration tests or service life is available. If the fatigue parameters could be obtained from accelerated vibration tests with random excitation only, the classic fatigue testing would become unnecessary. However, for the vibration fatigue testing the classical specimens cannot be used since they don't exhibit any dynamic response in the shaker's frequency range. Consequently, an alternative testing specimen with adequate dynamic response needs to be designed. In any case, when predicting the vibration fatigue life the vibration testing is an unavoidable step for two reasons: firstly, to verify the numerical model of the structure and, secondly, to obtain the damping ratios of the structure (damping greatly influences the magnitude of stress PSD and can only be identified experimentally). In terms of existing procedures for fatigue parameter assessment of welded joints [9] the new method developed in this study can be regarded as notch stress approach, relating the numerically obtained stress concentration at the fatigue zone to the actual fatigue life of specimen. Drawing further parallels to existing studies related to the presented research, Allegri and Zhang [23] applied inverse power law to obtain S-N curve from experimental results of fatigue testing with random loading [37]. However, the experimental data Allegri and Zhang analyzed was obtained by exciting the specimens with random loading in the frequency range, that was significantly below the specimens first natural frequency. Consequently, no structural dynamic amplification was involved in the stress response leading to the fatigue fracture. To the author's knowledge, the only research dealing with vibration fatigue of rivet joints under random vibration was done by Proso et al. [22] who focused their work on the nonlinear damage accumulation. Compared to the existing studies, the main contribution of the presented research is complete and comprehensive vibration fatigue analysis of randomly excited dynamic structure with complex joint, applied to real experimental data from vibration tests. With it, one can adequately evaluate the vibration testing results and obtain joint fatigue parameters. However, in order to obtain values of fatigue parameters, that can be used in future analyses, the following assumptions were made: linear model of dynamic structure, linear damage accumulation and Gaussian stress load distribution. In this study a Tovo-Benasciutti frequency-domain counting method (Eq. (17)) is implemented in the process of fatigue-parameter assessment. This enables an accurate estimation of the fatigue parameters in which the fatigue tests with different forms of stress-response PSDs can be evaluated altogether. With this, it is possible to establish or broaden a fatigue-parameter database from different accelerated fatigue tests and to predict the vibration fatigue life of a new product, based on the results of earlier vibration fatigue tests of a different structure made from the same material. From Eq. (15) which relates the number of cycles N to the stress amplitude a it is clear that for the case of harmonic load the fatigue parameters b and C directly relate to the slope of the Woehler curve and can easily be obtained with the linear least-squares method [5]. However, when assessing the fatigue parameters with the Tovo-Benasciutti frequency method (Eq. (17)) the damage intensity is not linearly related to different shapes of stress PSD. Therefore the numerical minimization of the sum of the squared differences between the estimated and experimental fatigue lives of the test specimens is used. Different shapes of stress PSD arise due to the uncertainties of riveting process (Sec. 2.1) and scatter of damping ratio between specimens, consequently requiring individualized numerical model for each particular specimen in order to obtain adequate stress concentration values. Here, the function being minimized is written as: AT (b, С ) = X ( log,,, (Tact. ) - log10 (Tes,. (b, C))}2. (20) In Eq. (20) i denotes a single test specimen, Tact is the measured fatigue life of the specimen and Test is the estimated fatigue life of the specimen, based on the numerical model and damage-intensity equation (Eq. (16)). In this study, a simplex algorithm [38] is used for the numerical minimization of AT (b, C). 2 EXPERIMENTAL WORK AND NUMERICAL ANALYSIS In order to show the feasibility of assessing the fatigue parameters with a random base-excitation a blindhole rivet joint was analyzed both experimentally and numerically. Firstly, a series of experiments with different excitation levels was performed on test specimens with a single blind-hole rivet joint, Fig. 2. Secondly, the individualized numerical model was constructed to obtain a valid stress response for each test specimen. Lastly, the fatigue parameters of the rivet joint were obtained with a numerical minimization of the difference between the actual and calculated vibration fatigue lives. 2.1 Random Vibration Test of a Single Rivet Joint In engineering practice the blind-hole rivet joints are used for assembling and attaching sheet metal; usually there are several rivet joints in the whole riveted assembly. From previous testing experience, for the case of a riveted assembly the rivet-joint failure occurs on the edge between the rivet head and the rivet stem, Fig. 3. Here, an isolated single rivet joint is used as a test specimen with the same location of the anticipated fatigue zone, as can be seen from the numerical model in Section 2.2. inertial accelerometer moving base Fig. 2. Experiment setup for the fatigue test of a single blind-hole rivet joint The riveting material is an aluminum alloy A95052 with 2.5 % Mg. The blind hole for the rivet's installation has a diameter of 5.6 mm, and the total height of the rivet is 11.5 mm. The riveted metal strip is made from stainless steel with dimensions of 2 mm X 20 mm X 60 mm. The weight, bolted to the metal strip, weighs 19.1 g and the distance between the weight's and rivet's central axis is 40 mm. Fig. 3. Cross-section of a rivet joint during riveting process For each vibration test the rivet joint was produced with a manual riveting tool by pulling the rivet mandrel. The riveting force was monitored during the riveting process; the maximum measured force was 1800 N. The fatigue test of the rivet joint was performed by exciting the experimental setup with a white noise defined in the frequency area of the first natural frequency from 150 Hz to 350 Hz; outside of defined frequency range the excitation was not present. The PSD amplitude of white noise acceleration profile was monitored and held constant through the whole fatigue test of individual specimen. A total of 16 specimens were tested on a LDS V555 electro-dynamical shaker with constant PSD amplitudes at levels from 30 (m/s2)2/Hz to 100 (m/s2)2/Hz. The fatigue crack initiated on the outer radii on the edge between the rivet stem and the rivet head, which then proceeded through to reach the inner radii of the rivet stem, as can be seen in Fig. 4. The results of the fatigue tests are presented in Fig. 5 in terms of a graph relating the acceleration PSD amplitude to the fatigue life. Fig. 4. Fatigue crack at the rivet-joint Fig. 5. Experimental fatigue lives for AlMg25 rivet joints Apart from the experimental results, presented in Fig. 5, it is also necessary to address the uncertainties occurring within the experiment and to appropriately account for them in the numerical evaluation of the experiment. For the proposed experimental setup the main source of uncertainty presents a process of forming a rivet joint. Here, rivet joint is formed by pulling the rivet mandrel (Fig. 3) through the rivet until the mandrel breaks. When connecting the metal sheet to the base by riveting in the blind hole the pulling length until breaking varied, from 3 mm to 7.5 mm. This alters the geometry and stiffness of the rivet joint, which then changes the initial natural frequency ffl01 and the damping £01 of the tested specimen. As can be seen from Fig. 6, the initial natural frequencies ranged from 266 Hz to 294 Hz. Furthermore, due to the scatter of the pull length the damping and the stress amplitudes at the fatigue zone between the specimens vary as well. The measured data for all the tested specimens are collected in Table 1. This should be taken into account during the numerical evaluation of the experimental results. In order to obtain accurate material parameters for the rivet joint these uncertainties were considered in the numerical evaluation by building several numerical models with different geometries of deformed rivet and modal damping, to obtain an appropriate stress response for each specimen. 300 0.0 0.2 0.4 0.6 0.8 1.0 Damage Fig. 6. Frequency shift of the first natural frequency for the 16 tested specimens Furthermore, during the fatigue test the values of the natural frequencies and the damping ratios change. These changes are a consequence of the growing fatigue crack and can be easily observed from the measurements shown in Figs. 6 and 7. Transmissibility measurements were made with the control and response accelerometer (see Fig. 2) every 20 seconds. By comparing individual frequency shifts (Fig. 6) it is obvious that the natural frequency of the specimens with a shorter fatigue life decreased steadily, as opposed to the specimens with a longer fatigue life, where the decrease of the natural frequency occurred in the last 10 % of the fatigue life. In any case, the natural frequency of the specimens remained within the excitation-frequency band of 150 Hz to 350 Hz; Table 1. Results of rivet-joint fatigue testing and numerical calculation Specimen Initial/си [Hz] Initial [Hz] PSD amplitude [(m/s2)2/Hz] Measured life [s] Calculated life [s] V01 286 0.0140 30 7100 10328 V02 285 0.0142 30 25260 11119 V03 281 0.0177 40 5060 5982 V04 266 0.0167 40 5480 2963 V05 288 0.0134 50 1580 1790 V06 294 0.0139 50 5140 3896 V07 280 0.0164 60 1600 982 V08 284 0.0200 60 1740 3112 V09 279 0.0193 70 820 853 V10 285 0.0165 70 800 809 V11 266 0.0267 80 540 1092 V12 270 0.0254 80 900 1068 V13 281 0.0230 90 680 790 V14 266 0.0251 90 580 683 V15 294 0.0137 100 160 273 V16 266 0.0175 100 260 114 • includes and predicts continuous changes of the structure's modal parameters, • predicts the response of a nonlinear dynamic structure with increasing nonlinearity. The main aim of this study is to present a feasible procedure to assess the fatigue parameters from a random base-excitation that can be easily applied to different MDOF structures without extensive analysis of the phenomena that occur during the fatigue process. Therefore, a linear numerical model of a dynamic structure is proposed, whose linear response matches the initial response that was measured for each specimen. The model's modal parameters are assumed to be constant during the fatigue life. Since the random excitation is broad-band and covering the whole area of the frequency shift, the assumption of a constant natural frequency can be made. Although the damping ratio of the numerical model is constant, a rise of the damping ratio is indirectly included in the final values of the fatigue parameters. The increase of the damping ratio reduces the actual stress response in the fatigue zone. Since the constant damping ratio is included in the numerical model, a more severe stress history is predicted with the numerical model when compared to the actual stress during the fatigue test. By equalizing the fatigue lives of the actual test with a lower stress response with the numerically obtained fatigue lives at a higher stress response (Eq. (20)) a higher value of the assessed fatigue strength is anticipated (usually, the fatigue strength equals the yield strength of the material [5]). therefore, the excitation in the specimen's resonant area was ensured throughout the fatigue test. 20 150 200 250 300 350 /[Hz] Fig. 7. Transmissibility during the fatigue testing of specimen V5 2.2 Numerical Evaluation of Experimental Results As discussed in Sec. 2.1, the specimens have different initial modal properties (natural frequency and damping), which tend to change in different ways during the fatigue test. Additionally, due to the nature of the rivet joint and the growing fatigue crack, the rivet-joint specimens are nonlinear dynamic structures. To accommodate all the stated considerations in the fatigue-parameter assessment a fatigue model should be constructed that: • includes different initial natural frequencies and damping ratios of the individual test specimen, 2.2.1 Numerical Model of a Rivet Joint In order to include the different initial modal parameters of the tested specimens into the numerical analysis a number of different numerical models with different rivet geometries were designed; the rivet geometry is based on the pull length of the rivet mandrel. A shorter pull length results in a thinner rivet wall and a larger distance between the rivet-base contact and the rivet head, which leads to a lower stiffness of the rivet joint. In contrast, a longer pull length results in a thicker rivet wall and a shorter distance between the rivet-base contact and the rivet head, making the rivet joint stiffer. When calculating the stem-wall thickness the volume of the riveting material in the rivet stem was considered to be constant. In this way the first natural frequencies for the numerical model ranged from 265 Hz to 300 Hz. In order to precisely tune the numerical model to the experiment every test specimen was analysed individually to identify the appropriate pull length distance and the damping ratio for the first natural frequency. The comparison between the experimental transmissibilities and the transmissibilities used in the numerical analysis is presented in Fig. 8 for four tested specimens. 20 240 260 280 300 320 Frequency [Hz] Fig. 8. Comparison of initial transmissibilities for selected specimens V02, V03, V06 and V16 (- - - measured, — from numerical model) In Fig. 9 the stress state of the rivet for the specimen's first mode-shape is shown. The maximum equivalent stress for the first mode-shape is located on the edge of the rivet stem and the rivet head, which complies with the actual fatigue-crack initiation (see Sec. 2.1, Fig. 4). The comparison between the actual and the numerically obtained stress state was not carried out in this research due to the inaccessibility of the critical fatigue zone for the test specimen. However, any differences between the values of the actual stress load and the numerically obtained stress load will be indirectly compensated in the assessed fatigue parameters, which hold for the fatigue analysis with the presented numerical model. Therefore, the measurement of the stress state in the fatigue zone is not necessary as long as the validation of the acceleration response of the numerical model to the measured acceleration response has been made. "ft m 1 [iviraj p +2.109e+05 - +1.875e+05 - +1.641 e+05 H \T\TF l \ i \ И f fl.406e+05 "Ml - +9.375e+04 4-7 m 1 f-t-04 - +4.687e+04 K2.344e+04 H.099e-ll III 1 -Ц- --U Fig. 9. Von Mises stress state of first modeshape at the rivet for the pull length of 5.75 mm 2.2.2 Fatigue-Parameter Assessment and Discussion Using the validated numerical model for each tested specimen the fatigue parameters can be assessed. However, since the characteristics of the stress response (Eqs. (11) to (13)) differ from specimen to specimen no linear correlation exists between the damage intensities of the specimens. Therefore, to retrieve the fatigue exponent and the fatigue strength of a rivet joint a numerical minimization with of the function AT (b, C) (Eq. (20)) was proposed using the simplex algorithm. In the process of AT (b, C) minimization the damage intensity for all numerical models was calculated with Tovo-Benasciutti method for given fatigue parameters at certain step of simplex algorithm and compared to the actual fatigue lives. The minimization method showed strong convergence. Using the identified rivet-joint fatigue parameters Basquin's equation (Eq. (15)) can be written as: o= 1748.3 -106 • N~1/752. (21) The value of the minimized function AT (b, C) equals 0.647. Fig. 10 presents the comparison between the actual fatigue lives and the fatigue lives calculated using the fatigue parameters displayed in Eq. (21). The data presented in 10 is also stated in Table 1. From Fig. 10 it is clear that the developed method for assessing the rivet-joint fatigue parameters gives feasible results. Here, it should be noted that the obtained fatigue parameters only hold for the presented numerical model of a rivet joint, which was also used for the minimization of AT (b, C). In the case of different types of rivet-joint modelling (i.e. simplified model, additional fillet between the rivet stem and the rivet head) the obtained fatigue parameters should not be used; hence, the presented procedure for the assessment of the rivet-joint fatigue parameters must be repeated on the modified numerical model to obtain the appropriate fatigue parameters. From Eq. (21) it is evident that the obtained value of fatigue strength C is reasonably high; generally, the fatigue strength C equals the yield strength of the material [5], which is about 200 MPa for aluminum alloys. The high value of the fatigue strength C can be assigned the high stress concentration of the equivalent stress aeq in the numerical model, Fig. 9. From the obtained fatigue-strength parameter it is reasonable to assume that the stress concentration in the numerical model is higher than in the actual specimen. However, as stated in Sec. 2.2.1, the actual stress in the fatigue zone could not be measured due to the inaccessibility of the fatigue zone. An additional and by author's experience the main reason for the high value of the fatigue strength C is the increase of the structural damping during the fatigue test [39]. The identified fatigue exponent b for the aluminum rivet joint shows good correlation to the values obtained by Allegri and Zhang [23] and Česnik in Slavič [39]. 10 К •........, .......; ........< 10 10 10 10 Predicted life [s] Fig. 10. Experimental and calculated fatigue lives of rivet-joint specimens From the results of the fatigue-parameter assessment we can conclude that the method gives relevant results and is feasible for the case of the random excitation of specimens with narrow-band stress response near the single natural frequency. However, the method introduced in this study is, by definition, feasible even if the tested specimens have different stress responses due to different geometries or forms of the excitation profiles. This increases the method's relevance for the case of the multi-modal stress-response spectrum [40] and [41]. The presented method gives an engineer the possibility to obtain the fatigue parameters of the chosen material or design element from the results of different vibration test on different structures. 3 CONCLUSIONS In this study a complete procedure for fatigue-parameter assessment in the case of random base-excitation was presented and applied to the linearized model of a rivet joint. The presented procedure and its experimental verification consist of several key steps: • calculation of the stress response to random base-excitation using the SMURF method, • implementation of a frequency cycle-counting method into the fatigue-parameter assessment, • execution of a series of vibration fatigue tests to obtain experimental data, • construction of a numerical model of the rivet joint that considers the uncertainties of the riveting process, • tuning of a numerical model for a particular test specimen, • minimization of the total deviation between the experimental and calculated fatigue lives as a function of unknown fatigue parameters. In the experimental part of the study a test specimen with an isolated rivet joint was designed in order to reflect the most common rivet-joint fatigue failure and to reduce the number of uncertainties that could arise during the fatigue test. The riveting process and its repeatability were recognized as the main source of uncertainty, which was also taken into consideration during the construction of the numerical model. Regardless of the variance of the riveting process the experimental results show the expected trend of the total fatigue life, similar to the Woehler curve. A detailed numerical model of the rivet joint was constructed that was individually tuned to reflect the actual response of the particular test specimen. The variations of the numerical model with different natural frequencies, damping ratios and geometries were than combined to assess the fatigue parameters as accurately as possible. The proposed method was implemented into the procedure of the fatigue-parameter assessment considering the full shape of the stress-response PSD. For the case of the proposed experiment the methods gives comparable values of fatigue exponent to the related studies. Considering the fatigue strength the method gives the effective value, that inherently includes real-life phenomenon of damping ratio increase during the damage accumulation process. The main contribution of this study is the implementation of the frequency based cycle-counting method into the fatigue parameter assessment and its application to the real experimental results. The obtained fatigue parameters of the rivet joint reflect the actual fatigue phenomenon of the rivet and can be used together with the proposed numerical model for the prediction of the fatigue lives for different structures with single or multiple rivet joints excited with different shapes of random vibration profiles. Additionally, with the presented method it is possible to obtain fatigue parameters that are valid for the proposed numerical model, without direct monitoring of the stress state in the fatigue zone. The obtained fatigue parameters indirectly compensate the rise of the damping ratio during the fatigue process and any discrepancies between the stress-concentration factors for the actual specimen and the numerical model. 4 REFERENCES [1] MIL-STD-810 (2008). Environmental Engineering Considerations and Laboratory Tests, revision G, Method 514.6: Vibration. US Department of the Air Force, Washington D.C. [2] Lalanne, C. (2002). Mechanical Vibration & Shock: Specification Development, Volume V. Taylor & Francis, New York. [3] Lalanne, C. (2002). Mechanical Vibration & Shock Fatigue Damage, Volume IV. Taylor & Francis, New York. [4] Česnik, M., Slavič, J., Boltežar, M. (2012). Uninterrupted and accelerated vibrational fatigue testing with simultaneous monitoring of the natural frequency and damping. Journal of Sound and Vibration, vol. 331, no. 24, p. 5370-5382, D0l:10.1016/j.jsv.2012.06.022. [5] Lee, Y.-L., Pan, J., Hathaway, R. B., Barkey, M.E. (2005). Fatigue Testing and Analysis: Theory and Practice. Elsevier, Oxford. [6] Dirlik, T. (1985). Application of Computers in Fatigue Analysis. Ph.D. thesis, University of Warwick, Warwick. [7] Benasciutti, D. (2004). Fatigue Analysis of Random Loadings. Ph.D. thesis, University of Ferrara, Ferrara. [8] Mršnik, M., Slavič, J., Boltežar, M. (2013). Frequency-domain methods for a vibration-fatigue-life estimation - Application to real data. International Journal of Fatigue, vol. 47, p. 8-17, D0I:10.1016/j.ijfatigue.2012.07.005. [9] Hobbacher, A. (2003). Recommendations for Fatigue Design of Welded Joints and Components. International Institute of Welding, Paris. [10] Grover, H.J. (1967). Fatigue of Aircraft Structures. US Government Printing Office, Washington D.C. [11] Schijve, J. (1994). Fatigue of aircraft materials and structures. International Journal of Fatigue, vol. 16, no. 1, p. 21-32, DOI:10.1016/0142-1123(94)90442-1. [12] ABS (2003). Guide for the Fatigue Assessment of Offshore Structures. American Bureau of Shipping, Houston. [13] Bishop, N.W.M. (1988). The Use of Frequency Domain Parameters to Predict Structural Fatigue. University of Warwick, Warwick. [14] Liu, X., Sooklal, V.K., Verges, M.A., Larson, M.C. (2006). Experimental study and life prediction on high cycle vibration fatigue in BGA packages. Microelectronics Reliability, vol. 46, no. 7, p. 1028-1138, D0I:10.1016/j.microrel.2005.09.011. [15] Chen, Y.S., Wang, C.S., Yang, Y.J. (2008). Combining vibration test with finite element analysis for the fatigue life estimation of PBGA components. Microelectronics Reliability, vol. 48, no. 4, p. 638-644, D0I:10.1016/j.microrel.2007.11.006. [16] Yu, D., Al-Yafawi, A., Nguyen, T.T., Park, S., Chung, S. (2011). High-cycle fatigue life prediction for Pb-free BGA under random vibration loading. Microelectronics Reliability, vol. 51, no. 3, p. 649-656, D0I:10.1016/j.microrel.2010.10.003. [17] George, T.J., Seidt, J., Shen, M.H.H., Nicholas, T., Cross, C.J. (2004). Development of a novel vibration-based fatigue testing methodology. International Journal of Fatigue, vol. 26, no. 5, p. 477-486, D0I:10.1016/j.ijfatigue.2003.10.012. [18] Pagnacco, E., Lambert, S., Khalij, L., Rade, D.A. (2012). Design optimisation of linear structures subjected to dynamic random loads with respect to fatigue life. International Journal of Fatigue, vol. 43, p. 168-177, D0I:10.1016/j. ijfatigue.2012.04.001. [19] Paulus, M., Dasgupta, A. (2012). Semi-empirical life model of a cantilevered beam subject to random vibration. International Journal of Fatigue, vol. 45, p. 82-90, D0I:10.1016/j. ijfatigue.2012.06.008. [20] Han, S.-H., An, D.-G., Kwak, S.-J., Kang, K.-W. (2013). Vibration fatigue analysis for multi-point spot-welded joints based on frequency response changes due to fatigue damage accumulation. International Journal of Fatigue, vol. 48, p. 170177, D0I:10.1016/j.ijfatigue.2012.10.017. [21] Braccesi, C., Cianetti, F., Lori, G., Pioli, D. (2015). Random multiaxial fatigue: A comparative analysis among selected frequency and time domain fatigue evaluation methods. International Journal of Fatigue, vol. 74, p. 107-118, D0I:10.1016/j.ijfatigue.2015.01.003. [22] Proso, U., Slavič, J., Boltežar, M. (2016). Vibration-fatigue damage accumulation for structural dynamics with non-linearities. International Journal of Mechanical Sciences, vol. 106, p. 72-77, D0I:10.1016/j.ijmecsci.2015.12.005. [23] Allegri, G., Zhang, X. (2008). On the inverse power laws for accelerated random fatigue testing. International Journal of Fatigue, vol. 30, no. 6, p. 967-977, D0l:10.1016/j. ijfatigue.2007.08.023. [24] Skorupa, M., Skorupa, A., Machniewicz, T., Korbel, A. (2010). Effect of production variables on the fatigue behaviour of riveted lap joints. International Journal of Fatigue, vol. 32, no. 7, p. 996-1003, DOI:10.1016/j.ijfatigue.2009.11.007. [25] Skorupa, A., Skorupa, M., Machniewicz, T., Korbel, A. (2014). Fatigue crack location and fatigue life for riveted lap joints in aircraft fuselage. International Journal of Fatigue, vol. 58, p. 209-217, DOI:10.1016/j.ijfatigue.2013.01.014. [26] Huang, W., Wang, T.-J., Garbatov, Y., Guedes Soares, C. (2012). Fatigue reliability assessment of riveted lap joint of aircraft structures. International Journal of Fatigue, vol. 43, p. 54-61, DOI:10.1016/j.ijfatigue.2012.02.005. [27] Sun, X., Stephens, E.V., Khaleel, M.A. (2007). Fatigue behaviors of self-piercing rivets joining similar and dissimilar sheet metals. International Journal of Fatigue, vol. 29, no. 2, p. 370-386, DOI:10.1016/j.ijfatigue.2006.02.054. [28] Fu, M., Mallick, P.K. (2003). Fatigue of self-piercing riveted joints in aluminum alloy 6111. International Journal of Fatigue, vol. 25, no. 3, p. 183-189, DOI:10.1016/S0142-1123(02)00115-9. [29] Ewins, D.J. (2000) Modal Testing: Theory, Practice and Application. ResearchStudies Press Ltd., Baldock. [30] Česnik, M., Slavič, J., Čermelj, P., Boltežar, M. (2013). Frequency-based structural modification for the case of base excitation. Journal of Sound and Vibration, vol. 332, no. 20, p. 5029-5039, DOI:10.1016/j.jsv.2013.04.038. [31] Maia, N.M.M., Silva, J.M.M. (1997). Theoretical and Experimental Modal Analysis. Research Studies Press Ltd., Baldock. [32] Newland, D.E. (1984). Random Vibration and Spectral Analysis. Longman Scientific & Technical, Essex. [33] Bonte, M., de Boer, A., Liebregts, R. (2004). Prediction of mechanical fatigue caused by multiple random excitations. Proceedings of International Conference on Noise and Vibration Engineering, p. 697-708. [34] Pitoiset, X., Preumont, A. (2000). Spectral methods for multiaxial random fatigue analysis of metallic structures. International Journal of Fatigue, vol. 22, no. 7, p. 541-550, DOI:10.1016/S0142-1123(00)00038-4. [35] Basquin, O.H. (1910). The exponential law of endurance tests. Proceedings of American Society of Testing Materials, vol. 10, p. 625-630. [36] Benasciutti, D., Tovo, R. (2005). Spectral methods for lifetime prediction under wideband stationary random processes. International Journal of Fatigue, vol. 27, no. 8, p. 867-877, DOI:10.1016/j.ijfatigue.2004.10.007. [37] Clevenson, S.A., Steiner, R. (1967). Fatigue life under random loading for several power spectral shapes, NASA-TR-266, NASA Langley Research Center, Hampton. [38] Nash, S.G. (1990). A History of Scientific Computing. ACM, New York. [39] Česnik, M., Slavič, J. (2014). Vibrational fatigue and structural dynamics for harmonic and random loads. Strojniški vestnik -Journal of Mechanical Engineering, vol. 60, no. 5, p. 339-348, DOI:10.5545/sv-jme.2014.1831. [40] Low, W.M. (2010). A method for accurate estimation of the fatigue damage induced by bimodal processes. Probabilistic Engineering Mechanics, vol. 25, no. 1, p. 75-85, DOI:10.1016/j.probengmech.2009.08.001. [41] Gao, Z., Moan, T.M. (2008). Frequency-domain fatigue analysis of wide-band stationary Gaussian processes using a trimodal spectral formulation. International Journal of Fatigue, vol. 30, no. 10-11, p. 1944-1955, DOI:10.1016/j. ijfatigue.2008.01.008. Vsebina Strojniški vestnik - Journal of Mechanical Engineering letnik 62, (2016), številka 7-8 Ljubljana, julij-avgust 2016 ISSN 0039-2480 Izhaja mesečno Razširjeni povzetki (extended abstracts) Xiaoning Liu, Gengkai Hu: Pregled uporabe kiralnosti pri elastičnih metamaterialih SI 71 Dezhi Wang, John E. Mottershead: Natančnost meritev in prostorska ločljivost metode Kriging korelacije digitalnih posnetkov SI 72 Yang Liu, Sheikh Islam, Ekaterina Pavlovskaia, Marian Wiercigroch: Optimizacija sistema vibroudarne kapsule SI 73 Ioannis P. Mitseas, Ioannis A. Kougioumtzoglou, Pol D. Spanos, Michael Beer: Določitev verjetnosti preživetja nelinearnega sistema MDOF v pogojih evolucijskega stohastičnega vzbujanja SI 74 Fabian M. Gruber, Daniel J. Rixen: Vrednotenje tehnik za redukcijo podstruktur s fiksnimi in prostimi stiki SI 75 Lionel Manin, Stefan Braun, Damien Hugues: Modeliranje interakcij med trakom in ovojnico pri transportu pošte skozi sortirni stroj SI 76 Martin Česnik, Janko Slavič, Miha Boltežar: Pridobitev parametrov utrujanja iz vibracijskih preizkusov z naključnim signalom: Aplikacija na kovičenem spoju SI 77 Osebne vesti SI 78 Pregled uporabe kiralnosti pri elastičnih metamaterialih Xiaoning Liu - Gengkai Hu Državni laboratorij za dinamiko in krmiljenje zračnih plovil, Ministrstvo za šolstvo, Šola za letalsko in vesoljsko tehniko, Institut za tehnologijo v Pekingu, Kitajska Metamateriali so umetno zasnovani kompozitni materiali z neobičajnimi fizikalnimi lastnostmi, kakršnih običajno ni mogoče najti v naravi. Takšne lastnosti se navadno pojavijo v okolju z valovanji, ko nastopi pojav lokalne resonance. Načela snovanja metamaterialov se torej zanašajo na ustvarjanje primernih lokalnih resonanc in na njihove interakcije z valovi v ozadju. V tem pregledu smo pokazali, da je s povezavo med vrtenjem in deformacijami dvodimenzionalnih kiralnih trdnih teles mogoče doseči edinstven resonančni mehanizem za snovanje elastičnih metamaterialov, podan pa je tudi pregled novosti na tem področju. Čeprav so že dolgo znani različni mehanizmi za pripravo elastičnih metamaterialov, pa je snovanje praktičnih elastičnih metamaterialov s kombiniranimi eksotičnimi lastnostmi še vedno težaven problem zaradi kompleksnosti elastičnih valov. Prispevek začnemo z metamaterialom negativne gostote, ki je ustvarjen z integracijo kiralne rešetke in resonančnih vključkov. Ta metamaterial ni le primeren za blaženje širokopasovnih vibracij, temveč zaradi kiralnosti izkazuje tudi mešano translacijsko-rotacijsko resonanco. Tovrstna mešana resonanca je še dodatno preučena za realizacijo dvakrat negativnega metamateriala na podlagi čistih trdnih materialov. Podan je tudi predlog zasnove praktičnega elastičnega metamateriala, ki ima istočasno negativno gostoto in modul elastičnosti ter je narejen iz enega samega materiala. Zasnova je bila eksperimentalno potrjena po tehniki plošč s pojavom negativnega loma. Za ustreznejšo karakterizacijo dvodimenzionalnih kiralnih elastičnih materialov in metamaterialov so predstavljene tudi novosti na področju kiralnih mikropolarnih konstitutivnih modelov. Za izračun pasovne strukture metamateriala s kompleksno konfiguracijo je bila uporabljena Blochova valovna analiza v povezavi z metodo končnih elementov. Dejanske dinamične lastnosti so bile ovrednotene s s tehniko makro-mikro prehoda, izposojeno iz tradicionalne mikromehanike, s harmonično obremenjenimi robnimi pogoji. Vsi demonstracijski diskretni modeli so formulirani analitično. Opravljene so bile tudi polne harmonične in tranzientne valovne simulacije negativnega loma s komercialno programsko opremo COMSOL in LS-DYNA za analizo po metodi končnih elementov. Vzorec z negativnim lomom je bil izdelan s tehniko laserskega rezanja velike plošče iz nerjavnega jekla, dvodimenzionalni val pa je bil posneman z najnižjo zvrstjo Lambovega vala. Valovni signal je bil ustvarjen s piezoelektričnim vzbujanjem in merjen z laserskim vibrometrom. Pri razvoju dvodimenzionalne kiralne mikropolarne konstitutivne relacije je bila uporabljena predstavitev s tenzorjem, ki ga ni moč reducirati. Z integracijo resonatorjev z negativno gostoto v kiralno rešetko smo prikazali zmožnost dušenja valov v širokem pasu z enkrat negativnim elastičnim metamaterialom. Zasnovali smo dvakrat negativen elastičen metamaterial s praktično preprosto enotsko celico ter prvič eksperimentalno prikazali negativen lom pri elastičnem metamaterialu. Razvit je mikropolarni konstitutivni model dvodimenzionalnih kiralnih trdnih snovi in uporabljen pri materialih s kiralno rešetko. Translacijsko-rotacijski resonančni mehanizem je trenutno še omejen pri dvodimenzionalnih kiralnih trdnih snoveh. Ostaja razširitev predlagane sheme na tridimenzionalni primer. Dobre možnosti snovanja pri kiralnih trdnih snoveh bi bilo mogoče dodatno raziskati za izboljšanje funkcije elastičnih metamaterialov. S pomočjo kiralnosti in povezanih resonanc bi lahko občutno poenostavili zasnovo dvakrat negativnih elastičnih metamaterialov. Uvedba kiralnosti pa terja še bolj dovršeno teorija kontinuuma, npr. teorijo mikropolarnosti, za boljše razumevanje in karakterizacijo tega materiala. Avtorji upajo, da bo predstavljena raziskava pomagala pri uporabi elastičnih metamaterialov z naprednimi valovi. Ključne besede: razširjanje valov, dejanska lastnost, elastični metamateriali, kiralne trdne snovi, negativen lom, teorija mikropolarnosti *Naslov avtorja za dopisovanje: 5 South Zhongguancun ulica, Haidian okrožje, Peking, Kitajska. liuxn@bit.edu.cn SI 71 Natančnost meritev in prostorska ločljivost metode Kriging korelacije digitalnih posnetkov Dezhi Wang1 - John E. Mottershead12* 1 Univerza v Liverpoolu, Središče za tehnično dinamiko, Združeno kraljestvo 2 Univerza v Liverpoolu, Institut za tveganja in negotovost, Združeno kraljestvo Korelacija digitalnih posnetkov (DIC) je danes verjetno najbolj priljubljena metoda za merjenje odmikov in deformacij na področju eksperimentalne mehanike. Algoritme DIC lahko razvrstimo v dve glavni kategoriji: na lokalne algoritme DIC, ki bazirajo na podmnožicah, in globalne algoritme DIC, ki delajo s polnim poljem. Vrednotenje učinkovitosti različnih algoritmov DIC z ozirom na natančnost meritev in prostorsko ločljivost po standardnem postopku je pomembno za praktično uporabo. Idealni algoritem DIC bi v istem trenutku dosegal tako odlično natančnost kakor tudi odlično prostorsko ločljivost, toda v praksi so vedno potrebni določeni kompromisi. Natančnost meritev je mogoče kvantificirati z vrednotenjem standardnega odklona merilnega šuma. Za preučevanje prostorskih fluktuacij polja odmikov je običajno treba določiti prostorsko ločljivost, kar pa je zelo težavna naloga, če k njej pristopamo eksperimentalno. V izogib tem omejitvam so bili uporabljeni sintetični posnetki peg z naloženimi ravninskimi sinusnimi odmiki. Za učinkovito statistično preiskavo prostorske ločljivosti so bile uporabljene različne prostorske frekvence. Na ta način je bila redefinirana tudi prostorska ločljivost lokalnih in globalnih algoritmov DIC. Z istočasnim upoštevanjem natančnosti in prostorske ločljivosti je zagotovljena verodostojna primerjava učinkovitosti različnih algoritmov DIC. Avtorja sta nedavno razvila nov globalni algoritem DIC, v katerega je vgrajena regresija Kriging. Regresijski model Kriging je integriran v globalni okvir DIC kot funkcija oblike v polnem polju za točnejše ocenjevanje neznanega polnega polja odmikov. Dejansko polje odmikov je bilo modelirano z Gaussovim naključnim procesom, rezultat pa je model Kriging kot najboljša linearna nepristranska napoved. Merilna negotovost v kontrolnih točkah Kriging je bila upoštevana s tehniko globalne regularizacije. Predlagani pristop Kriging DIC daje najverjetnejšo meritev namesto krivulje, ki se prilega merilnim točkam. Pristop se odlično izkaže pri ločljivosti odmikov, pri delu z odprtimi eksperimentalnimi podatki pa je učinkovitejši od mnogih drugih metod DIC v polnem polju. V predstavljeni študiji je po izpopolnjenem standardnem postopku ovrednotena učinkovitost metode Kriging DIC v primerjavi s klasično metodo DIC na podlagi podmnožic. Za kvantifikacijo natančnosti meritev odmika in prostorske ločljivosti algoritmov DIC so bili uporabljeni sintetični posnetki z naloženimi ravninskimi sinusnimi polji odmikov različnih prostorskih frekvenc. Predlagani sta dve novi merili neujemanja odmikov za redefiniranje natančnosti meritev in prostorske ločljivosti lokalnih in globalnih metod DIC. Meritev razlik med dejansko sinusoido in meritvijo je najbolj točna tedaj, ko poteka vzdolž normale v kontrolnih točkah Kriging. Tako se lahko bolje ocenijo relativne prednosti različnih algoritmov DIC. Prostorska ločljivost je polovica najnižje prostorske periode v slikovnih točkah, kjer neujemanje meritev DIC (RMSE) doseže prag, običajno 5 % ali 1,5 % amplitude sinusoide. Natančnost je definirana kot standardni odklon (pri neujemanju 5 % in 1,5 %), merjen v slikovnih točkah v smeri normale ali v smeri osi Y. Na koncu so rezultati za ovrednotenje relativne učinkovitosti različnih pristopov DIC prikazani še v grafični obliki. Izkaže se, da je pristop Kriging DIC robusten glede merilnega šuma ter mnogo učinkovitejši od klasičnega pristopa DIC na bazi podmnožic, tako glede natančnosti meritev odmika kot glede prostorske ločljivosti. Rezultati ocenjevanja učinkovitosti metode DIC so najboljši, če se neujemanje meri v normalni smeri in ne v smeri osi Y. Ključne besede: korelacija digitalnih posnetkov, regresija Kriging, merilna natančnost, prostorska ločljivost SI 72 *Naslov avtorja za dopisovanje: Univerza v Liverpoolu, Središče za tehnično dinamiko, L69 3GH, Združeno kraljestvo, j.e.mottershead@liv.ac.uk Optimizacija sistema vibroudarne kapsule Yang Liu1* - Sheikh Islam1 - Ekaterina Pavlovskaia2 - Marian Wiercigroch2 1 Univerza Roberta Gordona, Tehniška šola, Združeno kraljestvo 2Univerza v Aberdeenu, Center za aplikativne raziskave dinamike, Združeno kraljestvo Namen članka je opis optimizacije premočrtnega napredovanja sistema vibroudarne kapsule. Najprej so bile z eksperimentalno preiskavo na novem preizkuševališču ugotovljena najhitrejša in najbolj učinkovita napredovanja kapsule. Optimizirana sta bila krmilna parametra amplituda in frekvenca harmoničnega vzbujanja ter eden od sistemskih parametrov, t. j. razmerje togosti. V drugem koraku je bila preučena optimizacija oblike kapsule po metodi simulacije z računalniško dinamiko fluidov. Za najboljše napredovanje so bile zmanjšane na minimum sile upora in vzgona, ki delujejo na stacionarno kapsulo v cevi, skozi katero se pretaka fluid. Pri reševanju problema je bil uporabljen pristop eksperimentalne raziskave in simulacije z računalniško dinamiko fluidov (CFD). Tematsko področje članka je dinamika in krmiljenje za optimizacijo sistema vibroudarne kapsule. Primerjave eksperimentalnih rezultatov povprečnega napredovanja in učinkovitosti porabe energije potrjujejo, da krmilni parametri za najhitrejše napredovanje niso tudi najbolj učinkoviti. Sistem kapsule lahko tako deluje hitro ali pa energijsko učinkovito, odvisno od vsakokratnih zahtev. Rezultati CFD kažejo, da so sile upora in vzgon odvisni od dolžin kapsule in loka. Končno je bil pridobljen optimalen nabor dolžin kapsul in lokov. Optimalni geometrijski parametri kapsule za najboljše napredovanje so bili določeni z minimizacijo sil upora in vzgona v simulaciji CFD, pri čemer je bila uporabljena stacionarna »idealna« kapsula v cevi, skozi katero se pretaka fluid. Privzeli smo zaprto preizkuševališče in v simulacijah CFD je bila upoštevana samo zunanja kapsula. V nadaljnjih raziskavah bi bilo mogoče s parametričnimi preiskavami te optimalne zasnove s tridimenzionalnim modeliranjem in časovno odvisnimi simulacijami pridobiti celovito razumevanje interakcij med kapsulo in fluidom ter natančno napovedati delovanje kapsule v realnih scenarijih. Naslednji korak modeliranja za optimizacijo oblik bi bil upoštevanje premikajoče se kapsule v toku fluida. Članek obravnava dva optimizacijska problema: izbiro parametrov vzbujalne sile za pogon kapsule ter izbiro oblike kapsule za najmanjšo silo upora. Cilj v obeh primerih je zagotavljanje najboljšega napredovanja sistema kapsule. Novost pri tem sistemu vibroudarne kapsule je v tem, da ni potreben noben zunanji pogonski mehanizem, zato je lahko sistem zaprt in se premika neodvisno po kompleksnem okolju. Ugotovitve v tem članku so pomembne za snovanje in izdelavo prototipov v poljubnem povečanem ali pomanjšanem merilu. Ključne besede: sistem kapsule, vibroudarec, eksperiment, optimizacija, simulacija CFD, učinkovitost, pregled cevi *Naslov avtorja za dopisovanje: Robert Gordon University, School of Engineering, Garthdee Road, AB10 7GJ, Aberdeen, Združeno kraljestvo, y.liu8@rgu.ac.uk SI 73 Določitev verjetnosti preživetja nelinearnega sistema MDOF v pogojih evolucijskega stohastičnega vzbujanja Ioannis P. Mitseas14* - Ioannis A. Kougioumtzoglou2 - Pol D. Spanos3 - Michael Beer145 1 Leibnizova univerza v Hannovru, Fakulteta za gradbeništvo in geodezijo, Nemčija 2 Univerza Columbia, Oddelek za gradbeništvo in tehnično mehaniko, ZDA 3 Univerza Rice, Katedra L. B. Ryona za inženirstvo, Oddelek za strojništvo, ZDA 4 Univerza v Liverpoolu, Institut za tveganja in negotovost in Tehniška šola, Združeno kraljestvo 5 Univerza Tongji, Šola za gradbeništvo in Institut v Šanghaju za preprečevanje in pomoč ob nesrečah, Kitajska Naravne sile, ki vzbujajo dinamične sisteme, kot so npr. veter, valovi in seizmične obremenitve, so pogosto evolucijske narave. V takšnih pogojih se močno spreminja tako intenzivnost vzbujanj kakor tudi njihova frekvenčna vsebina. Ta razred strukturnih obremenitev je zato treba popisovati z nestacionarnimi stohastičnimi procesi. Močno vzbujani konstrukcijski sistemi lahko izkazujejo tudi značilno nelinearno vedenje s histerezami. Za skupnost, ki se ukvarja z dinamiko konstrukcij, je zato zelo pomemben razvoj tehnik za določanje odziva in vrednotenje zanesljivosti nelinearnih/histereznih sistemov, izpostavljenih evolucijskim stohastičnim vzbujanjem. Zanesljivostna analiza dinamičnega sistema je tesno povezana z določitvijo verjetnosti, da bo odziva sistem ostal znotraj predpisanih meja v določenem časovnem intervalu. Določitev omenjene časovno odvisne verjetnosti, ki je znana tudi kot verjetnost preživetja, že desetletja predstavlja izziv na področju stohastične dinamike. Alternativna definicija tega problema, znana tudi kot problem prvega prehoda, je vrednotenje verjetnosti, da bo odziv sistema prvič presegel vnaprej določen prag v danem časovnem intervalu. Eno najmočnejših orodij za vrednotenje zanesljivosti sistema so simulacije Monte Carlo. Obstajajo pa primeri, ko je računska zahtevnost teh tehnik prevelika, še posebej pri velikih in kompleksnih sistemih. Zato obstaja potreba po razvoju alternativnih učinkovitih aproksimativnih analitičnih/numeričnih tehnik za reševanje problema prvega prehoda. Razvita je bila aproksimativna analitična tehnika za določanje časovno spremenljive verjetnosti preživetja in povezane funkcije gostote verjetnosti prvega prehoda (PDF) za nelinearen/histerezni sistem z več prostostnimi stopnjami (MDOF) v pogojih evolucijskega stohastičnega vzbujanja. Na podlagi učinkovitega pristopa za zmanjšanje števila dimenzij ter konceptov stohastičnega povprečenja in statistične linearizacije je bil izvirni nelinearni sistem z n prostostnimi stopnjami razstavljen in preveden v n linearnih časovno spremenljivih (LTV) oscilatorjev s po eno efektivno eno prostostno stopnjo (SDOF), ki ustrezajo posameznim izhodiščnim prostostnim stopnjam. Definirana in izračunana sta časovno spremenljiva dejanska togost in dušenje za vsako posamezno prostostno stopnjo, amplitude PDF za nestacionarni marginalni, tranzitorni in skupni odziv pa so bile učinkovito določene z izrazi zaprte oblike. Končno so bile na računsko učinkovit način določene še PDF za verjetnost preživetja in prvi prehod sistema MDOF. Razvita tehnika je vsestranska, saj je z njo mogoče obravnavati razna nelinearna vedenja ter razna stohastična vzbujanja s poljubnimi neločljivimi oblikami evolucijskega spektra moči (EPS), ki izkazujejo močno spremenljivost intenzivnosti in frekvenčne vsebine. V numeričnih primerih je vključen tudi nelinearen sistem s tremi prostostnimi stopnjami po modelu Bouc-Wen. Zanesljivost tehnike dokazujejo primerjave z relevantnimi simulacijami Monte Carlo. Ključne besede: problem prvega prehoda, nelinearna stohastična dinamika, evolucijski stohastični procesi, nelinearni/histerezni sistemi, verjetnost preživetja, statistična linearizacija SI 74 *Naslov avtorja za dopisovanje: Leibnizova univerza v Hannovru, Fakulteta za gradbeništvo in geodezijo, 30167, Hannover, Nemčija, mitseas@irz.uni-hannover.de Vrednotenje tehnik za redukcijo podstruktur s fiksnimi in prostimi stiki Fabian M. Gruber* - Daniel J. Rixen Tehniška univerza v Münchnu, Institut za aplikativno mehaniko, Nemčija Tehnike za redukcijo podstruktur so učinkovita metoda za zmanjševanje velikih modelov, ki se uporabljajo pri analizi dinamičnega vedenja kompleksnih struktur. Najbolj priljubljen pristop je Craig-Bamptonova metoda (1968) s fiksnim stikom, ki temelji na nihajnih oblikah s fiksnim stikom in oblikah omejitev stika. Metode s prostim stikom, kot sta npr. MacNealova (1971) in Rubinova metoda (1975), nasprotno uporabljajo nihajne oblike prostega stika skupaj z oblikami pritrditve. Podstrukture se pri omenjenih metodah sestavljajo z odmiki stikov (glavni sestav). Dvojna Craig-Bamptonova metoda (2004) uporablja enake sestavine kot zgoraj omejeni metodi s prostimi stiki, le da se podstrukture sestavljajo s pomočjo sil na stiku (dvojni sestav). Ta metoda uveljavlja le šibko združljivost stikov med podstrukturami, s čimer se izognemo težavam z zaklepanjem stikov, ki se včasih pojavljajo pri pristopih s primarnimi sestavi in pri uporabi oblik prostega stika. Dvojna Craig-Bamptonova metoda daje v primerjavi z drugimi metodami prostega stika preprostejše reducirane matrike. Te reducirane matrike so redko posejane, podobne kot klasične Craig-Bamptonove matrike. V tem prispevku obravnavamo primarno (klasično) formulacijo Craig-Bamptonove metode, MacNealovo metodo, Rubinovo metodo in dvojno formulacijo Craig-Bamptonove metode. Izpeljava je podana na celovit in dosleden način. V članku so na kratko predstavljeni in primerjani splošni koncepti Craig-Bamptonove metode (CBM), MacNealove metode (MNM), Rubinove metode (RM) in dvojne Craig-Bamptonove metode (DCBM). Metode so obravnavane na treh primerih. Uporabljene so bile štiri različne metode podstrukturiranja. Predstavljena je formulacija primarne (klasične) Craig-Bamptonove metode, MacNealove metode, Rubinove metode in dvojna formulacija Craig-Bamptonove metode. Izpeljava je podana na celovit in dosleden način. Dvojna Craig-Bamptonova metoda daje boljše rezultate kot CBM z enakim številom normalnih oblik na podstrukturo kot osnovo za redukcijo, s primerljivim obsegom računanja in s podobno redko posejanostjo reduciranih matrik. Primerjava metod s prostim stikom pokaže, da se Rubinova metoda odreže nekoliko bolje kot dvojna Craig-Bamptonova metoda, nastale matrike pa so polne. Rubinova in dvojna Craig-Bamptonova metoda zagotavljata bistveno boljšo točnost aproksimacije kot MacNealova metoda, medtem ko slednja vedno ustvari popolnoma sklopljene reducirane matrike. Potrebne so dodatne raziskave negativnih lastnih vrednosti ter njihovega vpliva na nadaljnje korake po redukciji, če ni opravljeno filtriranje negativnih lastnih vrednosti. Opisane so lastnosti dvojne Craig-Bamptonove metode in ilustriran je dodaten korak filtriranja negativnih lastnih vrednosti med postopkom redukcije. Negativne lastne vrednosti reduciranega dvojnega sestavljenega problema, ki niso fizikalne narave, so neločljivo povezane s procesom redukcije po dvojni Craig-Bamptonovi metodi zaradi šibke združljivosti na stiku med podstrukturami. Filtriranje teh negativnih lastnih vrednostih je odločilni dejavnik odlične kakovosti aproksimacije dvojne Craig-Bamptonove metode. Dodatna računska obremenitev zaradi tega dodatnega koraka je zanemarljiva z ozirom na učinkovitost metode. Ključne besede: dinamično podstrukturiranje, sinteza oblik komponent, redukcija reda modela, dvojni sestav, Craig-Bamptonova metoda, metoda prostega stika *Naslov avtorja za dopisovanje: Tehniška univerza v Münchnu, Institut za aplikativno mehaniko, Boltzmannstr. 15, 85748 Garching, Nemčija, fabian.gruber@tum.de SI 75 Modeliranje interakcij med trakom in ovojnico pri transportu pošte skozi sortirni stroj Lionel Manin1* - Stefan Braun1 - Damien Hugues2 1 Univerza v Lyonu, INSA-Lyon, Francija 2 Solystic SAS, Francija Predstavljena raziskava obravnava mehaniko transporta pošte skozi sortirne stroje v poštnih centrih. Za robustno delovanje teh strojev mora biti zagotovljen pravilen transport pisem, to pa pomeni stabilen položaj pisemskih ovojnic glede na transportni trak. Stroj lahko obdela do 53.000 kosov pošte na uro, odvisno od njihove dolžine. V okviru te raziskave je bilo razvito simulacijsko orodje, ki napoveduje trajektorijo pisma pri transportu med dvema ravnima trakovoma in pri prečkanju več usmerjevalnih kretnic, kjer je stična površina med ovojnico in trakom zmanjšana. Vsaka od dveh strani ovojnice je v stiku z enim trakom. Koeficient trenja na obeh straneh se lahko razlikuje, npr. pri razglednicah ali pri ovojnicah s plastičnim okencem. Trakova naj bi potovala vodoravno in z enako hitrostjo, v resnici pa lahko prihaja do manjših hitrostnih razlik. Opaziti je mogoče tudi odstopanja v vodoravnosti trakov in vsi ti pojavi lahko privedejo do nepričakovanih situacij, ko se relativni položaj ovojnice glede na trak spreminja med transportom. Delo je osredotočeno na mehanske interakcije med trakom in ovojnico za napovedovanje položaja med transportom. Razvit je bil 2D-model, ki upošteva dinamično ravnotežje ovojnice in simulira prehodne pojave na različnih mestih v sortirnem stroju. Ovojnico opredeljujejo njena višina hp in širina lp, masa mp in vztrajnost Ip, položaj in velikost (hf x lf) morebitnega plastičnega okenca, koordinate masnega središča (xG, yG), koordinate levega spodnjega vogala (xOp,yOp) in začetni kot glede na vodoravno os qp. Trak je opredeljen s koordinatami dveh spodnjih vogalov in s širino. Iz vhodnih geometrijskih parametrov je izpeljan kot med trakovi in vodoravno osjo. V modelu dinamičnega trenja je treba določiti stično površino med stranicami ovojnice in trakovi za vsak časovni korak posebej. Sila trenja med trakom in ovojnico deluje v središču teh stičnih površin. Privzeto je, da ovojnico na eni strani drži na traku sila lepenja, druga stran ovojnice pa drsi po traku. Izračunan je položaj težišča ovojnice in kot glede na vodoravno os. Stična površina mora biti določena za vsak korak procesa simulacije, kar je ključno pri tem modelu. Kontaktne površine S predstavljajo mnogokotniki Pt. Gre za presek med trakovi in dvema stranema ovojnice, pri čemer je na eni strani plastično okence. Vsi ti deli so mnogokotniki, ki jih določajo koordinate njihovih štirih oglišč. Predlagana metoda uporablja presečne mnogokotnike in koordinate njihovih točk za določitev stičnih površin in njihovih masnih središč. Za določitev presečnih mnogokotnikov smo izkoristili tehniko prirezovanja mnogokotnikov iz računalniške grafike. V naši implementaciji je uporabljen algoritem prirezovanja Vatti, ki omogoča določitev preseka med mnogokotniki poljubnih oblik v dvorazsežnostnem prostoru. Za vrednotenje predlaganega dinamičnega modela je bilo opravljenih več simulacij. Pri tem so bila upoštevana opažanja iz realnih strojev za sortiranje pošte in intuitivno predvidevanje razvoja položaja ovojnice v različnih pogojih. Grafični vmesnik omogoča vizualizacijo zaporednih položajev ovojnice. Ključne besede: stroj za sortiranje pošte, ravni trakovi, ovojnica, stik med trakom in ovojnico, transport, simulacija SI 76 *Naslov avtorja za dopisovanje: Univerza v Lyonu, INSA LyonCNRS UMR5259, LaMCoS, F-69621, Francija, lionel.manin@insa-lyon.fr Pridobitev parametrov utrujanja iz vibracijskih preizkusov z naključnim signalom: Aplikacija na kovičenem spoju Martin Česnik - Janko Slavič - Miha Boltežar* Univerza v Ljubljani, Fakulteta za strojništvo, Slovenija Vibracijsko utrujanje dinamske strukture je, poleg neželenega pojava prekomernih vibracij in hrupa, poglavitni razlog za izvajanje numeričnih in eksperimentalnih dinamskih analiz izdelka. Dinamski odziv sistema vnese v izdelek določeno napetostno stanje; v primeru neugodnih dinamskih lastnosti izdelka lahko že nizke vibracijske obremenitve hitro privedejo do odpovedi izdelka zaradi utrujanja materiala. Eksperimentalno preverjanje ustreznosti izdelka je ustaljen postopek vibracijskega testiranja prototipov, v zadnjem času pa se veliko raziskovalnega dela usmerja v numerično analizo vibracijskega utrujanja. Namen predstavljene raziskave je vzpostavitev metodologije, ki bo na podlagi eksperimentalnih rezultatov vibracijskih preizkusov omogočila pridobitev parametrov utrujanja za obravnavano strukturo oziroma material. Za izračun pričakovane življenjske dobe izdelka pri vibracijskem utrujanju je potrebno poznati tip in velikost vibracijske obremenitve ter materialne parametre utrujanja. Slednje materialne karakteristike so redko na razpolago, saj se jih pridobi le z obširnim namenskim preizkušanjem. Z novo metodologijo, predstavljeno v članku, se je pridobilo materialne parametre zgolj iz rezultatov vibracijskega testiranja s kinematskim vzbujanjem z naključnim signalom. Zaradi takšnega vzbujanja je razvita metodologija primerna tudi za dinamske strukture z nelinearnimi spoji, saj je v tem primeru odziv strukture bolj stabilen v primerjavi z do sedaj uporabljanim harmonskim vzbujanjem. V sklopu raziskave se je metodologijo uporabilo na kovičenem spoju, s katerim je pločevina pritrjena na osnovni material. V primerjavi z obstoječimi pospešenimi preizkusi utrujanja je pri razviti metodologiji eksperimentalna izvedba bistveno enostavnejša, po drugi strani pa je potrebna obsežnejša numerična analiza za končno pridobitev parametrov utrujanja. V članku je najprej predstavljen pristop k analizi dinamskega odziva strukture pri kinematskem vzbujanju. Le-ta temelji na metodi SMURF, s katero se preko procesa strukturne modifikacije in principa superpozicije lastnih oblik izračuna odziv na predpisan pomik vpetja. V nadaljevanju se je izračunan odziv napetosti v frekvenčnem prostoru povezalo z akumulacijo poškodbe preko metode Tovo-Benasciutti, ki posledično poda tudi oceno o pričakovani življenjski dobi strukture. V naslednjem koraku se je izvedlo vibracijske preizkuse na pripravljenih vzorcih s kovičenim spojem. Iz eksperimentalnih rezultatov se je izkazalo, da dinamski odziv posameznega vzorca močno zavisi od končne geometrije kovičenega spoja, ki se, zaradi naključne narave procesa kovičenja, med posameznimi vzorci precej razlikuje. Posledično se je numerični model parametriziralo do te mere, da se je za vsak eksperimentalni vzorec izdelalo individualni numerični model, pri katerem je bilo ujemanje med izmerjenim in izračunanim odzivom zadovoljivo. Med samimi vibracijskim preizkusom je prišlo tudi do pričakovanega padca lastne frekvence in porasta dušenja, kar pa zaradi načina vzbujanja z naključnim signalom, ni predstavljalo bistvenega vpliva na stabilnost odziva vzorca. Parametri vzbujanja za kovičen spoj so bili pridobljeni z numerično minimizacijo razlike med dejanskimi in pričakovanimi življenjskimi dobami vzorcev. Razvita metoda predstavlja način, kako z uporabo naprednih števnih metod v frekvenčni domeni ustrezno ovrednotiti porušitve izdelkov pri vibracijskih preizkusih. Tak pristop omogoča razvojnim raziskovalcem vzpostavitev baze materialnih parametrov utrujanja za uporabljene materiale ter neposredno vrednotenje alternativnih konstrukcijskih rešitev in materialov iz vidika vibracijskega utrujanja. Ključne besede: preizkus vibracijskega utrujanja, parametri utrujanja, naključno kinematsko vzbujanje, frekvenčne števne metode DOKTORSKE DISERTACIJE, ZNANSTVENA MAGISTRSKA DELA DOKTORSKA DISERTACIJA Na Fakulteti za strojništvo Univerze v Ljubljani je obranil svojo doktorsko disertacijo: • dne 15. junija 2016 Urban TOMC z naslovom: »Izboljšave prenosa toplote aktivnega magnetnega regenerato rja« (mentor: izr. prof. dr. Andrej Kitanovski, somentor: prof. dr. Alojz Poredoš); V doktorskem delu so predstavljene izboljšave prenosa toplote aktivnega magnetnega regeneratorja iz stališča izboljšanja temperaturnega razpona ter obratovalne frekvence naprave. V namen preiskovanja izboljšanja obratovalnega temperaturnega razpona aktivnega magnetnega regeneratorja smo že obstoječ numerični model enoslojnega aktivnega magnetnega regeneratorja nadgradili v model večslojnega regeneratorja. Za potrebe numeričnega modeliranja večslojnega regeneratorja je bilo potrebno implementirati izmerjene magnetokalorične lastnosti La-Fe-Co-Si materialov. Magnetoaklorične lastnosti je bilo potrebno procesirati do takšne mere, da so bile primerne za nadaljno implementacijo v numerični model večslojnega regeneratorja. Z uporabo nadgrajenega numeričnega modela večslojnega regeneratorja smo nato izvedli obširno numerično analizo delovanja večslojnih regeneratorjev. Delovanje večslojnih regeneratorjev smo ovrednotili s stališča doseganja hladilnih moči in hladilnih števil pri različnih temperaturnih razponih. Poleg tega smo izvedli tudi eksperimentalno analizo treh različnih večslojnih regeneratorjev. Eksperimentalni rezultati so služili za validacijo numeričnega modela večslojnega regeneratorja. Drugi del doktorske naloge je namenjen analizi povečanja obratovalne frekvence z uporabo novega koncepta aktivnega magnetnega regeneratorja z implementiranimi termoelektričnimi toplotnimi stikali. V ta namen smo razvili nov numerični model aktivnega magnetnega regeneratorja z implementiranimi termoeletričnimi toplotnimi stikali. Razvit model omogoča simulacije dinamike prenosa ter transporta toplote takšne naprave. Rezultati simulacij nove magnetokalorične naprave s toplotnimi stikali so nam podali enega prvih vpogledov v delovanje takšne naprave. Poleg tega smo nov koncept delovanja termoelektričnih toplotnih stikal tudi eksperimentalno ovrednotili in sicer tako, da smo dokazali njihovo uporabnost v novem konceptu aktivnega magnetnega regeneratorja. * Na Fakulteti za strojništvo Univerze v Mariboru sta obranila svojo doktorsko disertacijo: • dne 7. junija 2016 Tjaša ZUPANČIČ HARTNER z naslovom: »Notranja oksidacija mikrolegiranega zlata za medicinske aplikacije« (mentor: prof. dr. Ivan Anžel); V doktorski disertaciji je obravnavana problematike izdelave mikrolegiranega disperzij sko utrjenega zlata z disperzijo oksidov (La2O3), redko zemeljskega elementa (La) s procesom notranje oksidacije ter študija biokompatibilnosti materiala.. V prvem delu raziskav smo izdelali zlitino Au - 0,5 ut.% La, s postopkom klasičnega litja. Rezultati karakterizacije so pokazali, da že zelo majhna količina mikrolegirnega elementa lantana v zlati zlitini, povzroči tvorbo intermetane faze Au6La, ki pomembno izboljša mehanske lastnosti mikrolegiranega zlata, v primerjavi s čistim zlatom. V drugem delu smo s tehnologijo hitrega strjevanja in valjanja, izdelali trakove z metastabilno mikrostrukturo. S hitrim strjevanjem zlitine Au - 0,5 ut.% La smo dosegli spremembo morfologije strjevanja. Mikrostruktura hitro strjenih trakov sestoji iz zelo finih zrn aAu, s premerom do nekaj mikrometrov ter homogeno razporejeno intermetalno fazo (Au6La), po celotnem preseku traku. Izmerjena trdota trakov je pokazala, da se je v primerjavi z litim stanjem, trdota povečala za skoraj 20 %. Valjane trakove smo rekristalizacijsko žarili in ugotovili, da že majhen dodatek 0,5 % La, zvišuje temperaturo rekristalizacije do 200 K in povzroči rekristalizacijski zastoj. Karakterizacija materiala, je pokazala, da intermetalna faza Au6La v valjanih trakovih tvori neprekinjeno mrežo po mejah zrn, ki dodatno niža hitrost rekristalizacije. Sledila je notranja oksidacija zlitine v litem, valjanem in hitro strjenem stanju. Rezultati notranje oksidacije zlitine v litem ter valjanem stanju so pokazale, da proces notranje oksidacije poteče le po mejah kristalnih zrn z oksidacijo lantana iz intermetalne faze (Au6La). Med tem ko notranja oksidacija hitro strjenih trakov poteče s oksidacijo lantana iz trdne raztopine in tvorbo CNO ter delno po mejah zrn s oksidacijo lantana iz intermetalne faze (Au6La). Na osnovi rezultatov smo identificirali mehanizem notranje oksidacije in postavili model notranje oksidacije zlitine Au- 0,5 ut.% La, ki rezultira v disperzijskem utrjanju mikrolegirane zlate zlitine v površinskem sloju, kljub temu, da je topnost kisika v kristalni mreži zlata zelo majhna. V zaključni fazi smo potrdili, da je material Au-0,5 ut.% La biokompatibilen, kar daje možnost njegove uporabe v obliki tankih nanosov na medicinske implantate; • dne 20. junija 2016 David MOČNIK z naslovom: »Inteligentno toleriranje sklopov glede na tehnološke zmožnosti obdelovalnih postopkov« (mentor: prof. dr. Jože Balič); Sodobna proizvodnjaje podvržena najrazličnejšim zahtevam, ki jih povezuje zahteva po učinkovitosti. Da zadostimo tej zahtevi, je treba tolerance sestavnih delov sklopa načrtovati s premislekom. Z ustreznim načrtovanjem toleranc lahko namreč zelo vplivamo na zmanjšanje proizvodnih stroškov. V ta namen je v doktorski disertaciji za reševanje kompleksnega problema načrtovanja toleranc razvit in predstavljen inteligentni sistem toleriranja, ki s pomočjo metod umetne inteligence, na podlagi vhodnih podatkov porazdeli tolerance sestavnih delov, tako da so stroški izdelave minimalni. Razvita sta dva različna modula za načrtovanje toleranc; modul z optimizacijo z rojem delcev in modul, ki temelji na gravitacijskem iskalnem algoritmu. Uporaba razvitega sistema je prikazana na dveh realnih primerih. Primerjane so vrednosti najnižjih doseženih stroškov s posamezno optimizacijsko metodo, vrednosti predlaganih toleranc in izbrani proizvodni procesi - stroji, ki jih je predlagala inteligenca. Skozi opravljena testiranja pri zasnovi inteligentnega sistema toleriranja se je optimizacija z rojem delcev izkazala za najučinkovitejšo metodo. Razvit je tudi uporabniški vmesnik, ki omogoča enostavno načrtovanje toleranc. V zaključku raziskave je potrjena teza doktorske disertacije, hkrati pa so tudi podane smernice za nadaljnje raziskave. * ZNANSTVENO MAGISTRSKO DELO Na Fakulteti za strojništvo Univerze v Ljubljani je z uspehom zagovarjala svoje magistrsko delo: • dne 22. junija 2016 Suzana DOMJAN z naslovom: »Parametrični modeli za večkriterijsko analizo skoraj nič energijskih stavb« (mentor: prof. dr. Sašo Medved) Na Fakulteti za strojništvo Univerze v Mariboru so z uspehom zagovarjali svoje magistrsko delo: • dne 14. junija 2016 Dean ČERNEC z naslovom: »Analiza posedanja magnetno občutljivih mikrodelcev in določitev magnetne občutljivosti« (mentor: prof. dr. Matjaž Hriberšek); • dne 15. junija 2016 Natalija GORJANC z naslovom: »Čiščenje odpadnih vod iz steklarske industrije« (mentor: prof. dr. Darinka Fakin); • dne 15. junija 2016 Borut SUHADOLNIK z naslovom: »Inovativna tehnična orodja za optimiranje procesne proizvodnje« (mentor: prof. dr. Borut Buchmeister); • dne 17. junija 2016 Borut TURIČNIK z naslovom: »Primerjava konceptov načrtovanja montažnih linij kuhalnih aparatov« (mentor: prof. dr. Borut Buchmeister); Information for Authors All manuscripts must be in English. Pages should be numbered sequentially. The manuscript should be composed in accordance with the Article Template given above. The maximum length of contributions is 10 pages. Longer contributions will only be accepted if authors provide juastification in a cover letter. For full instructions see the Information for Authors section on the journal's website: http://en.sv-jme.eu . SUBMISSION: Submission to SV-JME is made with the implicit understanding that neither the manuscript nor the essence of its content has been published previously either in whole or in part and that it is not being considered for publication elsewhere. All the listed authors should have agreed on the content and the corresponding (submitting) author is responsible for having ensured that this agreement has been reached. The acceptance of an article is based entirely on its scientific merit, as judged by peer review. Scientific articles comprising simulations only will not be accepted for publication; simulations must be accompanied by experimental results carried out to confirm or deny the accuracy of the simulation. Every manuscript submitted to the SV-JME undergoes a peer-review process. The authors are kindly invited to submit the paper through our web site: http://ojs.sv-jme.eu. The Author is able to track the submission through the editorial process - as well as participate in the copyediting and proofreading of submissions accepted for publication - by logging in, and using the username and password provided. SUBMISSION CONTENT: The typical submission material consists of: - A manuscript (A PDF file, with title, all authors with affiliations, abstract, keywords, highlights, inserted figures and tables and references), - Supplementary files: • a manuscript in a WORD file format • a cover letter (please see instructions for composing the cover letter) • a ZIP file containing figures in high resolution in one of the graphical formats (please see instructions for preparing the figure files) • possible appendicies (optional), cover materials, video materials, etc. Incomplete or improperly prepared submissions will be rejected with explanatory comments provided. In this case we will kindly ask the authors to carefully read the Information for Authors and to resubmit their manuscripts taking into consideration our comments. COVER LETTER INSTRUCTIONS: Please add a cover letter stating the following information about the submitted paper: 1. Paper title, list of authors and their affiliations. 2. Type of paper: original scientific paper (1.01), review scientific paper (1.02) or short scientific paper (1.03). 3. A declaration that neither the manuscript nor the essence of its content has been published in whole or in part previously and that it is not being considered for publication elsewhere. 4. State the value of the paper or its practical, theoretical and scientific implications. What is new in the paper with respect to the state-of-the-art in the published papers? Do not repeat the content of your abstract for this purpose. 5. We kindly ask you to suggest at least two reviewers for your paper and give us their names, their full affiliation and contact information, and their scientific research interest. The suggested reviewers should have at least two relevant references (with an impact factor) to the scientific field concerned; they should not be from the same country as the authors and should have no close connection with the authors. FORMAT OF THE MANUSCRIPT: The manuscript should be composed in accordance with the Article Template. The manuscript should be written in the following format: - A Title that adequately describes the content of the manuscript. - A list of Authors and their affiliations. - An Abstract that should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation, as well as the methodology employed. It should summarize the results and state the principal conclusions. - 4 to 6 significant key words should follow the abstract to aid indexing. - 4 to 6 highlights; a short collection of bullet points that convey the core findings and provide readers with a quick textual overview of the article. These four to six bullet points should describe the essence of the research (e.g. results or conclusions) and highlight what is distinctive about it. - An Introduction that should provide a review of recent literature and sufficient background information to allow the results of the article to be understood and evaluated. - A Methods section detailing the theoretical or experimental methods used. - An Experimental section that should provide details of the experimental set-up and the methods used to obtain the results. - A Results section that should clearly and concisely present the data, using figures and tables where appropriate. - A Discussion section that should describe the relationships and generalizations shown by the results and discuss the significance of the results, making comparisons with previously published work. (It may be appropriate to combine the Results and Discussion sections into a single section to improve clarity.) - A Conclusions section that should present one or more conclusions drawn from the results and subsequent discussion and should not duplicate the Abstract. - Acknowledgement (optional) of collaboration or preparation assistance may be included. Please note the source of funding for the research. - Nomenclature (optional). Papers with many symbols should have a nomenclature that defines all symbols with units, inserted above the references. If one is used, it must contain all the symbols used in the manuscript and the definitions should not be repeated in the text. In all cases, identify the symbols used if they are not widely recognized in the profession. Define acronyms in the text, not in the nomenclature. - References must be cited consecutively in the text using square brackets [1] and collected together in a reference list at the end of the manuscript. - Appendix(-icies) if any. SPECIAL NOTES Units: The SI system of units for nomenclature, symbols and abbreviations should be followed closely. Symbols for physical quantities in the text should be written in italics (e.g. v, T, n, etc.). Symbols for units that consist of letters should be in plain text (e.g. ms-1, K, min, mm, etc.). Please also see: http://physics.nist.gov/cuu/pdf/sp811.pdf . Abbreviations should be spelt out in full on first appearance followed by the abbreviation in parentheses, e.g. variable time geometry (VTG). The meaning of symbols and units belonging to symbols should be explained in each case or cited in a nomenclature section at the end of the manuscript before the References. Figures (figures, graphs, illustrations digital images, photographs) must be cited in consecutive numerical order in the text and referred to in both the text and the captions as Fig. 1, Fig. 2, etc. Figures should be prepared without borders and on white grounding and should be sent separately in their original formats. If a figure is composed of several parts, please mark each part with a), b), c), etc. and provide an explanation for each part in Figure caption. The caption should be self-explanatory. Letters and numbers should be readable (Arial or Times New Roman, min 6 pt with equal sizes and fonts in all figures). Graphics (submitted as supplementary files) may be exported in resolution good enough for printing (min. 300 dpi) in any common format, e.g. TIFF, BMP or JPG, PDF and should be named Fig1.jpg, Fig2.tif, etc. However, graphs and line drawings should be prepared as vector images, e.g. CDR, AI. Multi-curve graphs should have individual curves marked with a symbol or otherwise provide distinguishing differences using, for example, different thicknesses or dashing. Tables should carry separate titles and must be numbered in consecutive numerical order in the text and referred to in both the text and the captions as Table 1, Table 2, etc. In addition to the physical quantities, such as t (in italics), the units [s] (normal text) should be added in square brackets. Tables should not duplicate data found elsewhere in the manuscript. Tables should be prepared using a table editor and not inserted as a graphic. REFERENCES: A reference list must be included using the following information as a guide. Only cited text references are to be included. Each reference is to be referred to in the text by a number enclosed in a square bracket (i.e. [3] or [2] to [4] for more references; do not combine more than 3 references, explain each). No reference to the author is necessary. References must be numbered and ordered according to where they are first mentioned in the paper, not alphabetically. All references must be complete and accurate. Please add DOI code when available. Examples follow. Journal Papers: Surname 1, Initials, Surname 2, Initials (year). Title. Journal, volume, number, pages, DOI code. [1] Hackenschmidt, R., Alber-Laukant, B., Rieg, F. (2010). Simulating nonlinear materials under centrifugal forces by using intelligent cross-linked simulations. Strojniški vestnik - Journal of Mechanical Engineering, vol. 57, no. 7-8, p. 531-538, D0I:10.5545/sv-jme.2011.013. Journal titles should not be abbreviated. Note that journal title is set in italics. Books: Surname 1, Initials, Surname 2, Initials (year). Title. Publisher, place of publication. [2] Groover, M.P. (2007). Fundamentals of Modern Manufacturing. John Wiley & Sons, Hoboken. Note that the title of the book is italicized. Chapters in Books: Surname 1, Initials, Surname 2, Initials (year). Chapter title. Editor(s) of book, book title. Publisher, place of publication, pages. [3] Carbone, G., Ceccarelli, M. (2005). Legged robotic systems. Kordić, V., Lazinica, A., Merdan, M. (Eds.), Cutting Edge Robotics. Pro literatur Verlag, Mammendorf, p. 553576. Proceedings Papers: Surname 1, Initials, Surname 2, Initials (year). Paper title. Proceedings title, pages. [4] Štefanić, N., Martinčević-Mikić, S., Tošanović, N. (2009). Applied lean system in process industry. MOTSP Conference Proceedings, p. 422-427. Standards: Standard-Code (year). Title. Organisation. Place. [5] ISO/DIS 16000-6.2:2002. Indoor Air — Part 6: Determination of Volatile Organic Compounds in Indoor and Chamber Air by Active Sampling on TENAX TA Sorbent, Thermal Desorption and Gas Chromatography using MSD/FID. International Organization for Standardization. Geneva. WWW pages: Surname, Initials or Company name. Title, from http://address, date of access. [6] Rockwell Automation. Arena, from http://www.arenasimulation.com, accessed on 200909-07. EXTENDED ABSTRACT: When the paper is accepted for publishing, the authors will be requested to send an extended abstract (approx. one A4 page or 3500 to 4000 characters). The instruction for composing the extended abstract are published on-line: http://www.sv-jme.eu/information-for-authors/ . COPYRIGHT: Authors submitting a manuscript do so on the understanding that the work has not been published before, is not being considered for publication elsewhere and has been read and approved by all authors. The submission of the manuscript by the authors means that the authors automatically agree to transfer copyright to SV-JME when the manuscript is accepted for publication. All accepted manuscripts must be accompanied by a Copyright Transfer Agreement, which should be sent to the editor. The work should be original work by the authors and not be published elsewhere in any language without the written consent of the publisher. The proof will be sent to the author showing the final layout of the article. Proof correction must be minimal and executed quickly. Thus it is essential that manuscripts are accurate when submitted. Authors can track the status of their accepted articles on http://en.sv-jme.eu/. PUBLICATION FEE: Authors will be asked to pay a publication fee for each article prior to the article appearing in the journal. However, this fee only needs to be paid after the article has been accepted for publishing. The fee is 240.00 EUR (for articles with maximum of 6 pages), 300.00 EUR (for articles with maximum of 10 pages), plus 30.00 EUR for each additional page. The additional cost for a color page is 90.00 EUR. These fees do not include tax. Strojniški vestnik - Journal of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia, e-mail: info@sv-jme.eu http://www.sv-jme.eu Contents Papers Xiaoning Liu, Gengkai Hu: Elastic Metamaterials Making Use of Chirality: A Review Dezhi Wang, John E. Mottershead: Measurement Precision and Spatial Resolution with Kriging DIC Yang Liu, Sheikh Islam, Ekaterina Pavlovskaia, Marian Wiercigroch: Optimization of the Vibro-Impact Capsule System Ioannis P. Mitseas, Ioannis A. Kougioumtzoglou, Pol D. Spanos, Michael Beer: Nonlinear MDOF system Survival Probability Determination Subject to Evolutionary Stochastic Excitation Fabian M. Gruber, Daniel J. Rixen: Evaluation of Substructure Reduction Techniques with Fixed and Free Interfaces Lionel Manin, Stefan Braun, Damien Hugues: Modelling the Belt - Envelope Interactions during the Postal Mail Conveying by a Sorting Machine Martin Česnik, Janko SLavič, Miha BoLtežar: Assessment of the Fatigue Parameters from Random Vibration Testing: Application to a Rivet Joint 9770039248001