An equivalent approach for modelling butterfly-hysteresis passive variable friction damper

This research work proposes an equivalent approach modelling the peculiar hysteretic behavior of passive variable friction damper, for engineers who call for a simplified way to simulate such damper in anti-earthquake engineering practice with acceptable accuracy. This approach models the special nonlinear behavior of the damper via paralleling 3 basic link elements that commonly built-in commercially available structural analysis software, instead of developing direct hysteretic model that demands for both post development functionality on software side and programming skill on engineer side. The mathematical derivation and analytical geometric equating work are comprehensively conducted. The numerical accuracy of analysis output yielded using equivalent model is digitally verified against the one yielded using generic FEM simulation of the same damper subject to several loading cases. Code-automated optimization workflow for yielding product-wise modelling parameters is developed that is convenient and straightforward for application by engineers without requirement of high expertise, which is demonstrated by generating optimized parameters to equivalently model tested specimen. Passive variable friction damper can be modelled by assembling friction spring, elastic spring and wen model into parallel. Code-automated optimization process can locate the most ideal modelling parameters combination in predefined numeral fields. Hysteretic behaviors yielded from equivalent modelling approach match well corresponding ones from generic FEA simulation. Passive variable friction damper can be modelled by assembling friction spring, elastic spring and wen model into parallel. Code-automated optimization process can locate the most ideal modelling parameters combination in predefined numeral fields. Hysteretic behaviors yielded from equivalent modelling approach match well corresponding ones from generic FEA simulation.


Introduction
It is commonly recognized that the damage caused by earthquake induced ground motion to civil structure is adverse depending on its intensity, especially in those earthquake-prone region. Traditional approach addressing this impact is to either increase structure's strength or energy absorption capacity, while the latter one is usually achieved by energy dissipation contributed by plastic hysteresis at element level (well-known as 'plastic hinge zone'), which is demonstrated and observed as effective but with inevitable compensation as leaving localized structural damage [1].
According to the statistics, up to 2011 around 20,000 structures spreading in more than 30 countries have been equipped by either seismic isolation (SI) or supplemental energy dissipating (ED) techniques [2]. Drawing focus into ED system that this research relates to, which functions by adding additional damping capacity to the original structure, therefore reduces the demand of inelastic hysteretic energy dissipation by structural frame [3,4].
The prevalent and commercially available ED systems can be classified by controlling mechanisms as passively controlled, active controlled, and semi-active controlled [5]. Although with comparison to passive ED system, active ED system can be more event-adaptive, intelligent and can tackle wider band width of external excitation (mainly of earthquake, wind), its effectivity strongly relies on large external power supply which is probable to be inoperative during extreme earthquake event or energy-wasting if building is exposed to wind-frequent environment [6,7]. Semi-active controlled ED system has been researched by many scholars [8,9], which can be categorized as an intermediate portfolio between active and passive controlled ED systems. Differing to active ED, semi-active ED requires much less power supply in aspects of both life-time consumption and instantaneous power peak [10], which is only used for shifting specific mechanical parameters as an automated dynamic real-time optimization to improve efficiency of ED system [5,11]. Despite that high efficiency of energy dissipation plus considerable economy that semi-active ED solution can achieve [12,13], hindrances still exist, such as complexity in control system design and manufacturing [5], relatively low damping capacity and limited reliability due to the complexity of its constituent mechanics [10,14]. The last part forming the ground of ED systems that is right the one dominantly adopted in earthquake engineering by this far across worldwide, belongs to passive ED systems. Regarding energy-dissipating mechanisms, they can be categorized as displacement-dependent (metallic yield damper, friction device, buckle-restrained brace, etc.), velocity-dependent (viscous fluid damper, etc.) and mass tuning (tuned mass damper, tuned liquid damper, etc.) and hybrid type (viscoelastic damper, etc.) [1,15].
This research focuses on a newly developed branch of friction ED systems, also called friction damper (FD), which follows or quasi-follows Coulomb's friction law. By this far many prototypes of FD have been researched and developed [3]. Owing to its good economy, mechanical simplicity, easy applicability and large hysteretic energy dissipation capacity, FD has gained very wide practical application hence big market. However, without a high initial static force that builds up sufficiently large hysteretic loop at which FD starts sliding to dissipate energy, there will not be such large damping capacity it can supply to structure as an impactful advantage of FD, whereas large initial static force does stiffen global structure and consequently induces large floor acceleration [16]. Namely, under frequent dynamic loading case with low or moderate intensity, it rarely plays the role of dissipating energy, but inversely induces larger floor acceleration instead [17].
To overcome the shortcoming of FD outlined above, a new branch of FD portfolio has been developed in various patterns, called passive variable friction damper (PVFD). It is worth of clarification herein that, there is no configuration of response sensor nor real-time optimizing computing chip configured to vary the factors (mostly to shift normal force) that determines the friction force, which shall be categorized as semi-active damper. Instead, PVFD employs pure mechanical adaptivity to adjust friction force via varying normal force, as a response to either displacement or rotation [18,19]. The categorization of PVFD can hence be located in the middle of passive and semi-active ED system, which employs pure mechanical adaptivity to vary damping force whereas less intelligent comparing to semi-active ED system. Given below ( Fig. 1) are 3 typical patterns of hysteretic loop of PVFD devices, according to illustrations originally delivered by Wang et al. [18], Zhou and Peng [19], CSI [20].
A considerable quantity of research has been conducted on several prototypes of PFVD, most of which yields similar remarks on PVFD that, with comparison to traditional FD, PVFD manifests better improved control performance on structural lateral drift. However, application of PVFD in real engineering practice remains limited due to that very few commercially available structural FEA software provides PVFD link elements for modelling use. Well-known ones SAP2000 and ETABS, produced by CSI, America, provides a link element named 'friction spring' , which enables the simulatability of the hysteretic behavior of Dual-flag PVFD [ Fig. 1c], but lack of the former two. Barzegaret al. [17], Downeyet al. [21], Downeyet al. [22] employed the LuGre friction model to simulate the dynamic behavior of PVFD, which accounts the stick-slip motion along with Stribeck effect. The effectivity and simplicity of LuGre model for modelling behavior of PVFD is evidenced in research [22], which is also highly numerically adjustable to fit wider parametric range of PVFD products [17]. However, there are 3 reasons that the efficiency and practicability of LuGre model is not widely adopted in civil engineering practice. First, by this far, most commercially available structural FEA software build in no LuGre model, which nevertheless is possible to be enabled via the post-development interface but requires relatively higher expertise and provided post-development function by adopted software. Second, it requires additional experiment for determining the modelling parameter, which compromises by cost and time efficiency. Third, because LuGre model involves velocity and is evidenced as rate-dependent [23], thus requires an optimization process for yielding acceptable modelling parameters [21]. Even LuGre model is experimentally demonstrated as capable to simulate more explicit friction behavior considering memory-dependent mechanism beyond Coulomb's friction model [23], but the need of reaching such depth of analytical accuracy is seldom called in civil engineering practice unlike aerospace ones, namely, it is an expensive compromise to induce a series of unnecessary impacts such as rate-dependency, multiple experimental parameter-determining cycles and furtherly repetitive parametric optimization process, only just for simulating the dynamic behavior of PVFD. To view wider, Barzegar et al. [17] conducted numerical modelling and analysis via a way to highly simplify a benchmark building into a lumped-mass shear stand, so that enabled the application of LuGre model achieved in pure coding environment. Zhou and Peng [19] have made use of ANSYS (generic FEA software) to simulate the explicit physical behavior of PVFD, yielded hysteresis matches the one from a benchmark building equipped by PVFD but modelled using program IDARC that is developed in University of Buffalo, which is mainly for academic use.
This paper introduces an equivalent approach of modelling the peculiar hysteretic behavior of the butterfly-hysteresis PVFD [ Fig. 1a], which can be achieved in commercially available structural analysis software with practical ease, if which builds in required basic link elements. One of those capable ones, SAP2000 by CSI, America [20], is employed in this research. In the following section, the physical mechanism resulting butterfly-pattern nonlinear behavior is illustrated and mathematically derived by means of vector superposition. Section 3 introduces the formation of paralleled link assemblage that can equivalently model the nonlinear behavior of PVFD by means of pure analytical geometric equating. In Sect. 4, an optimization workflow along with coded iterative algorithm is developed to automate the numerical process determining the most optimized modelling parameters in predefined variable searching fields according to input control points extracted from experimental hysteresis of targeted PVFD product. To demonstrate the effectiveness of the parametric optimization workflow and algorithm, such a work towards tested specimen is conducted as an exemplification, following which, a comprehensive generic FEA model is configured in Abaqus to verify the numerical accuracy of using the equivalent model.

Mechanical characteristics
One PVFD (Fig. 2) product is developed and researched prior to this study that resembles general plate FD which is assembled by an intermediate friction panel sandwiched by 2 clamp boards, all 3 pieces are fastened together using screws and being pre-compressed to exert normal force  Step history of displacement loading for generating friction. Resultant damping force is being increased while the damper stretches outwards along the end upheaving ramp, which consistently tightens spring fasteners hence increases total clamping force towards maximum deformation. Schematic of 3D model of tested specimen of PVFD are presented in Fig. 3, along with experiment configuration (Fig. 3). External load is exerted periodically in displacement pattern (Fig. 4). Hysteresis of developed and tested specimen is in pattern of butterfly-hysteresis one (Fig. 1a). Yielded hysteretic loop is smoothed into piecewise form and fitted regarding theoretical derivations (Fig. 5) given in Eqs. 1-5. Reaction force represented by each colored piece on loop is derived based on Coulomb's friction theory when in sliding, and linear elastic assumption while deforming without sliding.
In Eqs. 1 ~ 5, F 1 -F 5 denote reaction forces represented by colored segments on hysteretic loop ( Fig. 5), respectively. Grayed segments denoting the central symmetric half can be obtained via rotating the colored pieces altogether around origin by 180 degrees in F-X plane. Symbols x 1 -x 4 denote the deformations of 4 joint points at which the reaction force changes. Symbols k i , f c , µ, k s , θ represent 5 constants of, respectively, stiffness of elastic-slip, total precompression force, friction coefficient of slipping faces, aggregated stiffness by all spring fasteners, angle of end sloped faces measured to the flattened middle face. Regarding Fig. 2, symbols f x1 , f x2 , d s denote 3 intermediate variables, which are respectively, projection of friction force to the spatial plane coplanar with the flattened middle face where initial friction F 2 acts on, projection of normal compression force to the same spatial plane, compressive deformation of spring fasteners. Particularly, regarding Fig. 5 and referring to Eqs. 3 and 5, F 3 has relatively larger magnitude than F 5 so that the slope of orange segment is steeper than that of purple, because f x2 plays a resistive role during offcentering motion, whereas plays a motivating role in centering motion.

Equivalent modelling approach
By this far, there is very few built-in link elements available for modelling the butterfly-hysteresis PFVD in structural FEA software, especially of those prevalent and commercially available ones, while such innovative and evidenced as efficient ED device calls for more promotion and wider application in real engineering practice. In this research, an equivalent modelling approach has been proposed, via assembling 3 built-in link elements that are dual-flag friction spring (DFFS), negative-stiffness multi-linear elastic (NSMLE), and plastic Wen (PW) dampers, into parallel as a 'super-link' (Fig. 6). 4 shoulders in stiffening stage of the super-link are not as one single slope but linked two, which can be reasoned from the perspective of geometry to detailed view where the hysteresis of DFFS begins to be cut through by NSMLE. If the slopes of these 4 shoulders in fitted hysteresis (Fig. 5) can well match the secant slopes of the numerically assembled one (Fig. 6), such distortion in shape can be deemed as acceptable.
The super-link equivalently modelling butterfly-hysteresis PFVD, has been numerically tested in SAP2000 environment using single-degree-of-freedom (SDoF) subjected to harmonic displacement loading with amplitude valued by maximum deformation capacity of tested PVFD device, yielded hysteresis of which is given in Fig. 7.

Decomposition of hysteretic loop
Unlike those built-in link elements, parameters of assembled super-link used to equivalently model the peculiar hysteretic pattern cannot be determined from the real mechanical property of PVFD device, neither from experiment nor standardized database of product portfolio. On the other hand, because of the geometric shift at 4 shoulders on the hysteresis of super-link comparing to the experimental one, there is not exact analytical solution of modelling parameters from given real mechanical parameters of PVFD that determines its unique hysteretic loop. Therefore, an optimization methodology better with automated workflow for practical ease, is called to be invented for conducting modelling parameters to fit various PVFD device as representable as possible. Figure 8a-c symbolized parameters determining complete nonlinear behavior in SAP2000 environment for link components DFFS, NSMLS and PW, respectively, which are all pending to be yielded as final output from optimization work. According to the complete coupling of deformation of 3 components assembled in parallel, analytical piecewise F-X function of the super-link can be derived via making arithmetic summation of 3 components that geometrically illustrated in Fig. 6.
The decomposition diagram given in (Fig. 9) directs the focus from the total hysteresis of super-link into the isolated shoulder region, because the contribution of PW is relatively the most simplistic one out of all 3 link components hence the most challenging geometric complexity of this synthetic hysteresis drops into the 4-shoulder region, which is completely dependent on contribution of DFFS and NSMLE components. Nonetheless PW's contribution will be accounted for optimization in terms of energy dissipation within one loading circle and will be addressed in following section.

Numerical optimization approach
The isolated region incorporating 2 shoulders in 1st and 4th quadrants that red-outlined in Fig. 9, has been symbolized into Fig. 10, symbols assigned to coordinates of inflection points are collected in Table 1, each of which can be mathematically determined based on input 7 parameters that extracted from the real hysteretic loop of targeted PVFD device. Derivation of tabularized coordinates in Table 1 are given by Eqs. 6-17.
Research Article SN Applied Sciences (2022) 4:234 | https://doi.org/10.1007/s42452-022-05098-6 In Eqs. 6-17, d 0 , d u , respectively denote the positive displacement value at which PVFD initiates its friction increasing mechanism during off-centering motion, and absolute value of its maximum deformation capacity, hence identical to x 2 , x 3 in Fig. 5. K f0 , k n2 , k f1 , k f2 are variables pending to be optimized into one optimal combination within predefined searching range (called optimization variables), k n1 is forced to equate k f0 for purpose of cancelling each other in inner constant friction platform. Shifts of hysteretic geometry of isolated part if arbitrarily varying these 4 optimization variables are exemplified in Fig. 11, each set is obtained via varying only one variable while keeping the rest 3 as constant. 5 indices quantifying fitting quality (Fig. 12) of optimization output, are set as yielding force (F y ), maximum positive reaction forces (F max ), minimum negative force (F min ), displacement where F min is reached (x 4 ), and maximum energy dissipation in one hysteretic cycle (E max ). F y , F max , and E max are commonly recognized as significant mechanical parameters of displacementdependent ED device in terms of anti-earthquake performance, while F min and x 4 make significant impact the geometric fidelity in 1 st and 4 th quadrants. Regarding illustration delivered by Fig. 12, F y is identical to F py (Fig. 8(c)), F max can be obtained by summation of F c (Table 1) and F py ,. F min is represented by F d (Table 1), x 4 is represented by X d (Table 1). In following content, superscript dot notated over these 5 indices is to distinguish the value of super-link with corresponding one from experimental PVFD hysteretic loop.
The area of hysteresis, denoted as E max , reflecting energy that the super-link can dissipate in one hysteretic cycle, can be calculated via making summation of its 3 decomposed regions given in Fig. 13. Analytical solution of E' max is derived and presented by Eq. 18. Particularly, real values F y , F max , F min , x 4 can be directly read from hysteretic loop of tested PVFD, E max should instead be calculated via integrating '(F n+1 + F n )·dx/2' over recorded data points series within one complete loading cycle.  Figure 14 presents the interactive workflow for engineers to generate optimized modelling parameters via simply inputting the mechanical parameters regarding targeted PVFD device that is to be applied into engineering practice.

Summary of optimization workflow
In particular, the second step requires the searched ranges of 4 variables to be defined, according to either analytical derivations employed for optimization, or visual patterns of hysteresis of equivalent super-link to experimental one, k f1 and k f2 , of DFFS component, must be less than largest secant value from 4 shoulder regions of experimental hysteresis, k n2 must be less than k f0 . K fo on pure analytical ground, can be arbitrary adopted up to infinite, because it will be cancelled out by k n1 (identical to k f0 ) of NSMLS component. However, larger searching range defined for k f0 will not only consume more computational efforts, but also cause convergence difficulty if being modelled in SAP2000. Based on above considerations, ground of all 4 searching ranges can be set as 0, upper bound of searching range of k f0 , is suggested to be set no more than twice of k e from experimental hysteresis; upper bound of k n2 can be set identical to k f0 ,; upper bound of k f1 can be set as secant value of 1 st quadrant shoulder region from experimental hysteresis, and lower bound of k f2 can be set as negative of value set to k f1 .

Optimization for modelling parameters of tested specimen
To demonstrate an application of using proposed optimization method for fast determining the optimal parameters of components assembled to equivalently simulate nonlinear behavior of butterfly-hysteresis PVFD, numerical modelling and analysis of such super-link connecting a SDoF, within SAP2000 environment, has been conducted. Property definition of 3 link components assembling into the super-link, are for instance displayed in Fig. 15. Displacement load imposed to SDoF model is set in harmonic pattern with amplitude of 0.03 m identical to the maximum displacement capacity of tested PVFD. Nonlinear direct integration method is employed for computation with no additional damping defined, time step for which is set as sufficiently fine for catching desired nonlinear behavior of super-link. Parameters from hysteresis of tested PVFD specimen are referred as optimization input (as constants) ( Table 2). For illustration purpose, in comparison to experimental one, 16 hysteresis graphically presented in "Appendix A", are yielded from SDoF assigned super-link in various parameter combinations, each model in collection differs from the others by shifted error ratio predefinitions hence completely different output modelling parameters. Table 3 gives the error ratio predefinition and output optimization variables for all trials that correspond to yielded hysteresis given in "Appendix A" by same order. Particularly, error ratio 'Err_Fmax' is employed to judge how well F`m ax can match F max , which is not predefined but output by each turn of running code instead. Because of the geometric difference between the real hysteresis of PVFD and equivalently modelled one, there will yield no solution if input 4 error ratio predefinitions (Err_F y , F min , x m , E max ) and 1 output judging ratio (Err_F max ) were all to be 0. In another words, the optimization process involves personal engineering judgement, any one larger of 4 input error ratio predefinitions will leads to smaller output of F`m ax hence larger absolute value of Err_F max , and vice versa. Therefore, engineering judgement should be employed to determine those relatively more interested variables according to the physical meanings they represent, then to diminish allowable error ratios for interested those, error ratios of the rest can hence be allowed to be adjusted incrementally until yielding a solution.
Subject to judgement from engineering perspective in this case, indices E max , F y , and the post-evaluated F max , are key mechanical parameters as for added-structuredamping use, in terms of its energy dissipation capacity and damping force magnitude. X 4 is relatively less significant but will impact small-displacement behavior if deviate too much from experimental value. Comparing to mentioned 3 indices, F min is the least impactful one out of all 5 ones. Therefore, the optimized parameter set for modelling super-link equivalent to tested PVFD product, is determined to be the one having relatively smaller error ratios in F y , x 4 , E max , while adjusting error ratio of F min to run code until yielded F`m ax satisfies to have small enough error ratio Err_F max . Error ratio predefinitions finally justified as optimal for tested PVFD, are given in Table 4, while output modelling parameters are given in Table 5. Hysteretic loop yielded by the super-link modelled in SAP2000 configured with output parameters (Table 5), is presented in Fig. 16 against experimental one.

Configuration of solid FEA model
The equivalent modelling approach for butterfly-hysteresis PVFD is digitally tested under SAP2000 environment, yielded F-X relationship meets expectation and sits strictly within error-ratio-confined numeral domain. However, whether the numerical accuracy of modelling PFVD from proposed super-link assembly drops in commonly acceptable level while conducting practical structural simulation using structural FEA (in this research, SAP2000), remains uncertain. To verify whether the equivalent modelling approach achieved in SAP2000 can give desired accuracy, another digital verification work has been conducted via setting a comprehensive solid FEA model simulating tested PVFD, utilizing Abaqus CAE software. As shown in Fig. 17, according to the geometry and physical property of tested PVFD product, a solid FEA model has been built in Abaqus CAE. Q355 steel is assigned to clamp board and middle slider while material property of a patented highly compacted layered polymer is assigned to friction panel part. Tie constraint is assigned onto the interface of steel clamp board to friction panel. Because the nonlinearity of PVFD is purely induced by friction mechanism, materials of both steel and polymer are defined with isotropic elasticity only. According to product specification, Young's modules are set as 2 × 10 11 Pa and 4 × 10 10 Pa for steel and polymer material, respectively. Poisson ratio for steel is set as 0.3 and 0.28 for polymer material. Mass density for both is assigned unity because the analysis is to be quasi-static and not accounting for gravity, so that mass induced inertia force and gravity can be excluded. Analytical section for all parts is defined as solid homogenous. In terms of the two-sided friction face, hard contact is defined for its normal behavior, penalty friction formulation is employed to catch its tangential sliding behavior with friction coefficient 0.62 and absolute elastic slipping distance 1 mm defined, while the elastic slipping distance cannot be determined directly from neither test output nor product specification hence is adjusted to match the tested non-slipping force-deformation relationship. According to specimen test results, the slipping behavior does not manifest rate-dependent manner within tested cycles, so the friction force is not assigned rate-dependency definition. To accurately catch the process of stiffness varying mechanism, structured mesh is employed and size of which is capped at 0.007 m (less than 0.01 m where the modelled PVFD varies stiffness). Screw fastener is equivalently modelled as linear spring connector with only axial degree of freedom enabled and constrained each end to the outermost surface edge of circular hole that preserved for fastener. Stiffness value for each connector section is defined identical to the one of precompressed stacked-spring disk.
Load is assigned in multi-stepped displacement form on to a reference point that being constrained to lead the motion of middle slider, in correspondence to SDoF simulated in SAP2000. Fix-end boundary condition is then assigned to clamp board's end on the side opposite to loading with U3 degree of freedom released. Pre-compression force acting to fastener is defined as connector force exerted in multi-stepped ramp pattern prior to and propagated in displacement loading process.
Like other categories of structural ED system, PVFD product aims to dissipate earthquake energy input to structure, so the uncertainty of ground motion should be considered while assessing performance of ED system. Inspired by approach for such digital testing 24 , to account such uncertainty, ground motions of 4 well-known big earthquake events occurring in 3 different regions, have been selected from PEER ground motion record database [24], as direct displacement loading function (integrated twice from original acceleration into displacement form) to compel PVFD modelled in both SAP2000 and Abaqus, instead of using simplistic harmonic function. Figure 18 displays these 4 waves, maximum amplitudes of which are all normalized and scaled to match the maximum displacement capacity of tested PVFD.

Result and discussion
Response time history extracted from analysis result is given (Fig. 19) in comparative pattern of SAP2000's SDoF against Abaqus's solid FEA.
The comparisons above manifest very close match in each pair. Therefore, the feasibility and accuracy of applying such equivalent modelling approach to simulate the peculiar nonlinear behavior of butterfly-hysteresis PVFD in engineering practice, is hereby validated in digital way.

Conclusions
PVFD has been developed by many researchers so far, which is demonstrated as possessing improved anti-earthquake performance comparing to traditional FD. To simulate its peculiar nonlinear behavior including negative stiffness, several methods have been proposed by researchers, the numerical accuracy and applicability however vary, nonetheless the threshold of applying of PVFD in real engineering practice remains higher than the other general link elements that are widely built in popular structural modelling and analysis software. In this research, an innovative equivalent modelling approach addressing such difficulty has been developed comprehensively. In advance, a product of PVFD has been developed, manufactured, and tested, which features butterfly-pattern hysteretic behavior.
An equivalent approach for modelling the peculiar hysteretic behavior of butterfly-hysteresis PVFD that assembles 3 general link elements in parallel into a super-link element, has been proposed and numerically tested in SAP2000. Hysteresis yielded from this super-link under displacement loading gives close match to the tested one.
An optimization method has been developed, via which to yield modelling parameters best fitting targeted butterfly-hysteresis PVFD product within predefined searching ranges and error ratios. To enable most engineers to apply this method in engineering practice with ease, a coded algorithm has been written, which can automate entire workflow producing product-wise modelling parameters for defining 'superlink' element, while being optimized to meet customizable criteria that judges fitting rationality and fidelity. An example of applying such workflow for tested butterfly-hysteresis PFVD specimen has been conducted, yielded super-link element meets predefined optimization criteria as expected. Only limited human intervention is required throughout whole process.
Numerical accuracy of using assembled super-link element has been verified by means of digital comparative work. In the work, generic 3D solid FEA model has been built in Abaqus that in real physical detail simulates tested PVFD device. According to comparison made to hysteresis yielded by super-link in SAP2000 and Abaqus model, under several displacement loadings accounting for random feature of natural ground motions, a very close match is achieved of using these 2 modelling approaches. Therefore, the numerical accuracy of modelling hysteretic behavior of butterfly-hysteresis PVFD using equivalent super-link has been verified.
Although the equivalent approach modelling the peculiar nonlinear behavior of PVFD is straightforward for both understanding and practicing, with acceptable error rate in terms of key parameters and digitally verified good accuracy, the work of validation and verification still need to be undertaken in other ways. In future research, it is recommended to use comprehensive experimental tests to sufficiently verify the accuracy and effectiveness of proposed equivalent modelling approach. As for more detailed comparison and validation purpose, this approach is suggested to be digitally tested within different structural analysis software with 3 basic link elements built-in, or conducting same work by means of investigating the direct formula of the hysteretic model then programing which into analysis software that supports post development functionality in terms of defining new hysteretic rule.
Funding The research leading to these results received funding from 'research for improving performance of disaster resistance of traditional bamboo-structured civil houses-a subsidiary work to Chinese nationally prioritized research plan' under Grant Agreement No 2020YFD1100703-04.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors have no relevant financial or nonfinancial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.