Nonlinear vibration of knitted spacer fabric under harmonic excitation

Knitted spacer fabrics can be an alternative material to typical rubber sponges and polyurethane foams for the protection of the human body from vibration exposure, such as automotive seat cushions and anti-vibration gloves. To provide a theoretical basis for the understanding of the nonlinear vibration behavior of the mass-spacer fabric system under harmonic excitation, experimental, analytical and numerical methods are used. Different from a linear mass-spring-damper vibration model, this study builds a phenomenological model with the asymmetric elastic force and the fractional derivative damping force to describe the periodic solution of the mass-spacer fabric system under harmonic excitation. Mathematical expression of the harmonic amplitude versus frequency response curve (FRC) is obtained using the harmonic balance method (HBM) to solve the equation of motion of the system. Parameter values in the model are estimated by performing curve fit between the modeled FRC and the experimental data of acceleration transmissibility. Theoretical analysis concerning the influence of varying excitation level on the FRCs is carried out, showing that nonlinear softening resonance turns into nonlinear hardening resonance with the increase of excitation level, due to the quadratic stiffness term and the cubic stiffness term in the model, respectively. The quadratic stiffness term also results in biased vibration response and causes an even order harmonic distortion. Besides, the increase of excitation level also results in elevated peak transmissibility at resonance.


Introduction
The human body is sensitive to vibrations of various frequency ranges. For example, internal human organs and the nervous system are negatively affected by the frequency ranges from 4 to 8 Hz, and from 20 Hz to below 250 Hz, respectively. 1 However, occupations dealing with electrical and pneumatic powered rotary tools and processes in farming, mining, quarrying, demolition, and road construction inevitably produce dangerous vibrations in specific frequency ranges. 2 Prolonged exposure to vibrational environments may cause discomfort and even mechanical injuries and other diseases to workers in these professions. For example, workers operating hand-held power tools may develop hand arm vibration syndrome (HAVS) such as vibration-induced white finger (VWF) and peripheral neurological disorders. 2 To minimize the risks of vibration injury, vibration dampening materials can be used to attenuate the magnitude of vibration. However, traditional materials such as rubber sponges and polyurethane foams may have problems in comfort to wear and recycling. As a potential substitute, knitted spacer fabric is breathable, reusable, versatile and has good compression property. Especially for weft-knitted spacer fabric, it exhibits good conformability due to the loop nature of weft-knitted stitches. In practical applications, the cushioning performance of knitted spacer fabric as functional seat cushions or anti-vibration gloves is tool or excitation spectrum specific, because different vibrating machines induce different frequency spectrum and intensity of the vibration. 3 It is also common knowledge that in a linear mass-spring-damper vibration system, a low stiffness results in a broadened range of excitation frequencies in which the vibration can be isolated. Besides, a high damping is beneficial when vibration at resonance is a concern. Nevertheless, reducing the stiffness of the cushioning material may become disadvantageous when the excitation frequencies are very low; similarly, increasing damping may also be disadvantageous when the excitation frequencies are much higher than the resonance frequency. As a result, to minimize adverse vibrations and slow down the onset of vibration induced diseases, fabric design should match the specific application scenario. 4 Knitted spacer fabrics are sandwiched textiles which usually consist of two outer layers connected but kept apart by a spacer layer. They have been made into compression bandage, 5,6 wound dressing [7][8][9] and impact protector 10-12 due to their particular mechanical and thermophysiological properties. For these uses, the compression behavior plays a very important role. To date, the compression behavior of knitted spacer fabrics has been extensively studied using experimental, numerical and analytical approaches. [13][14][15][16][17][18][19][20][21][22][23] As a type of three-dimensional textile structure made of polymeric materials, knitted spacer fabric exhibits nonlinear compression force. Liu and Hu 16 experimentally identified that the compression force versus displacement relationship of knitted spacer fabric could be divided into the initial stage, the linearly elastic stage, the plateau stage and the densification stage. They 19 also developed a finite element model to strengthen the understanding of the compression mechanism of spacer fabric structure, finding that the nonlinear compression behavior is due to the post-buckling, torsion, shear, rotation, contacts between spacer monofilament, and contacts between spacer monofilament and outer layers. Parametric study also showed that smaller monofilament inclination, coarser monofilament and smaller fabric thickness result in higher compression resistance. Liu and Hu 17 also proposed a constitutive model consisting of seven parameters to describe the compressive stress-strain relationship of knitted spacer fabric. The proposed model outperforms three existing constitutive models by having the smallest fitting errors.
The quasi-static compression behavior is closely related with the pressure relief 24 and the impact protection 10,25 performance of knitted spacer fabric as a cushioning material. Relevant research has already shown its energy absorption capacity from the perspectives of quasi-static compression and impulsive loading. Studies on the vibration properties and vibration dampening performance of knitted spacer fabrics have also been carried out in recently years. Blaga et al. [26][27][28][29] studied the impact response of knitted spacer fabrics by employing the free vibration method. It was found that knitted spacer fabric should have both the capacity of absorbing vibration energy and sufficient stiffness to avoid its collapse. 28 Frydrysiak and Pawliczak 1 made a comparative study on the vibro-insulation properties of knitted spacer fabric and typical upholstery materials as vibration damping inserts in workplace seating, and concluded that knitted spacer fabrics can be a viable alternative to typical sponges, with the unquestionable advantage of having much lower thickness. Krumm et al. 30 compared the vertical seat transmissibility of warp knitted spacer fabrics and standard foam cushion, and found that the design of the seat backrest cushion should prefer warp knitted spacer fabrics, while the design of the seat pan cushion should prefer foams. Liu and Hu 31 experimentally investigated the vibration isolation properties of warp-knitted spacer fabrics under harmonic excitation, and found that thicker fabric possessed lower resonant and isolation frequencies, and thus showed better vibration isolation performance. Chen et al. 32,33 studied the vibration transmission property of warp-knitted spacer fabrics under harmonic excitation and under damped free vibration conditions, respectively, based on a linear mass-spring-damper model. They found that the structural and material properties of warp-knitted spacer fabrics have significant effects on vibration transmission properties. However, taking knitted spacer fabric with top-loaded mass as a linear vibration system may be inadequate. Chen et al. 2 has experimentally studied the vibration isolation performance of weft-knitted spacer fabrics under harmonic excitation, and found that the curve shape in the transmissibility curve is asymmetric which bends to left at resonance, indicating the mass-spacer fabric vibration system is nonlinear. The study also showed that changing the load mass and the vibration intensity changes the loading conditions of the fabric structure, and thus also changes fabric stiffness and vibration isolation performance.
Based on the current research background, this study seeks to find a mathematical model to describe the vibration behavior of knitted spacer fabric under harmonic excitation. As a linear model only applies to very small vibration, it is necessary to take nonlinearities into account in building the equation of motion when the vibration is large. Nevertheless, if adopting the complex nonlinear elastic force-displacement relationships suitable for describing the quasi-static compression behavior of polymeric materials, such as the abovementioned Liu's seven-parameter model 17 and the LS model, 34 which divides the typical compression stress-strain response of polymer foam into three regions, the vibration model will be difficult to arrive at an analytical solution. To circumvent this issue, force nonlinearity in the vibration model can be expressed in the form of a cubic polynomial consisting of linear-quadratic-cubic stiffness terms, instead of complex nonlinear force expressions. The quadratic stiffness term indicates force asymmetry.
Moreover, due to the viscoelasticity of polymer materials, a fractional derivative term is adopted in the vibration model to describe the time-dependent property of spacer fabric. The fractional-order derivative was raised by Leibniz more than 300 years ago. 35 Its nonlocal property indicates that the history states of a system have an influence on its current state. Moreover, it has a fading memory in that it weighs the recent past more heavily than the distant past. 36,37 This property of having an unlimited and decaying memory makes fractional differential equations applicable to a description of complicated dynamical systems in the real world. 38,39 Fractional calculus has become a popular instrument in many scientific and engineering fields such as viscoelasticity, hereditary physics, structural hysteresis, rheology, electrochemistry, bioengineering, mechanics, automatic control, signal and image processing, quantum evolution of complex system, etc. 35,40 Bagley and Torvik 41 has used fractional derivate models to describe the frequency-dependent damping behavior of materials and systems very well. Deng et al. 36 and Deng 42 have also demonstrated the fractional derivative model is applicable in describing the viscoelasticity of flexible polyurethane foam during vibration.

Fabric sample preparation
Weft-knitted spacer fabric was fabricated on a 14-gauge STOLL CMS 822 computerized flat knitting machine. Front and back needle systems produced two separate outer-layer pieces, which were knitted with single jersey structure using one 100D nylon multifilament yarn and one 30D Spandex/70D nylon covering yarn together. Simultaneously, polyester monofilaments with 0.12 mm diameter made tuck stitches every six-needle distance on two needle systems in turn. In this way, spacer monofilaments connected two outer-layer fabrics, and thus a sandwiched structure was obtained. After steaming treatment, fabric thickness was increased as the shrinkage of elastic outer layers caused the rotation of spacer monofilaments. Fabric samples were further relaxed under an environmental condition of 20°C and 65% relative humidity for 5 h. The diagram of fabric structure and structural characteristics of fabric including surface loop density, areal mass and fabric thickness are listed in Table 1. It is worth noting that only half numbers of needles participated in making the tucking texture, in which the tuck stitches were evenly distributed over the whole outer layers. Assuming using all needles to knit tuck stitches, the amounts of monofilaments consumed will be doubled, limiting the shrinkage of outer layers during the steaming treatment. Besides, monofilaments inside the spacer will become compacted and collapsed rather than standing up, resulting in much lower fabric thickness and less supporting capacity. Therefore, the structure with all needles knitting tuck stitches was not preferred.
The cross-sections of the fabric along the course direction and the wale direction are shown in Figure 1(a) and (b), respectively. Spacer monofilaments have a cross-over structure along the course direction, and a curved shape along the wale direction. As shown in Figure 1(b), as the linking points A and B of each monofilament on two outer layers are not aligned vertically, relative slip between two outer layers can happen along the wale direction when a vibration force or a compression force is applied. To avoid the transverse instability, two identical spacer fabrics were bonded together using a double-sided adhesive tape as shown in Figure 1(c). In this way, the topmost layer and the base layer of the laminated fabric can maintain right opposite to each other when bearing load. Figure 1(d) shows the quasi-static compression curve of the laminated spacer fabric. Clearly, force nonlinearity takes place, as is also indicated by the stiffness curve.

Vibration test
An electromagnetic vibration shaker equipped with a 350 mm × 350 mm platform was used to measure the acceleration transmissibility T of the laminated weftknitted spacer fabric. The acceleration transmissibility T is defined as the ratio of the absolute value of the acceleration of the mass to the absolute value of the acceleration of the base platform. T > 1 means vibrations is amplified, while T < 1 means vibration is reduced. Spacer fabric of a flat area of 150 mm × 150 mm was top-loaded with a 2 kg metallic block of a flat area of 90 mm × 90 mm, and they were placed on the center of the shaker platform. The shaker was excited by sinusoidal sweeps from 4 Hz to 500 Hz with a sweep rate of 1.0 Oct/min. Each sweep event had a constant excitation level G , namely, 0.1 g, 0.2 g, or 0.3 g, where g is the gravitational acceleration (9.81 m/s 2 ). To maintain a constant excitation level, the excitation amplitude will decrease as the driving frequency increases. Acceleration transmissibility curve at a frequency range was obtained during one sweep event.
Three replicates were carried out for each testing condition. In addition to the acceleration transmissibility T , the harmonic amplitude 2 A was also calculated from T and the corresponding phase angle ϕ using the equation , which is based on the reference work, 43 where ω is the angular frequency. The frequency response curves (FRCs) of the fabric from the experiment, that is T and 2 A in response to the excitation frequency, are shown in Figure 2. The curve shape near resonance peak is close to symmetric for 0.1 g excitation level, in which case the dynamic load and dynamic deformation for spacer fabric is relatively small, so it can be treated as a linear vibration system. For 0.2 g and 0.3 g excitation levels, however, the vibration behavior of the mass-spacer fabric system exhibits softening nonlinearity as the curves bend to the left side around resonance peaks. It is also shown that as the excitation level increases, the harmonic amplitude increases as well, but at the same time, the degree of nonlinear softening increases, that is, the peak acceleration transmissibility and the resonance frequency decreases, which implies the decrease of dynamic stiffness of the mass-spacer fabric system, referring to the equation f k m r = 1 2π . In brief, due to the nonlinearity of the mass-spacer fabric system, it is necessary to use a nonlinear model to describe its vibration behavior.

Building the equation of motion
Since the vibration takes place in the vertical direction, the mass-spacer fabric vibration system is simplified as unidimensional. Due to polymeric viscoelasticity, the dynamical mechanical behavior of spacer fabric is different from the static mechanical behavior in that the history states also affect the current state of the dynamical system. As a result, the stiffness coefficients identified from the quasi-static compression curve in Figure 1(d) cannot well describe the vibration behavior of spacer fabric. Therefore, a phenomenological model is proposed in this study based on the nonlinear softening phenomenon observed in the vibration experiment. In this case, spacer fabric is treated as a nonlinear spring with linear, quadratic and cubic stiffness coefficients k , k 2 , and k 3 , plus a linear viscous damper with damping coefficient c , as shown in Figure 3 When the quadratic stiffness coefficient k 2 is zero, the nonlinear spring consisting of linear and cubic stiffness terms is called a Duffing isolator. The Duffing equation was first introduced by the German electrical engineer Georg Duffing 44 in 1918, and was widely used for describing many physical, engineering and even biological problems. 45 In order to describe the asymmetric elastic force-deflection relationship for spacer fabric, indicated by the nonlinear softening phenomenon observed in the vibration experiment as shown in Figure 2, an additional quadratic stiffness term k x 2 2 is introduced, resulting in the Helmholtz-Duffing equation, which has been used to describe the vibration behavior in circumstances such as (1) a mistuned load mass using Euler buckled beams as negative stiffness corrector, 46 (2) an asymmetric excitation force, comprised of a harmonic component and a static component, applied to an isolator, 47,48 and (3) buckled beams under transverse excitation. 49,50 For instance, as shown by Rega et al. 49 and by Mayoof and Hawwa, 50 for buckled beams under transverse excitation, the initial curvature gives rise to the quadratic term and the mid-plane stretching gives rise to the cubic term. Increasing the curvature weakens the cubic nonlinearity, while its effect on the quadratic nonlinearity is non-monotonic. For the forcedeflection relationship of knitted spacer fabric in this study, the quadratic stiffness term represents force asymmetry. Using the Helmholtz-Duffing equation is a more . The pre-compression displacement x s is balanced out by the gravity force of the load mass. During vibration, the force transmitted by the spacer fabric is where x 2 is the sinusoidal displacement of the forced oscillation, and x 1 is the displacement of the load mass. A positive x 2 or x 1 means the motion is upward. At the steady state of vibration, the dynamic force balance leads to Let where x is the dynamic deformation of spacer fabric. The harmonic displacement of the forced oscillation can be expressed as x X t where ω is the angular frequency, so the excitation acceleration is Let the amplitude of excitation acceleration be G X = −ω 2 2 , which is also called the excitation level. Then, we have Furthermore, in order to account for the viscoelasticity of weft-knitted spacer fabric, a fractional derivative term aD x α is added into equation (4). Finally, to describe the vibration behavior of the mass-spacer fabric system under harmonic excitation, the proposed equation has the form of mx +cx +kx +k x +k x aD x = mG t 3 3  

Approximate analytical solution
The solution of the above equation mainly referred to the research work by Deng 42 and by Abolfathi. 51 The dynamic deformation of spacer fabric x depends on the excitation frequency, so that the frequency-domain solution of x is obtained by using the harmonic balance method 52 (HBM) with the first-order approximation here. To use the HBM, the coordinate system is transformed in order to change the asymmetric elastic force caused by the quadratic term into a symmetric elastic force containing a linear and cubic term, plus a constant force. 51,53 The transformation is shown as below.
The polynomial has a relation of equivalence as follows: which results in Let the new coordinate be z t x t ( ) ( ) = +δ , then after coordinate transformation, equation (5) becomes On the other hand, the fractional derivative of constant δ is where α ≥ 0 , and t −α is a decaying function. 35,41,54 Assume the time is adequately long to reach the steady At last, the equation of motion is transformed into Solving the unknowns k k k c a , , , , , 2 3 α in the original equation (5) is equivalent to solving the unknowns β κ α , , , , , k c a 3 in equation (11). The advantage of using the transformed form is that it can be treated as an isolator with symmetric elastic force excited by a harmonic force -mG t cos( ) ω and a constant force β . The resulting steady-state amplitude z t ( ) contains a harmonic term and a bias term. Assume the steady-state amplitude has the form in which A A jA r i = + . The harmonic amplitude is 2 A and the static displacement is A 0 . We use the harmonic amplitude versus excitation frequency curve to perform curve fit and identify model parameters. The steady-state response can also be written as where Using the HBM, to equate the constants on both sides of equation (11), we have The static displacement A 0 has only one real root, which is associated with A by Equation (19) can also be squared to obtain which is used to recover the harmonic amplitude ( 2 A ) versus excitation frequency curve after model parameters have been identified by curve fit. The static displacement A 0 can then be expressed using A referring to equation (16). Hence, the approximate analytical solution of the amplitude-frequency relationship obtained will be composed of a harmonic component and a static component.

Fitness function and goodness of fit
Given the expression of harmonic amplitude versus excitation frequency relationship, the next procedure is to find the optimal parameter estimates that best describe the vibration behavior of the mass-spacer fabric system by optimization strategy, which can be stated as where f x n ( ) is the fitness function that characterizes the harmonic amplitude versus excitation frequency relationship, and N is the volume of experimental data. For nonlinear vibrations, the harmonic amplitude may be a multivalued function for a certain range of excitation frequencies. Therefore, the fitness function can be built from equation (20), which has the form of The value of min ( ) x f x 2 2 is also the sum of square error (SSE). As a statistical indicator of the goodness of fit for a model while at the same time considering the effect of data volume, the root mean square error (RMSE) is given as Experimental data for the mass-spacer fabric system under 0.1 g, 0.2 g, and 0.3 g excitation level conditions were used for fitting. The raw data on a logarithmic scale in the frequency range of 10-25 Hz were converted into the linear scale format by interpolation. If the selected frequency range is too wide, the nonlinear resonance peak will have a lowered weighing factor during fitting, so the nonlinear feature at resonance may not be well captured. By optimizing the fitness function in equation (24), parameter estimates of β κ α , , , , , k c a 3 in the model can be solved. For optimization using MATLAB, a differential evolution (jDE) algorithm developed by Brest et al., 55 and Zhang and Sanderson 56 was adopted. The jDE algorithm has advantages over the commonly used method of least squares for nonlinearly constrained optimization. The method of least squares is a local optimization method, which depends on the choices of initial values and bounds of unknown model parameters. In contrast, the different evolution algorithm achieves solutions more efficiently and is free of choosing an initial value, and more importantly, it results in a higher level of the goodness of fit. were also considered, in contrast with using the model structure with parameters β κ α , , , , , k c a 3 . Table 2 summarizes their RMSEs to compare their fitting performances. Firstly, it is shown that the six-parameter model with β κ α , , , , , k c a 3 achieves the best curve fit by giving the smallest RMSEs. Secondly, using the six-parameter model with β κ α , , , , , k c a 3 results in reduced RMSEs for all three excitation level conditions, as compared with using the model without the fractional derivative term aD x α , that is, the four-parameter model with β κ , , , k c 3 . This makes clear the significance of the fractional derivative term for the model. Thirdly, using the model with the fractional derivative term aD x α alone, that is, the five-parameter model with β κ α , , , , k a 3 , gives better curve fit compared with using the model with the viscous damping coefficient c alone, that is, the four-parameter model with β κ , , , k c 3 . This again reveals the advantage of using the fractional derivative term to describe the vibration behavior of the mass-spacer fabric system. Last but not least, for 0.1 g and 0.3 g excitation level conditions, the viscous damping coefficient is redundant, since the five-parameter model with β κ α , , , , k a 3 gives the same RMSE as the six-parameter model with β κ α , , , , , k c a 3 ; however, for the 0.2 g excitation level condition, cx  and aD x α both helps improve the model structure. This is also demonstrated by the parameter estimates in Figure 4.

Results of curve fit
The physical significance of the fractional derivative term aD x α is some combination of the linear elastic force and the viscous damping force when 0 1 ≤ ≤ α . Thus, it contributes to both the elastic force and the damping force. When α = 0 , aD x α evolves into a linear spring ax ; when α = 1 , aD x α evolves into a viscous damper ax  . For the 0.3 g excitation level condition, however, the bestfit value of the fractional order α is larger than one. The extension of the interval ( 0 2 ≤ ≤ α ) makes the physical meaning of α more difficult to be defined. Despite this, it helps improve the model structure to some degree by giving a reduced RMSE.

Reconstructed force-displacement curve
To reconstruct the elastic force-displacement curves in the form of polynomials, the fractional derivative term aD x α is not considered, otherwise it would also contribute to the elastic force. The parameter estimates by curve fit using the four-parameter model with β κ , , , k c 3 are shown in Table 3, for the mass-spacer fabric system under the condition of G 0 = 0.1 g excitation level. The parameter estimates of k and k 2 from the original equation of motion are also shown. The SI units are used for all of parameters.
Spacer fabric is pre-stressed by a load mass. During harmonic excitation, the compressive displacement for spacer fabric changes dynamically. The emergence of the bias force β in equation (11) is due to the existence of an asymmetric stiffness term k x 2 2 in the system. The forcedisplacement relationship is y kx k x k x = + + 2 2 3

, that is
, where y represents the force and x represents the displacement which passes the origin ( , ) 0 0 of the coordinate plane, that is, the staticallyloaded position, and the point ( , ) − − δ β . In order to trace the frequency response curves (FRCs) to the nonlinear stiffness terms of the system, we let X x = +δ and Y y = + β , so the force-displacement relationship is recast into Y X k X = + κ 3 3 , as shown in Figure  6(a), which passes the point ( , ) δ β , that is the staticallyloaded position which is marked with a solid dot, and the origin ( , ) 0 0 , that is, the center of 180° rotational symmetry. The instantaneous stiffness at the statically-loaded position equals k , which is marked with a dashed line. It     the excitation frequency, as previously shown in Figure 5. It has removed the effect of the constant displacement δ caused by coordinate transformation. The negative value of the static displacement A 0 − δ in Figure 5 indicates that the center of oscillation shifts to the further compression direction for spacer fabric ( ∆ < X 0 and ∆ < Y 0 ), that is, the stiffness softening part. This is reasonable since a softening system undergoes a larger displacement than a hardening system under an identical level of force. The center of oscillation will shift to the softening part to balance out the difference. To comply with practice, the reconstructed elastic force-displacement curve in Figure  6(a) is also converted into the compression force versus compression displacement curve in Figure 6(b), when the load mass is 2 kg.

Backbone curve and the loci of peaks of the FRCs
To study the effects of excitation level on the FRCs, the loci of peaks of the FRCs concerning the harmonic amplitude, the static displacement and the acceleration transmissibility, and the backbone curves concerning the harmonic amplitude and the static displacement were used in order to track the trends of peak values and the resonance frequency. The derivations are shown as below, which referred to the work by Abolfathi. 51 At resonance, we have d d A φ = 0 . So the phase angle at resonance is φ = −90  . Using and Q c a = + ω ω π α α sin( ) 2 from equations (14), (20), and (21), at resonance, we will have is also noted that the tangent at the origin equals κ . The physical significance of ∆X is the change of compressive displacement for spacer fabric. To the positive direction ( ∆ > X 0 ) of the reference coordinate ( , ) δ β is the pressure relaxation ( ∆ > Y 0 ) process, while to the negative direction ( ∆ < X 0 ) is the further compression ( ∆ < Y 0 ) process. Due to the asymmetric force-displacement relationship about the statically-loaded position, displacements to the left and to the right exhibit stiffness softening and stiffness hardening characteristics separately, as is shown in Figure 6(a). Consequently, the center of harmonic displacement is no longer at the statically-loaded position. We consider the offset distance of the oscillation center from the statically-loaded position to be called the static displacement A 0 − δ , the value of which depends on As a result, the locus of harmonic amplitude at resonance is It is observed that the locus of harmonic response is not only a function of the excitation frequency, but also a function of the excitation level, load mass, viscous damping coefficient and the fractional derivative term. It does not depend on the stiffness of the system. Substituting equation (28) into equation (16), the locus of static displacement is solved accordingly. The transmitted force is composed of a harmonic component and a static component Since we have now If we let G = 0 , we obtain that P = 0 . By solving the equation and substituting A 0 in the form of A , the relation between the harmonic amplitude 2 A and the excitation frequency ω is obtained. This is the backbone curve of harmonic amplitude.

Effect of excitation level G
When the excitation level is very small, the system oscillates with a very limited amplitude around the statically loaded position. The tangential stiffness at the statically loaded position is the linear stiffness k . When the harmonic amplitude becomes infinitesimal, the backbone curve intercepts with the horizontal axis at the frequency , which equates to a linear system. The trend of the backbone curve leans to the left first and then to the right as the harmonic amplitude rise, in which case the nonlinear stiffness terms exhibit their influences. Positive cubic stiffness causes hardening behavior in the FRCs, characterized by peak bending to the right. Cubic stiffness result in symmetrical vibration response for displacements from the statically-loaded position to two opposite directions. However, the quadratic stiffness result in biased vibration response, and it is also the cause of the softening resonance peak. Aside from the primary harmonic amplitude, a static displacement component also exists for such a system. As shown in Figure 7(a) and (b), as the excitation level increases from 0.01 g to 0.1 g, that is, from 0 1 0 . G to G 0 , the amplitude curves bend to the left, an indication of stiffness softening. The equivalent stiffness in a softening system is lower than the tangential stiffness at the statically loaded position (i.e. the linear stiffness k ), as a result, the resonance frequency falls below f r . As the excitation level further increases from 0.1 g to 0.3 g, that is, from G 0 to 3 0 G , the amplitude curves bend to the right, an indication of stiffness hardening, and the corresponding resonance frequency is above f r . Figure 7(c) shows the corresponding transmissibility curves. Peak transmissibility is positively related with resonance frequency, which is in accordance with equation (33) for the expression of the locus of peak transmissibility.

Numerical simulation
Given the model parameters β κ α , , , , , k c a 3 identified under the 0.1 g excitation level condition as displayed in Figure 4(a), β and κ can be converted into the original unknowns k and k 2 in the governing equation (5), by using the transformation relationship in equation (7). The parameter values of k k k c a , , , , , 2 3 α obtained by the analytical method are brought into the MATLAB/Simulink block diagram as shown in Figure 8. Numerical simulation is carried out using the ode3 solver with a fixed step size 0.001 s. The simulation time is 100 s to achieve the steady state response. Besides, the fractional-order operator is approximated using the fourth-order Oustaloup filter. 57,58 Thus, the time-domain steady state response of the governing equation (5) is obtained. Time series of the dynamic displacement of the mass-spacer fabric system is then converted into the harmonic amplitude versus excitation frequency curve by using Fourier transform. Figure 9 shows the frequency response curve (FRC) obtained by the numerical method, in comparison with the experimental data and the analytical approximation of the harmonic amplitude 2 A versus frequency curve. Before evaluating the difference between the numerical result and the analytical approximation, it is also necessary to examine whether the analytical model has a good predictability. For this, model parameters identified by fitting the analytical model with experimental data under the 0.1 g excitation level condition are used to predict the FRCs for the 0.2 g and 0.3 g excitation level conditions. In the predicted analytical curves (in dotted "•" line) for 0.2 g and 0.3 g excitation levels as shown in Figure 9, the resonance peaks bend to the right exhibiting a hardening nonlinearity, which is different from the experimental results data (in circled "○" line) in which the resonance peaks bend to the left exhibiting a softening nonlinearity. On the other hand, in regions away from resonance, especially for high excitation frequencies, analytical predictions agree with the experimental results well. In these non-resonant regions, the mass-spacer fabric system undergoes small oscillations. In this case, the fitted nonlinear stiffness coefficients can correspond well with the elastic forcedisplacement relationship in this localized range. However, at resonance, the system undergoes the largest oscillations, especially for large excitation level conditions. In this case, the current analytical model becomes unreliable readily. To effectively predict the nonlinear feature at resonance, the model structure needs to be improved in future studies. level at resonance, an evident difference occurs between the numerical result and the analytical approximation. However, this difference is absent for the model eliminating the fractional derivative term. For excitation frequencies larger than 10 Hz, a discrepancy between the numerical results and the analytical approximations occurs. To explain this discrepancy, periodic solutions for the dynamic displacement of spacer fabric in the time domain are obtained. Time-domain numerical representations for spacer fabric under the 0.1 g excitation level are shown in Figure 10, which includes the time series, phase portrait and Fourier amplitude spectrum for the steadystate solutions when the excitation frequency varies from 8 to 14 Hz. The time series diagram contains the displacement of vibration platform and the displacement of fabric deformation. Phase portrait depicts the trajectories of a dynamic system in the state space, in which the horizontal and vertical axes represent the state variables of displacement and velocity. The closed trajectory in the phase portrait is a limit cycle. For a linear system under sinusoidal excitation force, the limit cycle appears as an oval. However, the limit cycle becomes distorted in a nonlinear system as is shown here. The peak at the driving frequency Ω and the peak at zero frequency in the Fourier amplitude spectrum represents the primary harmonic and the static displacement, respectively. The magnitudes of peak h and peak s are both marked in the spectrum. The peak s being an absolute value does not reflect the real sign of the static displacement. In fact, the static displacement is negative in this case, as is observed in Figure 7(a), indicating the center of oscillation shifts to the further compression direction for spacer fabric, that is, to the stiffness softening region.

Numerical versus analytical
The asymmetric period shape shown in the time series and the symmetry-breaking trajectory in the phase portrait suggest an even order harmonic distortion, which is caused by the quadratic stiffness term in the elastic force. This even order harmonic is represented by a peak at frequency 2Ω, that is, twice the driving frequency, in the Fourier amplitude spectrum. As the frequency increases from 8 to 12 Hz, this peak becomes stronger. Then, it fades with further increase of excitation frequency. This may cause the discrepancy between the approximate analytical solution and the numerical solution occurring at around 10 Hz excitation frequencies in Figure 9. Assume the governing equation of motion is changed from equation (5) to mx +cx +kx +k x +k x aD x = mG t 3 3   It results in a symmetric model with parameters k k k c a , , , , , 3 5 α instead of the original asymmetric model with parameters k k k c a , , , , , 2 3 α . Due to the fact that the stiffness terms in the equation of motion are all odd, the steady-state solutions will only present harmonics of odd order in the Fourier spectrum. Besides, the DC component in the asymmetric model, identified as a peak at zero frequency in the Fourier amplitude spectrum, will disappear for the symmetric model, in which case the numerical result will be in consistence with the approximate analytical solution at excitation frequencies around 10 Hz.

Conclusion
To promote the use of knitted spacer fabric for the protection of human body from vibration exposure, this study aims at a comprehensive understanding of the vibration behavior of the mass-spacer fabric system under harmonic excitation, with the use of experimental, analytical and numerical methods.
Experimental results of the acceleration transmissibility curve of the mass-spacer fabric system showed that increasing the vibration level gives rise to a nonlinear softening type of transmissibility curve, and it also results in a broadened isolation region. The nonlinear softening phenomenon in vibration is correlated with the nonlinear compressive force-displacement relationship of spacer fabric. As a result, to build an analytical model to describe the periodic response of the mass-spacer fabric system under harmonic excitation, a linear mass-spring-damper vibration model can only be acceptable for very small excitation levels. For large excitation levels, however, force nonlinearity should be considered. In our study, the analytical model uses an asymmetric polynomial composed of linear-quadratic-cubic stiffness terms to describe the nonlinear force of knitted spacer fabric. Besides, a fractional derivative damping term is used to account for the  viscoelasticity, which outperforms viscous damping by giving a higher level of the goodness of fit. The frequencydomain solution to the governing equation of motion was obtained using harmonic balance method (HBM).
Parameter analysis shows that frequency response curves (FRCs) can be softening, hardening and mixed softening-hardening depending on the magnitude of excitation level. Under small excitation levels, the linear stiffness is responsible for system dynamics. Under larger excitation levels, nonlinear stiffness terms start to exert their influences. With the increase of excitation level, resonance peak becomes softening first and then hardening, due to the quadratic and cubic stiffness terms, respectively. On one hand, the cubic stiffness causes hardening behavior in the FRCs, characterized by peak bending to the right. Cubic stiffness results in symmetrical vibration response for displacements away from the statically-loaded position to two opposite directions. On the other hand, the quadratic stiffness causes softening behavior in the FRCs, characterized by peak bending to the left. Quadratic stiffness also results in biased vibration response and causes an even order harmonic distortion in the numerical simulation.
Stemming from an empirical observation of nonlinear vibration by experiment, the proposed model cannot directly reveal the inherent mechanism of how material and structural properties of knitted spacer fabric act on the vibration dynamics. However, this could be done in future by performing a correlation analysis between model parameters and experimental variables, including material properties and structural properties such as the initial curvature of the three-dimensionally buckled monofilament. The current mathematical model provides a theoretical background for the development of weft-knitted spacer fabric as vibration isolator in the future. Using this model, parameter estimates from one excitation level condition can be used to predict the vibration behavior for a different excitation level condition when the excitation frequency is not at resonance. However, for the prediction of the nonlinear vibration at resonance, this model needs to be improved in future studies.