Jacek Tejchman
email: tejchmk@pg.gda.pl

Wei Wu

Modeling of textural anisotropy in granular materials with stochastic micro-polar hypoplasticity

Keywords: finite element method, granular material, micro-polar hypoplasticity, random field, shear zone, textural anisotropy

The paper deals with a numerical analysis of the effect of textural anisotropy on the behaviour of cohesionless granular materials with consideration of shear localization. For a simulation of the mechanical behavior of a granular material during a monotonic deformation path, an isotropic micropolar hypoplastic constitutive model was used. To describe textural effects, spatially correlated random fields of the initial void ratio were subject to rotation against the horizontal axis. The 2D random fields were generated using a conditional rejection method. The results were compared with those obtained with an anisotropic micro-polar constitutive model for an uniform distribution of the initial void ratio. The calculations were carried out with a dense granular specimen during plane strain compression under constant lateral pressure.

INTRODUCTION

Granular materials are heterogeneous and discrete systems composed of grains with different shape, size and roughness. Thus, their behaviour is influenced by the orientation of grains with respect to the direction of sedimentation [START_REF] Boehler | On yielding of oriented solids[END_REF][START_REF] Kanatani | Distribution of directional data and fabric tensor[END_REF][START_REF] Oda | Stress induced anisotropy in granular masses[END_REF][START_REF] Satake | A theory on the stress-induced anisotropy[END_REF], [START_REF] Khidas | Acoustic measurements of anisotropic elasticity in glass beads packings[END_REF]. This inherent anisotropy due to texture (fabric) is called a transverse isotropy since the material has a rotational symmetry with respect to one of the co-ordinates axes. The plane perpendicular to the orientation direction is called the bedding plane and it is a plane of isotropy.

The laboratory experiments in Fig. 1 show evidently that the orientation of the bedding plane relative to the principal stress directions has a pronounced effect on the stress-strain behaviour [START_REF] Arthur | Homogeneous and layered sand in triaxial compression[END_REF], [START_REF] Lam | Effect of initial anisotropic fabric and on strength and deformation characteristics of sand[END_REF][START_REF] Tatsuoka | Strength anisotropy and shear band direction in plane strain tests of sand[END_REF][START_REF] De Borst | Simulation of strain localization: a reappraisal of the Cosserat continuum[END_REF], 1994, 1997[START_REF] Abelev | Effects of cross anisotropy on 3-dimensional behaviour of sand. II: Stress-strain behaviour and shear banding[END_REF]. The shear stiffness, peak friction angle and average volume change are usually larger and strain corresponding to the peak value are smaller for loading perpendicular to the bedding plane than for loading parallel to it. The inclination of the shear zone with respect to the bottom becomes smaller. For large monotonic shearing, the stress ratio approaches a stationary value (Yamada and Ishihara 1979, [START_REF] Tatsuoka | Testing methods and results of element tests and testing conditions of plane strain model bearing capacity tests using air-dried dense Silver Leighton Buzzard sand[END_REF], i.e. anisotropy vanishes at residual state (critical state) at large shear deformation due to the so-called SOM-effect (swept out of memory effect) [START_REF] Gudehus | Attractors, percolation thresholds and phase limits of granular soils[END_REF]). The DEM simulations seem to confirm also this SOM-effect in granular bodies [START_REF] Pena | Investigation of the asymptotic states of granular materials using a discrete model of anisotropic particles[END_REF].

For describing the behaviour of granular materials within continuum mechanics, mainly elasto-plastic models [START_REF] Lade | Elasto-plastic stress-strain theory for cohesionless soil with curved yield surfaces[END_REF], Vermeer 1982, Pestana and Whittle 1999) and hypoplastic models [START_REF] Kolymbas | A rate-dependent constitutive equation for soils[END_REF], Wu 1992[START_REF] Darve | Yield surfaces and principle of superposition revisited by incrementally non-linear constitutive relations[END_REF], von Wolffersdorff 1996[START_REF] Gudehus | A comprehensive constitutive equation for granular materials[END_REF][START_REF] Bauer | Calibration of a comprehensive hypoplastic model for granular materials[END_REF][START_REF] Chambon | Incremental behaviour of a simple deviatoric constitutive CLoE model[END_REF][START_REF] Lanier | A general formulation of hypoplasticity[END_REF]) are applied. To describe the shear zone formation (thickness, inclination and spacing), however, all these models have to be enriched by a characteristic length of microstructure. This can be done by the following approaches: micro-polar theory (Mühlhaus1990, de Borst 1991, Tejchman and Wu 1993, Tejchman and Gudehus 2001, [START_REF] Maier | Numerische Modellierung der Entfestigung im Rahmen der Hypoplastizität[END_REF][START_REF] Huang | A polar extension of hypoplastic model for granular material with shear localization[END_REF][START_REF] Gudehus | Evolution of shear bands in sand[END_REF]), non-local theory [START_REF] Bazant | Yield limit degradation: non-local continuum model with local strain[END_REF], Pijaudier-Cabot and Bazant 1987, [START_REF] Brinkgreve | Geomaterial models and numerical analysis of softening[END_REF][START_REF] Chen | Non-local effects on dynamic damage accumulation in brittle solids[END_REF], [START_REF] Maier | Numerische Modellierung der Entfestigung im Rahmen der Hypoplastizität[END_REF][START_REF] Tejchman | Influence of a characteristic length on shear zone thickness in hypoplasticity with different enhancements[END_REF]), second-gradient theory [START_REF] Aifantis | On the microstructural origin of certain inelastic models[END_REF][START_REF] Aifantis | Update on class of gradient theories[END_REF], de Borst and Mühlhaus 1992, [START_REF] Sluys | Wave propagation, localisation and dispersion in softening solids[END_REF][START_REF] Chen | A non-local damage model for elastoplastic materials based on gradient plasticity theory[END_REF][START_REF] Pamin | Gradient-dependent plasticity in numerical simulation of localisation phenomena[END_REF]) and viscosity models [START_REF] Neddleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Loret | Dynamic strain localisation in elasto-visco-plastic solids, Part 1. General formulation and one-dimensional examples[END_REF][START_REF] Sluys | Wave propagation, localisation and dispersion in softening solids[END_REF], Belytschko et al. 1994[START_REF] Lodygowski | Numerical modelling of localized fracture of inelastic solids in dynamic loading process[END_REF]. Due to the presence of a characteristic length, the enhanced models regularize the ill-posedness of the underlying incremental boundary value problem caused by strain-softening and localization. The differential equations of motion of a regularized problem do not change their elliptic type during quasi-static calculations and hyperbolic type during dynamic calculations. Moreover, the numerical solutions of regularized problems do not suffer from the pathological mesh sensitivity (de Borst et al. 1992). Another numerical technique which enables to obtain mesh-independent results with localized deformation is the so-called strong discontinuity approach allowing finite elements with displacement discontinuities [START_REF] Larsson | Finite-element analysis of localization of deformation and fluid pressure in elastoplastic porous medium[END_REF][START_REF] Regueiro | Plane strain finite element analysis of pressure sensitive plasticity with strong discontinuity[END_REF][START_REF] Simone | Continous-discountinous modeling of mode-I and mode-II failure[END_REF]).

In enhanced hypoplastic models, two approaches of the inherent textural anisotropy have been already proposed. The first one was suggested by Tejchman et al. (2007) for second-gradient hypoplasticity wherein a structure tensor was included to take into account the spatial orientation of the bedding plane. The model was based an idea of [START_REF] Boehler | On yielding of oriented solids[END_REF]. The structural tensor was introduced into the nonlinear part of the hypoplastic equation. In this way, both strength and dilatancy were increased for the major normal stress perpendicular to the bedding plane. However, since the effect of anisotropy diminished in critical states at large shear deformation, this model was unable to describe the anisotropic

A c c e p t e d m a n u s c r i p t

3 behaviour if the initial void ratio of granulates was equal to the critical void ratio. The second approach was proposed by [START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF] within the framework of micro-polar hypoplasticity. They introduced a scalar density function, which depends on the orientation of the major principal stress with reference to the bedding plane. This formulation can be used to describe the anisotropic behaviour at all densities. In both enhanced models, only two additional material parameters were needed.

In this paper, a novel approach is attempted to describe textural anisotropy with usual isotropic micropolar hypoplasticity. First, the initial void ratio in the granular specimen was distributed stochastically by using a random correlated Gaussian field. Next, this field was rotated by different angles to simulate the specimen preparation process in laboratory experiments, which was characterized by the angle between the filling and loading direction (Fig. 1). An isotropic micro-polar hypoplastic constitutive model was used, which is able to describe the essential properties of granular materials during shear localization in a wide range of pressures and densities. The effect of this rotation on the spontaneous shear zone formation in dense sand during plane strain compression under constant lateral pressure was numerically investigated with the finite element method. To our knowledge, such FE-investigations have not been performed before. In addition, the numerical results from the isotropic model were compared with three different anisotropic micro-polar constitutive models for an uniform distribution of the initial void ratio.

The present paper is organized as follows. In Section 2, a micro-polar hypoplastic model is briefly described. Section 3 deals with the simulation of discrete random fields. The information about the finite element discretisation and boundary conditions is given in Section 4. The numerical results with an isotropic and anisotropic hypoplastic constitutive law modelling textural anisotropy are presented in Section 6 and 7. Finally conclusions are given in Section 8. Throughout this paper compressive stress and shortening strain were taken as negative. A superposed circle indicates objective time derivation and a superposed dot indicates material time derivation of a variable.

MICRO-POLAR HYPOPLASTIC CONSTITUTIVE MODEL

Non-polar hypoplastic constitutive models [START_REF] Gudehus | A comprehensive constitutive equation for granular materials[END_REF][START_REF] Bauer | Calibration of a comprehensive hypoplastic model for granular materials[END_REF], von Wolffersdorf 1996) describe the evolution of the effective stress tensor depending on the current void ratio, stress state and rate of deformation by isotropic non-linear tensorial functions. The tensorial functions can be obtained according to the representation theorem by [START_REF] Wang | A new representation theorem for isotropic functions[END_REF]. The constitutive models were formulated by a heuristic process considering the essential mechanical properties of granular materials undergoing homogeneous deformations. A striking feature pertinent to hypoplasticity is that the constitutive equation is incrementally non-linear in deformation rate. The hypoplastic models are capable of describing a number of significant properties of granular materials, e.g. non-linear stress-strain relationship, dilatant and contractant volumetric change, stress level dependence, density dependence and material softening. A further feature of hypoplastic models is the inclusion of critical states, in which grain aggregate can deform continuously at constant stress and constant volume. In contrast to plastic models, a decomposition of deformation into elastic and plastic parts, the formulation of yield surface, plastic potential, flow rule and hardening rule are not needed. The hallmark of these models are their simple formulation and procedure for determining the material parameters with standard laboratory experiments. The material parameters are related to the granulometric properties of granular materials, such as grain size distribution curve, shape, angularity and hardness of grains [START_REF] Herle | Determination of parameters of a hypoplastic constitutive model from properties of grain assemblies[END_REF]. A further advantage lies in the fact that one single set of material parameters is valid for a wide range of pressures and densities. A characteristic length can be introduced into hypoplasticity by means of micro-polar, non-local and second-gradient theory [START_REF] Maier | Numerische Modellierung der Entfestigung im Rahmen der Hypoplastizität[END_REF][START_REF] Tejchman | Influence of a characteristic length on shear zone thickness in hypoplasticity with different enhancements[END_REF]). In this paper, we adopted the micro-polar theory. A micro-polar model makes use of rotations and couple stresses, which have clear physical meaning for granular materials. The rotations can be observed during shearing and remain negligible during homogeneous deformation [START_REF] Oda | Micro-fabric and couple stress in shear bands of granular materials[END_REF]). [START_REF] Pasternak | Cosserat continuum modelling of granulate materials[END_REF] have demonstrated that the additional rotational degree of freedom of a micropolar continuum arises naturally by mathematical homogenization of an originally discrete system of spherical grains with contact forces and contact moments.

A micro-polar continuum considers deformations at two different levels, i.e. micro-rotation at the particle level and macro-deformation at the structural level [START_REF] Schäfer | Versuch einer Elastizitätstheorie des zweidimensionalen ebenen Cosserat-Kontinuums[END_REF][START_REF] Mühlhaus | Continuum models for layered and blocky rock[END_REF]). For the case of plane strain each material point has three degrees of freedom: two translations and one independent rotation (Fig. 2). The gradients of the rotation are related to the curvatures, which are associated with the couple stresses. The presence of the couple stresses gives rise to a non-symmetric stress tensor and a characteristic length.

The constitutive relationship between the rate of stress rate, the rate of couple stress, the strain rate and the curvature rate can be generally expressed by the following two equations (

σ σ = (1)
and

o c i i k l i k l i 5 0 m G(e, ,m ,d ,k ,d ), σ = (2)
The Jaumann stress rate and Jaumann couple stress rate in the above equations are defined by

o ij ij ik kj ik kj w w σ σ σ σ • = - + (3)
and

o i i i k k k k i m m 0.5w m 0.5m w • = - + . ( 4
)
The functions F ij and G i in Eqs. ( 1) and ( 2) represent isotropic tensor-valued functions of their arguments; σ ij is the Cauchy stress tensor, m i is the couple stress vector, e denotes the current void ratio, d kl c is the polar strain rate and k i denotes the rate of curvature vector:

c c ij ij ij ij d d w w , = + - and c i , i k w = .
(

) 5
The rate of deformation tensor d ij and the spin tensor w ij are related to the deformation velocity v i as follows:

, 2 / ) v v ( d i , j j , i ij + = , 2 / ) v v ( w i , j j , i ij - = . x / () () i i , ∂ ∂ = (6)
For moderate stress level, the grains of granular materials can be reasonably assumed to be incompressible. In this case, the change of void ratio depends only on the strain rate:

kk e (1 e )d . * = + (8)
For the numerical calculations, the following micro-polar hypoplastic constitutive equation are considered [START_REF] Tejchman | Shearing of a narrow granular strip with polar quantities[END_REF]:

o ^^c ij kl k s i j k l k 5 0 f [ L ( ,m ,d ,k d ) σ σ = + ^c c 2 ij d ij kl kl k k 50 f N ( ) d d k k d ] σ + (9)
and

o ^ĉ c i k l k 50 s i kl k 50 m / d f [ L ( ,m ,d ,k d ) σ = + ĉ c c 2 i d i kl kl k k 50 f N ( m ) d d k k d ] + , ( 10
)
wherein the normalized stress tensor ^ij σ is defined by

^ij ij kk σ σ σ = (11)
and the normalized couple stress vector ^i m is defined by

i i kk 50 m m d σ ∧ = , (12)
wherein d 50 is the mean grain diameter of soil. The scalar factors f s =f s (e, σ kk ) and f d =f d (e, σ kk ) in Eqs. ( 9) and ( 10) describe the influence of density and stress level on the incremental stiffness. The factor f s depends on the granulate hardness h s , the mean stress σ kk , the maximum void ratio e i and the current void ratio e:

1 n s i i k k s i i s h 1 e e f ( )( ) ( ) nh e e h β σ - + = - (13)
with

i0 d 0 i 2 c0 d 0 1 1 e e 1 1 1 h ( ) . 3 e e c 3 c α - = + - - (14)
The relative density in the above expression involves the void ratio in critical state e c , the minimum void ratio e d (the densest packing) and the maximum void ratio e i (the loosest packing). In a critical state, granular material experiences continuous deformation while the void ratio remains unchanged. The current void ratio e is bounded by the two extreme void ratios e i and e d . Based on experimental observations, the void ratios e i , e d and e c are assumed to depend on the stress level σ kk (see Fig. 3):

n i i 0 k k s e e exp[ ( / h ) ], σ = -- (16
)

n d d 0 k k s e e exp[ ( / h ) ], σ = -- (17)
n c c 0 k k s e e exp[ ( / h ) ], σ = -- (18)
wherein e i0 , e d0 and e c0 are the values of e i , e d and e c at σ kk =0, respectively. For the functions L ij , N ij , L i c and N i c , the following specific expressions are used (Tejchman and Gudehus 1999):

^^2 c c kl k ij 1 ij ij kl k 50 L a d ( d m k d ), σ σ = + + (19) ^^ĉ 2 2 c kl k i 1 i 5 0 1 i k l k 5 0 L a k d a m ( d m k d ) σ = + + , ( 20
) ^* ij ij 1 ij N a ( ), σ σ = + (21) ĉ 2 i i 1 c
N a a m , =

where The constitutive relationship requires the following ten material parameters: e i0 , e d0 , e c0 , φ c , h s , β, n ,α, a c and d 50 . The parameters h s and n are estimated from a single oedometric compression test with an initially loose specimen (h s reflects the slope of the curve in a semi-logarithmic representation, and n its curvature).

^1 * * 1 1 2 k l l k a c c [1 cos(3 )], σ σ θ -= + + (23) ^^* * * ^^kl lm mk * * 1 . 5 pq pq 6 cos( 3 ) ( ) [ ] θ σ σ σ σ σ = - (24)
The parameters α and β can be determined from a triaxial or plane strain test with a dense specimen. The critical friction angle φ c can be determined from the angle of repose or measured in a triaxial test with a loose specimen. The parameters e i0 , e d0 , e c0 and d 50 are obtained from conventional index tests (e c0 ≈e max , e d0 ≈e min , e i0 ≈(1.1-1.5)e max ). The calibration procedure was given in detail by Herle and Gudehus (1999).

SIMULATION OF DISCRETE RANDOM FIELDS

As granular materials consist of a large number of discrete grains deposited in a random manner, the system inevitably shows heterogeneous behavior to some extent. A reasonable way to account for this heterogeneity in numerical analyses is to create specimens of granular materials as a random field.

Considering two-dimensional or spatial random problems the input data is a set of random fields describing the initial parameters of the engineering model. In this work, a spatially correlated distribution of the initial void ratio e o was assumed in the form of two-dimensional truncated Gaussian random field:

( )

1 2 1 ( , ) o o e e x x νβ = + (26)
where the mean value of the initial void ratio is 0 e , ν is the coefficient of variation, and The midpoint method was applied. The method approximated the random field in each finite element by a single random variable defined as the value of the field at its centre. Randomness of the initial void ratio was described by the following homogeneous correlation function [START_REF] Walukiewicz | Simulation of nonhomogeneous random fields for structural applications[END_REF][START_REF] Górski | Non-linear models of structures with random geometric and material imperfections simulation-based approach[END_REF])

1 2 2 1 1 2 2 1 2 1 2 ( , ) (1 )
(1 ), , ,

x x x x d x x K x x s e x e x λ λ λ λ -∆ -∆ = × + ∆ + ∆ (27)
u u k k n n n p p p n p       = = =             x x K K x K x x x K K . ( 28
)
The unknown vector u x was estimated from the following conditional truncated distribution

( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 2 T 1 2 1 1 det 2 exp 2 1 m m t u k c u c c u c f t t π - - - -   = - - - -   -   x x K x x K x x , (29)
where c K and c x are described as the conditional covariance matrix and conditional expected value vector:

1 11 12 22 21 c - = - K K K K K , (30)
( )

1 12 22 c u k k - = + - x x K K x x .
To generate a discrete scalar random field, we make use of a base scheme with random values. The scheme covers all nodes i, i.e. 1≤i≤MN with MN=M×N (see Fig. 3). The numerical procedure consists of three stages. First, the four-corner random values are generated using an unconditional method (Fig. 3a). Next, a propagation scheme is defined (dotted rectangle in Fig. 3) and all random variables in the 9 area are generated one after another using the conditional method (Fig. 3b). In the third stage, the base scheme is shifted appropriately and the next group of unknown random values is simulated (Fig. 3c).

The base scheme is translated to cover all the field nodes (Fig. 3d). Using this algorithm, some 2000 random samples of the initial void ratio are generated. For details, the reader is referred to [START_REF] Walukiewicz | Simulation of nonhomogeneous random fields for structural applications[END_REF] and [START_REF] Górski | Non-linear models of structures with random geometric and material imperfections simulation-based approach[END_REF].

For the sake of simplicity and for a significant reduction of the computation time, only single realizations of the random correlated fields were performed in the medium dense specimen (of the height of h=14 cm and the width of b=4 cm) with o e -=0. [START_REF] Wang | A new representation theorem for isotropic functions[END_REF]. The standard deviations should be assumed based on experimental data. Due to the lack of experimental data, however, we assumed two different standard deviations s d =0.05 (small scatter) and s d =0.10 (large scatter) (see Fig. 4). We assumed strongly correlated fields in the horizontal direction (λ x1 =1) and weakly correlated fields in the vertical one (λ x2 =3). A strongly correlated field implied a small variation of e o . Therefore, the density in the horizontal direction was assumed to be close to an uniform distribution. Next, the generated random fields were rotated by the angle θ=0 o , θ=45 o and θ=90 o against the horizontal axis x 1 (Fig. 4) to simulate the filling process in laboratory tests (θ=90 o -δ (Figs.1), θ -bedding plane inclination). Afterwards, the FE-simulations of shear localizations during plane strain compression were performed.

INPUT DATA AND FE-IMPLEMENTATION

FE-calculations of plane strain compression tests were performed with a sand specimen with initial dimensions of h o =0.14 m (height), b=0.04 m (width) and l=1.0 m (depth -due to plane strain calculations). The specimen size was similar to those in the experiments by Vardoulakis (1980). In total, 896 quadrilateral elements (0.25×0.25 cm 2 ) divided into 3584 triangular elements were used. The quadrilateral elements composed of four diagonally crossed triangles were used to avoid volumetric locking due to dilatancy effects [START_REF] Groen | Three-dimensional elasto-plastic analysis of soils[END_REF]). Linear shape functions for displacements and the Cosserat rotation were used. The height of the finite elements was not larger than five times mean grain diameter to properly capture shear localization [START_REF] Tejchman | Numerical simulation of shear band formation with a polar hypoplastic model[END_REF]. The integration was performed with one sampling point placed in the middle of each element.

A quasi-static deformation was imposed through a constant vertical displacement increment ∆u prescribed at nodes along the upper edge of the specimen. The boundary conditions of null shear stress are imposed at the top and bottom of the specimen. To preserve the stability of the specimen against horizontal sliding, the node in the middle of the top edge was kept fixed. The vertical displacement increments were chosen as ∆u/h=0.00005. Some 8000 steps were performed.

As the initial stress state, a K 0 -state with σ 22 =γ d x 2 and σ 11 =K 0 γ d x 2 was assumed in the specimen; x 2 is the vertical coordinate measured from the top of the specimen, γ d =16.5 kN/m 3 denotes the initial volume weight and K 0 =0.50 is the earth pressure coefficient at rest (σ 11 -horizontal normal stress, σ 22 -vertical normal stress). Next, a confining pressure of σ c =200 kPa was prescribed.

For the solution of a non-linear equation system, a modified Newton-Raphson scheme with line search was used. The global stiffness matrix was calculated with only linear terms of the constitutive equations (Eqs.1 and 2). The stiffness matrix was updated every 100 steps. In order to accelerate the convergence in the softening regime, the initial increments of displacements and rotations in each calculation step were assumed to be equal to those in the previous step. The procedure was found to yield sufficiently accurate solutions with fast convergence. The magnitude of the maximum out-ofbalance force at the end of each calculation step was found to be smaller than 1 kN of the applied total

A c c e p t e d m a n u s c r i p t

10 vertical force along the top of the granular specimen. Due to the non-linear terms in deformation rate and material softening, this procedure turned out to be more efficient than the full Newton-Raphson method. The iteration steps were performed using translational and rotational convergence criteria. For the time integration of stresses in finite elements, a one-step Euler forward scheme was applied. The calculations were carried out with large deformations and curvatures using the so-called "Updated Lagrangian" formulation [START_REF] Tejchman | Influence of a characteristic length on shear zone thickness in hypoplasticity with different enhancements[END_REF]).

FE-RESULTS

WITHIN ISOTROPIC MICRO-POLAR HYPOPLASTICITY

The FE-analyses were carried out with the material constants for the so-called Karlsruhe sand. A predominantly quartz sand of medium grain size with the following index properties: e i0 =1.30, e d0 =0.51, The normalized overall vertical force in the granular specimen is the largest for the rotation angle θ==0 o , smaller for θ=45 o and the smallest for the angle θ=90 o (as in the experiments, Fig. 1). The same normalized vertical force was obtained in the residual state at large deformation for all rotation angles. The friction angle at the peak calculated from the principal stresses σ 1 =P/(bl) and σ 2 =σ c One shear zone occurs inside of the specimen which crosses the specimen and whose location is caused by a stochastic distribution of the initial void ratio. The shear zone thickness and its inclination can be determined by using one of the following indicators: shear deformation, void ratio, modulus of deformation and Cosserat rotation [START_REF] Tejchman | FE-studies on the influence of initial void ratio, pressure level and mean grain diameter on shear localisation[END_REF]). The shear zone thickness is insignificantly influenced by the rotation angle θ (it slightly increases with increasing θ). By making use of Cosserat rotations the shear zone thickness is about 8-9 mm [ [START_REF] Darve | Yield surfaces and principle of superposition revisited by incrementally non-linear constitutive relations[END_REF][START_REF] Górski | Non-linear models of structures with random geometric and material imperfections simulation-based approach[END_REF][START_REF] Górski | Deterministic and statistical size effect during shearing of granular layer within a micro-polar hypoplasticity[END_REF]×d 50 ]. The thickness of the shear zone is in good agreement with the experiments by Vardoulakis (1980). The shear zone inclination against the bottom always decreases slightly with increasing angle θ from υ=54 o (θ=0 o ) down to υ=50 o (θ =90 o ).

e c0 =0.82, φ c =30 o , h s =190 MPa,
The evolution of shear zones during deformation is shown in Fig. 7. As can be ascertained from Fig. 7, a pattern of shear zones can be observed in the sand specimen at the early stage of deformation. With As compared with the anisotropic models, where anisotropy is characterized by a structural tensor, the present approach is physics motivated. Indeed, the deposition process of granular materials usually gives rise to stratified media with numerous fine strata. Although the granular material of each stratum is isotropic, the overall behavior of the granular body can be strongly anisotropic.

FE-RESULTS

WITHIN ANISOTROPIC MICRO-POLAR HYPOPLASTICITY

As can be seen from the last section, even an isotropic model can describe the effect of anisotropy properly, provided the specimen is generated by random fields and rotated by an angle corresponding to laboratory experiments. In what follows, the numerical results of three anisotropic hypoplastic models are presented (assuming a uniform distribution of the initial void ratio).

In the approach by Tejchman et al. (2007), in order to take into account anisotropic properties, Equation 1 was enhanced by a structure tensor B kl [START_REF] Bauer | Investigations of shear banding in an anisotropic hypoplastic material[END_REF]):

o ^^c ij kl k s i j k l k 5 0 f [ L ( ,m ,d ,k d ) σ σ = + ^c c 2 kl d ijkl kl kl kl kl k k 50 f B ( s )N ( ) d d k k d ] σ + , (33)
where B ijkl is the fourth order tensor: In turn, in the approach by Niemunis (2003aNiemunis ( , 2003b) ) and [START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF], the critical void ratio in Eq.( 20) is replaced by

1 3 2 3 2 3 ( 2 ) ( )( ) ijkl ij kl ik jl ik jl jl ik B ss s s η η η η δ δ η η δ δ = + - + + - + θ θ θ   -   -       . (36)
n c0 c k k s e e exp[ ( / h ) ] σ - = -- (37)
where

* 0 0 0 c i j c c i j e e e s σ ∧ - → = + ∆ . ( 38
)
In the above relationship, the normalized deviator

* ij σ ∧ = ij σ * / kk σ is multiplied by the unit dyadic tensor / ij ij kl kl s s s s → = . (39)
Obviously, Eq.38 takes into account the inclination of the deviatoric stress with reference to the bedding plane. The parameter Another possibility is to modify the granular hardness h s (similarly to co e -) which affects directly the stiffness factor in Eq.18: The FE-analyses were carried out with the material constants for so-called Karlsruhe sand: e i0 =1.30, e d0 =0.51, e c0 =0.82, φ c =30 o , h s =190 MPa, β=1, n=0.40, α=0.20, a c =a 1 -1 and d 50 =0.5 mm. The parameters n and α were slightly reduced to obtain a similar shear resistance as in the isotropic micro-polar model.

* s i j s s i j h h h s σ ∧ - → = + ∆ . ( 41
The numerical results for 3 modified anisotropic micro-polar hypoplastic models are depicted in Figs. 11 and12 (with λ x1 =1, λ x2 =1 and s d =0.05). The results are qualitatively similar and close to those presented by [START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF] and [START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF]. They are also similar as the results with an isotropic micro-polar constitutive law (Figs. 5 and6). The shear resistance decreases with increasing bedding angle, and the shear zone width increases and its inclination υ decreases with increasing bedding angle. The shear zone width is slightly different for each model (for the same angle θ) due to the various rate of softening influencing the shear zone formation (the smaller rate of softening, the wider the shear zone width, [START_REF] Tejchman | FE-studies on the influence of initial void ratio, pressure level and mean grain diameter on shear localisation[END_REF]). The difference is about 1-2 mm. The evolution of one load-displacement curve after the peak in the anisotropic model by [START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF] at θ=90 o (curve 'a' in Fig. 11A) significantly differs from the remaining curves after the peak (the calculation was broken at u 2 t /h o =0.10) due to the fact that shear localization bumps into the top boundary at the left corner (Fig. 12A,a). Thus, it directly influences the vertical force acting on the top edge. In the remaining cases, shear localization occurs inside of the material far from the horizontal boundaries and the load-displacement diagrams are qualitatively the same.

CONCLUSIONS

The following conclusions can be derived on the basis of FE-calculations of plane strain compression with an isotropic micro-polar hypoplastic constitutive model using rotated random spatially correlated fields of the initial void ratio:

• The anisotropic effect is realistically described. The larger the bedding plane inclination, the smaller the peak internal friction angle. The vertical strain corresponding to the peak increases with increasing bedding plane inclination. The anisotropic effect vanishes in the residual state.

• The shear zone thickness slightly decreases with increasing bedding plane inclination. The shear zone inclination decreases with increasing bedding plane inclination.

• The larger the scatter of the initial void ratio, the larger the effect of anisotropy. Therein, the shear zone width slightly decreases and the its inclination with the base slightly increases.

• The position of the shear zone depends strongly of the stochastic distribution of the initial void ratio.

• The results with isotropic micro-polar constitutive model and stochastic distribution of the initial void ratio are qualitatively similar to those with anisotropic micro-polar constitutive models and uniform distribution of the initial void ratio.

The numerical results based on isotropic and anisotropic models show that anisotropy can be described realistically by both approaches. Although both the isotropic and the anisotropic models are of phenomenological nature, the approach based on the isotropic model is to be favored for its clear The isotropic approach mimics the sedimentation process. In fact, both approaches represent some homogenization of an anisotropic medium. In the anisotropic models, the homogenization is assumed at the element level via the constitutive equation. In the approach based on the isotropic model, the homogenization is performed numerically at the structure level.

The FE-analyses will be continued. In the next working steps, selected representative random fields will be chosen using the idea of stratified sampling or Latin hypercube sampling method (Górski and Tejchman 2007b) whose feature is that the joint probability distribution functions for random variables are divided into intervals with equal probability. A careful calibration procedure of material constants will be performed to obtain a quantitative comparison with experiments carried out by [START_REF] Tatsuoka | Testing methods and results of element tests and testing conditions of plane strain model bearing capacity tests using air-dried dense Silver Leighton Buzzard sand[END_REF].

: a) λ x1 =1, λ x2 =3, s d =0.05, θ=0 o ,b) λ x1 =1, λ x2 =3, s d =0.05, θ=45 o , c) λ x1 =1, λ x2 =3, s d =0.05, θ=90 o , d) λ x1 =1, λ x2 =3, s d =0.10, θ=0 o , e) λ x1 =1, λ x2 =3, s d =0.

A c c e p t e d m a n u s c r i p t 4 It

 4 has been shown that hypoplastic constitutive models without a characteristic length cannot describe realistically shear localization (Tejchman and Wu 1993, Maier 2002, Huang et al. 2002).

A c c e p t e d m a n u s c r i p t 5 The

 5 rate of Cosserat rotation w c is defined by

A c c e p t e d m a n u s c r i p t 6 In

 6 the above equations, the granulate hardness h s represents a reference pressure similar to the atmospheric pressure, the coefficients α and β express the dependence on density and stress level respectively, and n denotes the compression coefficient. The multiplier f d represents the dependence on the relative density:

  with

  parameter φ c is the friction angle in critical state and the parameter θ denotes the Lode angle in the deviatoric plane at ii σ ∧ =1, and * ij σ ∧ denotes the deviatoric part of ij σ ∧ . The micro-polar parameter a c in A c c e p t e d m a n u s c r i p t 7 Eq. (22) can be correlated with the grain roughness. This correlation can be established by studying the shearing of a narrow granular strip between two rough boundaries (Tejchman and Gudehus 2001). It can be represented by a constant, e.g. a c =1-5, or connected to the parameter a 1 -1 , e.g. a c =(0.5-1.5)×a 1 -1 . The parameter a 1 -1 lies in the range of 3.0-4.3 for the usually critical friction angle.

A c c e p t e d m a n u s c r i p t 8 where ∆x 1

 81 and ∆x 2 are the distances between two field points along the horizontal axis x 1 and vertical axis x 2 , λ x1 and λ x2 are the decay coefficients characterizing a spatial variability of the specimen properties (i.e. describe the correlation between the random field points) while the standard deviation s d represents their scattering. The random fields were generated using a conditional rejection method proposed by[START_REF] Walukiewicz | Simulation of nonhomogeneous random fields for structural applications[END_REF]. A discrete random field was described by multidimensional random variables defined at mesh nodes. The field was represented by the random vector[START_REF] Abelev | Effects of cross anisotropy on 3-dimensional behaviour of sand. II: Stress-strain behaviour and shear banding[END_REF] m × x , and its mean value[START_REF] Abelev | Effects of cross anisotropy on 3-dimensional behaviour of sand. II: Stress-strain behaviour and shear banding[END_REF] m × x . The covariance function was replaced by the symmetric defined covariance matrix ( ) m m × K (which is always positive when a correlation function is properly defined). The random variable vector[START_REF] Abelev | Effects of cross anisotropy on 3-dimensional behaviour of sand. II: Stress-strain behaviour and shear banding[END_REF] m × x was divided into blocks consisting of the unknown[START_REF] Abelev | Effects of cross anisotropy on 3-dimensional behaviour of sand. II: Stress-strain behaviour and shear banding[END_REF]

A c c e p t e d m a n u s c r i p t

- 1

 1 β=1, n=0.50, α=0.30, a c =a 1 and d 50 =0.5 mm (Tejchman and Gudehus 2001).The numerical results for three different spatially correlated random fields of the initial void ratio (with λ x1 =1, λ x2 =3 and s d =0.05) rotated by θ=0 o , θ=45 o and θ=90 o are depicted in Figs.5 and 6. The normalized load-displacement curves are shown in Fig.5(P -resultant vertical force on the top, b=0.04 m -specimen width, u 2 -vertical displacement of the top, h=0.14 m -initial height of the sand body, l=1.0 m -specimen depth). The deformed FE-meshes with the distribution of the equivalent total strain Fig.6(ε ij -strain tensor). The darker the region, the higher equivalent total the strain ε -. The entire range of strain was divided into 20 different grey scales.

A c c e p t e d m a n u s c r i p t 11 further

 11 deformation, strain localization continues to localize within a single zone. The shear zone becomes well visible immediately after the peak.The effect of the standard deviation s d is demonstrated in Figs.8 and 9. The larger the parameter s d , the larger the effect of anisotropy. The internal friction angle decreases from 43.81 o at u 2 /h=0.017 (θ=0 o ) down to 41.21 o (θ=90 o ) at u 2 /h=0.021. The difference between the peak internal friction angles is about 5%. The shear zone thickness decreases by about 1 mm and the its inclination against the bottom increases by about 2 o as compared to results with s d =0.05.

A c c e p t e d m a n u s c r i p t 12 Herein

 12 η i (i=1, 2, 3) are material parameters related to anisotropy. If η i =1, the constitutive relation for an initially isotropic material (1) is recovered. To describe diminishing anisotropy in critical state due to the SOM-effect, the additional anisotropic parameters were assumed to depend on the density factor f d in (20) as follows: η 1 =η 10 (f d -1), η 2 =η 20 (f d -1), η 3 =η 30 (f d -1). The following anisotropic parameters are assumed: η 10 =0.8, η 20 =0.25 and η 30 =0[START_REF] Tejchman | FE-studies on shear localization in an anistropic micro-polar hypoplastic granular material[END_REF]. Due to the assumption of η 30 =0, the micropolar anisotropic hypoplastic model requires 2 additional constants.

  describes the degree of anisotropy of the critical void ratio 0 c e from the isotropic reference at p=0. Assuming that anisotropic effects diminish at large deformations, at t=0 and the rate of decay B are the only additional material constants. The evolution equation for 0 ∆ c e is objective since d ij c and k i are objective. The anisotropic parameters were assumed to be 0 ∆ c e (t=0)=0.25 and B=5 (Tejchman and Niemunis 2006).

) A c c e p t e d m a n u s c r i p t 13 After

 13 some trial calculations, the following values were found suitable: 0 ∆ c e (t=0)=1.0 and B=5.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t 18 66Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :

 181234 Fig.1: Method used to prepare a granular specimen (SLB sand) and relationships among the stress ratio σ 1 /σ 3 , the average shear strain γ=ε 1 -ε 3 and the average volumetric strain ε v =ε 1 +ε 3 for different angles of δ between 90 and 0 degrees at σ 3 =80 kPa[START_REF] Tatsuoka | Testing methods and results of element tests and testing conditions of plane strain model bearing capacity tests using air-dried dense Silver Leighton Buzzard sand[END_REF]

Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :

 5678 Fig.4: Two-dimensional distribution of initial void ratio e o in the granular specimen for the following paramters: a)λ x1 =1, λ x2 =3, s d =0.05, θ=0 o ,b) λ x1 =1, λ x2 =3, s d =0.05, θ=45 o , c) λ x1 =1, λ x2 =3, s d =0.05, θ=90 o , d) λ x1 =1, λ x2 =3, s d =0.10, θ=0 o , e) λ x1 =1, λ x2 =3, s d =0.10, θ=90 o (h o =14 cm -specimen initial height, b=4 cm -specimen initial width, θ -bedding plane inclination, s d -standard deviation, λ xidecay coefficients)

Fig. 9 :

 9 Fig.9: Deformed meshes with the distribution of equivalent total strain

Fig. 10 :A c c e p t e d m a n u s c r i p t 20 Fig. 11 :Fig. 12 :A c c e p t e d m a n u s c r i p t 21 FIGURE 1 9 A c c e p t e d m a n u s c r i p t 32 FIGURE 10 Θ

 1020111221193210 Fig.10: Normal vector s of the bedding plane (a) and its inclination θ (b) with respect to a fixed coordinate system (θ -bedding plane inclination, υshear zone inclination)


Z K Peng
email: z.peng@sheffield.ac.uk

Z Q Lang
email: z.lang@sheffield.ac.uk

S A Billings
email: s.billings@sheffield.ac.uk

Nonlinear output frequency response functions of MDOF systems with multiple nonlinear components

In engineering practice, most mechanical and structural systems are modeled as Multi-Degree-of-Freedom (MDOF) systems such as, e.g., the periodic structures.

When some components within the systems have nonlinear characteristics, the whole system will behave nonlinearly. The concept of Nonlinear Output Frequency Response Functions (NOFRFs) was proposed by the authors recently and provides a simple way to investigate nonlinear systems in the frequency domain. The present study is concerned with investigating the inherent relationships between the NOFRFs for any two masses of nonlinear MDOF systems with multiple nonlinear components. The results reveal very important properties of the nonlinear systems. These properties clearly indicate how the system linear characteristic parameters govern the propagation of the nonlinear effect induced by nonlinear components in the system. One potential application of the results is to detect and locate faults in engineering structures which make the structures behave nonlinearly.

Nomenclature

x(t), u(t) the output and input of the nonlinear system ( ) X jω , )

( ω j U the spectrum of the system output and input ) ,..., ( the n th order NOFRF M, C, K the system mass, damping and stiffness matrices m i , c i , k i the i th mass, damping and stiffness parameter

) (i L FS and ) (i L FD
the restoring forces of L(i) th nonlinear damper and stiffness ( )

l i L r , ) (

Γ

the term introduced by the nonlinear force

Introduction

In engineering practice, for many mechanical and structural systems, more than one set of coordinates are needed to describe the system behaviour. This implies a MDOF model is needed to represent the system. In addition, these systems may also behave nonlinearly due to nonlinear characteristics of some components within the systems. For example, a beam with breathing cracks behaves nonlinearly because of the cracked elements inside the beam [START_REF] Chondros | Vibration of a beam with breathing crack[END_REF]. For nonlinear systems, the classical Frequency Response Function (FRF) cannot achieve a comprehensive description for the system dynamical characteristics, which, however, can be fulfilled using the Generalised Frequency Response Functions (GFRFs) [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF]. The GFRFs, which are extension of the FRFs to the nonlinear case, are defined as the Fourier transforms of the kernels of the Volterra series [START_REF] Worden | A harmonic probing algorithm for the multi-input Volterra series[END_REF]. The Volterra series and its derivative GFRFs are powerful tools for the analysis of nonlinear systems and have been widely studied in the past two decades [START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF][5] [START_REF] Bartos | Sun, Characterization of Abrupt Nonlinearity by the Volterra-Fourier Method[END_REF]. The applications of the Volterra series range from the electrical engineering [START_REF] Fard | Synchronous generator model identification for control application using volterra series[END_REF]~ [START_REF] Nam | Volterra series representation of time-frequency distributions[END_REF], communications [START_REF] Salgado | Experimental validation of Volterra series nonlinear modelling for microwave subcarrier optical systems[END_REF]~ [START_REF] Lang | Evaluation of output frequency responses of nonlinear systems under multiple inputs[END_REF], network theory [START_REF] Stegmayer | Volterra series and neural networks to model an electronic device nonlinear behavior[END_REF][14] to structure dynamics [START_REF] Zheng | Non-linear frequency-domain analysis of jack-up platforms[END_REF]~ [START_REF] Peng | Crack detection using nonlinear output frequency response functions[END_REF].

If a differential equation or discrete-time model is available for a nonlinear system, the GFRFs can be determined using the algorithm in [START_REF] Billings | Spectral analysis for nonlinear system, part I: parametric non-linear spectral analysis[END_REF]~ [START_REF] Peyton Jones | A recursive algorithm for the computing the frequency response of a class of nonlinear difference equation models[END_REF]. However, the GFRFs are much more complicated than the FRF. GFRFs are multidimensional functions [START_REF] Zhang | Analysing non-linear systems in the frequency domain, I: the transfer function[END_REF] [START_REF] Zhang | Analysing nonlinear systems in the frequency domain, II: the phase response[END_REF], which can be difficult to measure, display and interpret in practice. Recently, the novel concept known as Nonlinear Output Frequency Response Functions (NOFRFs) was proposed by the authors [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF]. The concept can be considered to be an alternative

A c c e p t e d m a n u s c r i p t

3 extension of the FRF to the nonlinear case. NOFRFs are one dimensional functions of frequency, which allow the analysis of nonlinear systems in the frequency domain to be implemented in a manner similar to the frequency domain analysis of linear systems and which provide great insight into the mechanisms which dominate important nonlinear behaviours. Using the NOFRF, the authors have investigated the resonance phenomena for a class of nonlinear systems [START_REF] Peng | Resonances and resonant frequencies for a class of nonlinear systems[END_REF]. Most recently, the concept of the NOFRF has been extended from the SISO case to the MIMO case by the authors [START_REF] Peng | Nonlinear Output Frequency Response Functions for Multi-Input Nonlinear Volterra Systems[END_REF].

The present study is concerned with the analysis of the inherent relationships between the NOFRFs for any two masses of MDOF systems with multiple nonlinear components. The results reveal, for the first time, very important properties of the nonlinear systems. These properties clearly indicate how the system linear characteristic parameters govern the propagation of the nonlinear effect induced by nonlinear components in the system. One potential application of the results is to detect and locate faults in engineering structures which make the structures behave nonlinearly. This study will focus on the derivation and verification of nine important properties of nonlinear MDOF systems using a NOFRF based analysis.

MDOF Systems with Multiple Nonlinear Components

Figure 1, a multi-degree freedom oscillator A typical multi-degree-of-freedom oscillator is shown as Figure 1, the input force is added on the Jth mass.

If all springs and damping have linear characteristics, then this oscillator is a MDOF linear system, and the governing motion equation can be written as

) (t F Kx x C x M = + +
(1) where

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = n m m m M 0 0 0 0 0 0 2 1 m n k n m n-1 m 2 m 1 k 2 k 1 x 1 x 2 x n-1 x n u(t) c 1 c 2 c n
A c c e p t e d m a n u s c r i p t 4 is the system mass matrix, and

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - + - - + - - + = - - n n n n n n c c ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - + - - + - - + = - - n n n n n n k k k k k k k k k k k k k K 0 0 0 0 0 0 1 1 3 3 2 2 2 2 1
are the system damping and stiffness matrix respectively. ' 1 ) , , (

n x x x =
is the displacement vector, and

' 1 ) 0 , , 0 ), ( , 0 , , 0 ( ) ( J n J t u t F - - =
is the external force vector acting on the oscillator. Equation ( 2) is the basis of the modal analysis method, which is a well-established approach for determining dynamic characteristics of engineering structures [START_REF] He | Modal Analysis[END_REF]. In the linear case, the displacements

) (t x i ( n i , , 1
=
) can be expressed as

∫ +∞ ∞ - - = τ τ τ d u t h t x i i ) ( ) ( ) ( ) ( (2)
where ) (

) ( t h i ( n i , , 1 =
) are the impulse response functions that are determined by equation [START_REF] Chondros | Vibration of a beam with breathing crack[END_REF], and the Fourier transform of ) (

) ( t h i
is the well-known FRF.

Assume there are L nonlinear components, which have nonlinear stiffness and damping, in the MDOF system, and they are the

) (i L th ( L i , , 1
=
) components respectively, and the corresponding restoring forces ) (

) ( ∆ i L FS and ) ( ) ( ∆ i L FD
are the polynomial functions of the deformation ∆ and ∆ , i.e.,

( ) ∑ = ∆ = ∆ P l l l i L i L r FS 1 , ) ( ) ( ) ( , ∑ = ∆ = ∆ P l l l i L i L w FD 1 ) ), ( ( ) ( ) (
where P is the degree of the polynomial. Without loss of generality, assume 1 ) (i L and

n J i L , , 1 ) ( ≠ ( L i ≤ ≤ 1
) and

( ) 1 ), ( ) ( i L i L r k = and ( ) 1 ), ( ) ( i L i L w c = . Denote ( ) ' ) ( ) 1 ( n nf nf NF = (3) where ( ) 1 ), ( if 1 , 1 ) ( if 1 ), ( , 1 ) ( if 0 ) ( ) ( L i i L l L i i L l L i i L i L l NonF NonF l nf i L i L ≤ ≤ = ≤ ≤ - = ≤ ≤ - ≠ ⎪ ⎩ ⎪ ⎨ ⎧ - = ( n l ≤ ≤ 1 ) (4)
and

∑ ∑ = - = - - + - = P l l i L i L l i L P l l i L i L l i L i L x x w x x r NonF 2 ) ( 1 ) ( ) ), ( ( 2 )
( 1 ) ( ) ), ( ( ) ( ) ( ) ( (5)
Then the motion of the MDOF oscillator in Figure 1 can be determined by below equation 6) are the motion governing equations of nonlinear MDOF systems with multiple nonlinear components. The L nonlinear components can lead the whole system to behave nonlinearly. In this case, the Volterra series [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF] can be used to describe the relationships between the displacements

) (t F NF Kx x C x M + = + + (6)
) (t x i ( n i , , 1 =
) and the input force ) (t u as below

i j i i N j j j i i d t u h t x τ τ τ τ ) ( ) ,..., ( ) ( 1 1 1 ) , ( ∏ ∑ ∫ ∫ = = ∞ ∞ - ∞ ∞ - - = (7)
under quite general conditions [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF]. In equation ( 7), ) ,..., ( 1) , (

j j i h τ τ , ( n i , , 1 = , N j , , 1 =
), represents the jth order Volterra kernel for the relationship between u(t) and the displacement of m i .

When a system is linear, its dynamical properties can be easily analyzed using the FRFs defined as the Fourier transform of ) (

) ( t h i ( n i , , 1 =
) in equation ( 2). However, as equation [START_REF] Fard | Synchronous generator model identification for control application using volterra series[END_REF] shows, the dynamical properties of a nonlinear system are determined by a series of Volterra kernels, such as ) ,..., ( 1) , (

j j i h τ τ , ( n i , , 1 = , N j , , 1
=
) for the MDOF nonlinear systems considered in the present study. The objective of this paper is to study the nonlinear MDOF systems using the concept of Nonlinear Output Frequency Response Functions (NOFRFs), which is an alternative extension of the FRF to the nonlinear case and is derived based on the Volterra series approach of nonlinear systems.

Nonlinear Output Frequency Response Functions

The definition of NOFRFs is based on the Volterra series theory of nonlinear systems. The Volterra series extends the well-known convolution integral description for linear systems to a series of multi-dimensional convolution integrals, which can be used to represent a wide class of nonlinear systems [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF].

Consider the class of nonlinear systems which are stable at zero equilibrium and which can be described in the neighbourhood of the equilibrium by the Volterra series

i n i i n N n n d t u h t x τ τ τ τ ) ( ) ,..., ( ) ( 1 1 1 ∏ ∑ ∫ ∫ = = ∞ ∞ - ∞ ∞ - - = (8)
where x(t) and u(t) are the output and input of the system, ) ,..., (

n 1 n h τ τ
is the nth order Volterra kernel, and N denotes the maximum order of the system nonlinearity. Lang and Billings [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF] derived an expression for the output frequency response of this class of nonlinear systems to a general input. The result is This expression reveals how nonlinear mechanisms operate on the input spectra to produce the system output frequency response. In (9), )

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = ∀ = ∫ ∏ ∑ = + + = - = ω ω ω ω σ ω ω ω π ω ω ω ω n n n i i n n n n N n n d j U j j H n j X j X j X ,..., 1 1 1 1 1 ) ( ) ,..., ( ) 2 ( 1 ) ( for ) ( ) ( (9)
( ω j X is the spectrum of the system output, )

( ω j X n represents the nth order output frequency response of the system,

1 1 + + - ∞ ∞ - ∞ ∞ - ∫ ∫ = ( 10
)
is the nth order Generalised Frequency Response Function (GFRF) [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF], and . Equation ( 9) is a natural extension of the well-known linear relationship

∫ ∏ = + + = ω ω ω ω σ ω ω ω n n n i i n n d j U j j H ,..., 1 1 1 )
) ( ) ( ) ( ω ω ω j U j H j X = , where ) ( ω j H
is the frequency response function, to the nonlinear case.

For linear systems, the possible output frequencies are the same as the frequencies in the input. For nonlinear systems described by equation ( 9), however, the relationship between the input and output frequencies is more complicated. Given the frequency range of an input, the output frequencies of system (9) can be determined using the explicit expression derived by Lang and Billings in [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF].

Based on the above results for the output frequency response of nonlinear systems, a new concept known as the Nonlinear Output Frequency Response Function (NOFRF) was recently introduced by Lang and Billings [START_REF] Peng | Resonances and resonant frequencies for a class of nonlinear systems[END_REF]. The NOFRF is defined as

∫ ∏ ∫ ∏ = + + = = + + = = ω ω ω ω ω ω ω ω σ ω σ ω ω ω ω n n n n i i n n i i n n n d j U d j U j j H j G ,..., 1 ,..., 1 1 1 1
) ( ) ( ) ,..., ( ) ( [START_REF] Peddanarappagari | Volterra series approach for optimizing fiber-optic communications system designs[END_REF] under the condition that 0

) ( ) ( ,..., 1 1 ≠ = ∫ ∏ = + + = ω ω ω ω σ ω ω n n n i i n d j U j U ( 12
)
Notice that )

( ω j G n
is valid over the frequency range of )

( ω j U n
, which can be determined using the algorithm in [START_REF] Lang | Output frequency characteristics of nonlinear system[END_REF].

By introducing the NOFRFs )

( ω j G n , N n , 1 =
, equation ( 9) can be written as

∑ ∑ = = = = N n n n N n n j U j G j X j X 1 1 ) ( ) ( ) ( ) ( ω ω ω ω (13)
which is similar to the description of the output frequency response for linear systems. The NOFRFs reflect a combined contribution of the system and input to the system output frequency response behaviour. It can be seen from equation [START_REF] Peddanarappagari | Volterra series approach for optimizing fiber-optic communications system designs[END_REF] that ) . For a nonlinear system, the dynamical properties are determined by the GFRFs n H (n= 1,…,N). However, from equation [START_REF] Salgado | Experimental validation of Volterra series nonlinear modelling for microwave subcarrier optical systems[END_REF] it can be seen that the GFRF is multidimensional [21][22] which makes the GFRFs difficult to measure, display and interpret in practice. According to equation [START_REF] Peddanarappagari | Volterra series approach for optimizing fiber-optic communications system designs[END_REF] is that it is one dimensional, and thus allows the analysis of nonlinear systems to be implemented in a convenient manner similar to the analysis of linear systems. Moreover, there is an effective algorithm [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF] available which allows the estimation of the NOFRFs to be implemented directly using system input output data.

( ω j G n

Analysis of MDOF Systems with Multiple Nonlinear Components Using NOFRFS

GFRFs of MDOF Systems with Multiple Nonlinear Components

From equation ( 6), the GFRFs ) ,...,

( j j i j j H ω ω , ( n i , , 1 = , N j , , ( 1 ) ,
= 1
) can be determined using the harmonic probing method [START_REF] Billings | Spectral analysis for nonlinear system, part I: parametric non-linear spectral analysis[END_REF] [START_REF] Billings | Mapping nonlinear integro-differential equations into the frequency domain[END_REF].

First consider the input ) (t u is of a single harmonic

t j e t u ω = ) ( (14)
Substituting [START_REF] Hussain | A new neural network structure for temporal signal processing[END_REF] and

t j i i e j H t x ω ω) ( ) ( ) 1 , ( = ( n i , , 1 = ) (15)
into equation [START_REF] Bartos | Sun, Characterization of Abrupt Nonlinearity by the Volterra-Fourier Method[END_REF] and extracting the coefficients of t j e ω yields, ( )

T J n J j H K jC M ) 0 0 1 0 0 ( ) ( 1 1 2 - - = + + - ω ω ω (16) where ( ) T n j H j H j H ) ( ) ( ) ( ) 1 , ( ) 1 , 1 ( 1 ω ω ω = (17)
From equation ( 16), it is known that ( ) It is obtained from equations ( 18)~(20

T J n J K jC M j H ) 0 0 1 0 0 ( ) ( 1 1 2 1 - - - + + - = ω ω ω (18) Denote K jC M j + + - = Θ ω ω ω 2 ) ( (19) and
⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = Θ - ) ( ) ( ) ( ) ( ) ( ) , ( ) 1 , ( ) , 1 ( ) 1 , 1 ( 1 ω ω ω ω ω j Q j Q j Q j Q j n n n n (20)
) that ) ( ) ( ) , ( ) 1 , ( ω ω j Q j H J i i = ) , , 1 ( n i = (21
) Thus, for any two consecutive masses, the relationship between the first order GFRFs can be expressed as

) ( ) ( ) ( ) ( ) ( 1 , 1 ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( ω λ ω ω ω ω + + + = = i i J i J i i i j Q j Q j H j H ) 1 , , 1 ( - = n i ( 22
)
The above procedure used to analyze the relationships between the first order GFRFs can be extended to investigate the relationship between the N th order GFRFs with 2 ≥ N . To achieve this, consider the input

∑ = = N k t j k e t u 1 ) ( ω

Substituting this input and

+ + + + + = + + t j N N i t j N i t j i i N N e j j H N e j H e j H t x ) ( 1 ) , ( ) 1 , ( 1 ) 1 , ( 1 1 ) , , ( ! )
( ) ( ) ( ω ω ω ω ω ω ω ω ( n i , , 1 = ) (23)
into equation ( 6), and, for the first row of equation ( 6), extracting the coefficients of

t j N e ) ( 1 ω ω + + yields ( ) ( ) 0 ) , , ( ) ( ) , , ( ) ( ) )( ( ) ( 1 ) , 2 ( 2 1 2 1 ) , 1 ( 2 1 1 2 1 2 1 1 = + + + - + + + + + + + + - N N N N N N N j j H k jc j j H k k c c j m ω ω ω ω ω ω ω ω ω ω (24)
Similarly, it can be easily deduced that, for the masses that are not connected to the th

) (i L ( L i , , 1
=
) spring, the GFRFs satisfy the following relationships

( ) ( ) 0 ) , , ( ) ( ) , , ( ) ( ) ( 1 ) , 1 ( 1 1 ) , ( 1 2 1
= + + + - + + + + + + - - N N n n N n N N n n N n N n j j H k jc j j H k jc m ω ω ω ω ω ω ω ω ω ω (25) ( ) (
)

) , , ( ) ( ) , , ( ) )( ( ) ( 1 ) , 1 ( 1 1
) , ( 1 1 1 2 1 N N i i N i N N i i i N i i N i j j H k jc j j H k k c c j m ω ω ω ω ω ω ω ω ω ω - + + + + + - + + + + + + + + - ( ) 0 ) , , ( ) ( 1 ) , 1 ( 1 1 1
= + + + - + + + N N i i N i j j H k jc ω ω ω ω ( n l L l L i ), ( , 1 ) ( , 1 - ≠ , L l ≤ ≤ 1
) (26) For the mass that is connected to the left of the ) (i L th spring, the GFRFs satisfy the following relationships ( )

) , , ( ) ( ) , , ( ) ( ) )( ( ) ( 1 ) , 2 ) ( ( 1 ) ( 1 1 ) ( 1 ) , 1 ) ( ( 1 ) ( 1 ) ( 1 ) ( 2 1 1 ) ( N N i L i L N i L N N i L L i L N i L i L N i L j j H k jc j j H i k k c c j m ω ω ω ω ω ω ω ω ω ω - - - - - - - + + + - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + + + + + + - ( ) 0 ) , , ( ) , , ( ) ( 1 ) ( , 1 ) ( 1 ) ), ( ( ) ( 1 )
( = Λ + + + + - - N i L i L N N N i L i L N i L j j j j H k jc ω ω ω ω ω ω

A c c e p t e d m a n u s c r i p t

9

( L i ≤ ≤ 1
) ( 27) For the mass that is connected to the right of the ) (i L th spring, the GFRFs satisfy the following relationships ( )

) , , ( ) ( ) , , ( ) )( ( ) ( 1 ) , 1 ) ( ( ) ( 1 ) ( 1 ) ), ( ( 1 )
( ) ( 1 1 ) ( ) ( 2 1 ) ( N N i L i L N i L N N i L i L i L N i L i L N i L j j H k jc j j H k k c c j m ω ω ω ω ω ω ω ω ω ω - + + + + + - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + + + + + + - ( ) 0 ) , , ( ) , , ( ) ( 1
) ( , 1 ) ( 1 ) , 1 ) ( ( 1 ) ( 1 1 ) ( = Λ - + + + - - + + + N i L i L N N N i L i L N i L j j j j H k jc ω ω ω ω ω ω ( L i ≤ ≤ 1 ) (28)
In equations ( 27) and ( 28

), ) , , ( 1 )
( , 1 ) ( N i L i L N j j ω ω - Λ
represents the extra terms introduced by

) (i L
NonF for the N th order GFRFs, for example, for the second order GFRFs, ( ) ( )

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , ( 1 ) 1 ), ( ( 2 ) 1 , 1 ) ( ( 2 ) 1 ), ( ( 1 ) 1 , 1 ) ( ( 2 ) 1 ), ( ( 1 ) 1 ), ( (
2 ) 1 , 1 ) ( ( 1 ) 1 , 1 ) ( ( ) 2 ), ( ( 2 1 ) 2 ), ( ( 2 1 ) ( , 1 ) ( 2 ω ω ω ω ω ω ω ω ω ω ω ω j H j H j H j H j H j H j H j H r w j j i L i L i L i L i L i L i L i L i L i L i L i L - - - - - - - + + - = Λ ( L i ≤ ≤ 1 ) (29)
Denote

( ) T N N n N N N N j j H j j H j j H ) , , ( ) , , ( ) , , ( 1 ) , ( 1 ) , 1 ( 1 ω
ω ω ω ω ω = (30)
and

[ ] T N N N N n A A j j A ) ( ) 1 ( ) , , ( 1 = ω ω (31)
where

( ) 1 ), ( if 1 , 1 ) ( if 1 ), ( , 1 ) ( if ) , , ( ) , , ( 0 1 ) ( , 1 ) ( 1
) ( , 1 ) ( L i i L l L i i L l L i i L i L l j j j j l A N i L i L N N i L i L N N ≤ ≤ = ≤ ≤ - = ≤ ≤ - ≠ ⎪ ⎩ ⎪ ⎨ ⎧ Λ Λ - = - - ω ω ω ω ( n l ≤ ≤ 1 ) (32)
then equations ( 25)~( 28) can be written in a matrix form as ) , , ( ) , , ( )) ( (

1 1 1 N N N N N j j A j j H j ω ω ω ω ω ω = + + Θ (33) so that ) , , ( )) ( ( ) , , ( 1 1 1 1 N N N N N j j A j j j H ω ω ω ω ω ω + + Θ = - (34)
Therefore, for each mass, the N th order GFRF can be calculated as

∑ = - - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Λ Λ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + + = L l N l L l L N N l L l L N T N l L i N l L i N N i j j j j j Q j Q j j H 1 1 ) ( , 1 ) ( 1 ) ( , 1 ) ( 1 ) ( , 1 1 ) ( , 1 ) , ( ) , , ( ) , , ( ) ( ( )
) ( ( ) , , ( ω ω ω ω ω ω ω ω ω ω ) , , 1 ( n i = (35) Define ) , , ( ) , , ( ) , , ( 1 ) , 1 ( 1 ) , ( 1 1 ,
N N i N N i N i i N j j H j j H ω ω ω ω ω ω λ + + = ) 1 , , 1 ( - = n i ( 36
)
then from equation [START_REF] Royston | Periodic response of mechanical systems with local nonlinearities using an enhanced Galerkin technique[END_REF], it can be known that, for two consecutive masses, the N th order

GFRFs have the following relationships

= - - - + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Λ Λ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Λ Λ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + + = L l N l L l L N N l L l L N T N l L i N l L i L l N l L l L N N l L l L N T N l L i N l L i N i i N j j j j j Q j Q j j j j j Q j Q 1 1 ) ( , 1 ) ( 1 ) ( , 1 ) ( 1 ) ( , 1 1 1 ) ( , 1 1 1 ) ( , 1 ) ( 1 ) ( , 1 ) ( 1 ) ( , 1 1 ) ( , 1 1 , ) , , ( ) , , ( ) ( ( )) ( ( ) , , ( ) , , ( ) ( ( )) ( ( ) , , ( ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω λ ) 1 , , 1 ( - = n i (37)
Equations ( 22) and [START_REF] Chakraborty | Dynamics of a weakly non-linear periodic chain[END_REF] give a comprehensive description for the relationships between the GFRFs of any two consecutive masses for the nonlinear MDOF system (6).

In addition, denote 0 ) , , ( 1

1 , 0 = N N j j ω ω λ ( N N , , 1 = ) (38) 0 ) , , ( 1 )
, 1 ( = Λ - N N J j j ω ω ( N N , , 1 = ) (39) N N N j j N N J , , 2 if 1 if 0 1 ) , , ( 1 ) ,
( = = ⎩ ⎨ ⎧ = Λ ω ω (40) ) , , ( ) , , ( 1 )
( , 1 ) ( 1 ) , 1 ) ( ( N i L i L N N N i L j j j j ω ω ω ω - - Λ - = Λ ( N N , , 1 = , L i ≤ ≤ 1 ) (41) ) , , ( ) , , ( 1 ) ( , 1 ) ( 1 ) ), ( ( N i
L i L N N N i L j j j j ω ω ω ω - Λ = Λ ( N N , , 1 = , L i ≤ ≤ 1 ) (42)
and

J L = ) 0 (
, then, for the first two masses, from equations ( 16) and ( 33), it can be known that

( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + - + + + + + + + - + + + = = 1 1 1 1 1 , 0 1 2 2 2 1 1 2 1 2 1 ) , 2 ( 1 ) , 1 ( 1 2 , 1 ) ( ) , , ( 1 )
( ) ( ) ( ) , , ( ) , , ( ) , , (
k jc jc k m k jc j j H j j H N N N N N N N N N N N N ω ω ω ω λ ω ω ω ω ω ω ω ω ω ω ω ω λ ( N N , , 1 = ) (43)
Starting with equation ( 43), and iteratively using equations ( 16) and ( 33) from the first mass, it can be deduce that, for the masses that aren't connected to nonlinear components and the Jth spring, the following relationships hold for the GFRFs.

( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + - + - + + + + - + + + = = + - + - + + + + 2 1 1 1 , 1 1 1 1 , 1 1 1 1 1 ) , 1 ( 1 ) , ( 1 1 , ) ( ) , , ( 1 )
( ) , , ( 1 ) ( ) , , ( ) , , ( ) , ,
(

N i i i N i i N N i i N i i N i N i N N i N N i N i i N m k k c c j k jc j j H j j H ω ω ω ω λ ω ω ω ω λ ω ω ω ω ω ω ω ω λ ( n i < ≤ 1 , ) ( , 1 ) ( , l L l L J i - ≠ , L l , , 0 = , N N , , 1 = ) (44)
For the masses that are connected to nonlinear components and the Jth spring, from equations ( 16), ( 27) and ( 28), it can be known that the following relationships hold for the GFRFs.

, 1 ) , 1 ( 1 ) , ( ( 1 ) , 1 ( 1 ) , ( 1 1 1
1 , N N i N N i N i i N i i N N N i N N i N i i N j j H j j jc k j j H j j H ω ω ω ω ω ω ω ω λ ω ω ω ω ω ω λ ( ) ( , 1 ) ( l L l L i - = , L l , , 0 = , N N , , 1 = ) (45) where ( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + - + - + + + + - + + + = + - + - + + + 2 1 1 1 , 1 1 1 1 , 1 1 1 1 1 1 , ) ( ) , , 1
( ) , , ( 1 )
( ) , , ( 1 )
N i i i N i i N N i i N i i N i N i N i i N m k k c c j k jc ω ω ω ω λ ω ω ω ω λ ω ω ω ω λ ( ) ( , 1 ) ( l L l L i - = , L l , , 0 = , N N , , 1 = ) (46) Moreover, denote 1 ) , , (
+ N n n N ω ω λ , ( N N , , 1 = ), 0 1 = + n c and 0 1 = + n k ( 1 , 1 =
. Then, for the last two masses, from equations ( 16) and ( 33) it is can be deduced that

) , , ( ) , , ( ) , , ( 1 ) , , ( 1 ) , 1 ( 1 ) , ( 1 , 1 1 1
, N N n N N n N n n N N n n N j j H j j H ω ω ω ω ω ω λ ω ω λ - - - = = ( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + - + + - + + + - + + + = + + + + ) ( ) , , ( 1 ) , , ( 1 ) ( ) ( 1 1 1 , 1 1 1
, 1 2 1 1 N n n N n n N n n N n n N N n n N n c c j k k m k jc ω ω ω ω λ ω ω λ ω ω ω ω ( N N , , 1 = ) (47)
Starting with equation (47), and iteratively using equations ( 16) and [START_REF] Yuan | On the determination of phase constants for the study of the free vibration of periodic structures[END_REF], it can be deduce that, for the masses that aren't connected to nonlinear components and the Jth spring, the following relationships hold for the GFRFs.

) , , (

N N i N N i N i i N N i i N j j H j j H ω ω ω ω ω ω λ ω ω λ - - - = = ( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + - + + - + + + - + + + = + + + + ) ( ) , , ) , , ( ) , , ( 1 ) , , ( 1 ) , 1 ( 1 ) , ( 1 , 1 1 1 ,
( 1 1 1 , 1 1 1 , 1 2 1 1 N i i N i i N i i N i i N N i i N i c c j k k m k jc ω ω ω ω λ ω ω λ ω ω ω ω ( n i ≤ ≤ 2 , ) ( , 1 ) ( l L l L i - ≠ , L l , , 0 = , N N , , 1 = ) (48) ( 1 ) , , ( 1 ) ( )
For the masses that are connected to nonlinear components and the Jth spring, from equations ( 16), ( 27) and ( 28), it can be known that the following relationships hold for the GFRFs.

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Λ + + + + = = - - - - ) , , ( ) , , ( ) ( 1 1 ) , , ( ) , , ( ) , , ( ) , , ( 1 ) , 1 ( 1 ) , ( 1 1 1 , 1 ) , 1 ( 1 ) , ( 1 1 ,
N N i N N i N i i N i i N N N i N N i N i i N j j H j j jc k j j H j j H ω ω ω ω ω ω ω ω λ ω ω ω ω ω ω λ A c c e p t e d m a n u s c r i p t 12 ( ) ( , 1 ) ( l L l L i - = , L l , , 0 = , N N , , 1 = ) (49) where ( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + + + - + + - + + + - + + + = + + + + - ) ( ) , , ( 1 ) , , ( 1 ) (
) ( ) , , ( 1 1 1 , 1 1 1 , 1 2 1 1 1 1 , N i i N i i N i i N i i N N i i N i N i i N c c j k k m k jc ω ω ω ω λ ω ω λ ω ω ω ω ω ω λ (50)
From different perspectives, equations ( 43)~(46) and equations (47)~(50) give two alternative descriptions for the relationships between the GFRFs of any two consecutive masses for the nonlinear MDOF system (6).

NOFRFs of MDOF Systems with Multiple Nonlinear Components

According to the definition of NOFRF in equation ( 11), the N th order NOFRF of the ith mass can be expressed as

∫ ∏ ∫ ∏ = + + = = + + = = ω ω ω ω ω ω ω ω σ ω σ ω ω ω ω N N N N q q N N q q N N i N i d j U d j U j j H j G ,..., 1 ,..., 1 1 ) , ( ) , ( 1 1 ) ( )
( ) ,..., ( ) ( ) 1 , 1 ( n i N N ≤ ≤ ≤ ≤ (51) where )
(

ω j U is the Fourier transform of ) (t u .
According to equation (44), for the masses that aren't connected to nonlinear components and the Jth spring, equation (51) can be rewritten as

∫ ∏ ∫ ∏ = + + = = + + = + + = ω ω ω ω ω ω ω ω σ ω σ ω ω ω ω ω λ ω N N N N q q N N q q N N i N i i N N i d j U d j U j j H j G ,..., 1 ,..., 1 1 ) , 1 ( 1 1 , ) , ( 1 1 )
( ) ( ) ,..., ( ) , , ( ) ( ( )( ) [ ] ) ( ) ( 1 ) , 1 ( 1 1
, 1 2 1 1 ω ω ω ω λ ω ω j G k jc k jc m k jc N i i i i i i i N i i i + + + - + + + + + - + - + = ( n i < ≤ 1 , ) ( , 1 ) ( l L l L i - ≠ , L l , , 0 = , N N , , 1 = ) (52) Denote ω ω ω ω ω λ ω ω ω λ = + + + + + = = ) ( 1 1 , ) , 1 ( ) , ( 1 , 1 ) , , ( ) ( )
( ) ( N N i i N N i N i i i N j G j G (53)
Therefore, for two consecutive masses that aren't connected to these nonlinear components and the Jth spring, the NOFRFs have the following relationship

( )( ) [ ] 1 1 , 1 2 1 1 1 , ) ( 1 )
(

+ + - + + + + + + - + - + = i i i i i i N i i i i i N k jc k jc m k jc ω ω ω λ ω ω ω λ ( n i < ≤ 1 , ) ( , 1 ) ( l L l L i - ≠ , L l , , 0 = , N N , , 1 = ) (54) where 0 ) ( 1 , 0 = ω λ N .

A c c e p t e d m a n u s c r i p t 13

Similarly, for the masses that are connected to nonlinear components and the Jth spring, from equations ( 44) and ( 45), it can be deduced that

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Γ + + = + + + ) ( ) ( 1 1 ) ( ) ( ) , 1 ( ) , ( 1 , 1 , ω ω ω ω λ ω λ j G j jc k N i N i i i i i N i i N ( ) ( , 1 ) ( l L l L i - = , L l , , 0 = , N N , , 1 = ) (55) where ( ) ( ) [ ] 1 1 , 1 2 1 1 1 , ) ( 1 )
(

+ + - + + + + + + - + - + = i i i i i i N i i i i i N k jc k jc m k jc ω ω ω λ ω ω ω λ (56)
and

∫ ∏ ∫ ∏ = + + = = + + = Λ = Γ ω ω ω ω ω ω ω ω σ ω σ ω ω ω ω N N N N q q N N q q N N i N i d j U d j U j j j ,..., 1 ,..., 1 1 ) , ( ) , ( 1 1 ) ( ) ( ) , , ( ) ( (57)
Equations ( 54)~(57) give a comprehensive description for the relationships between the NOFRFs of two consecutive masses of the nonlinear MDOF system [START_REF] Bartos | Sun, Characterization of Abrupt Nonlinearity by the Volterra-Fourier Method[END_REF].

Using the same procedure, from equations (47)~(50), an alternative description can be established for the following relationships between the NOFRFs of two consecutive masses. For the masses that aren't connected to nonlinear components and the Jth spring

( )( ) [ ] i i i i i i N i i i N i N i i i N i i N k jc k jc m k jc j G j G + + + - + - + = = = + + + - - - ω ω ω λ ω ω ω ω ω λ ω λ 1 1 , 1 2 ) , 1 ( ) , ( , 1 1 , ) ( 1 ) ( ) ( ) ( 1 )
(

( n i ≤ ≤ 2 , ) ( , 1 ) ( l L l L i - ≠ , L l , , 0 = , N N , , 1 =
) (58) For the masses that are connected to nonlinear components and the Jth spring

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Γ + + = = = - - - - - ) ( ) ( 1 1 ) ( ) ( ) ( ) ( 1 ) ( ) , 1 ( ) , ( 1 , ) , 1 ( ) , ( , 1 1
, ω ω ω ω λ ω ω ω λ ω λ j G j jc k j G j G N i N i i i i i N N i N i i i N i i N ( ) ( , 1 ) ( l L l L i - = , L l , , 0 = , N N , , 1 = ) (59) where 1 ) ( , 1 = + ω λ n n N ( N N , , 1
=
) and

( )( ) [ ] i i i i i i N i i i i i N c k jc k m k jc ω ω ω λ ω ω ω λ + + + - + - + = + + + - 1 1 , 1 2 1 , ) ( 1 ) ( (60)
From different perspectives, both equations (54)~(57) and equations (58)~(60) give a comprehensive description for the relationships between the NOFRFs of any two consecutive masses of the nonlinear MDOF system [START_REF] Bartos | Sun, Characterization of Abrupt Nonlinearity by the Volterra-Fourier Method[END_REF]. . Then, from equations (54)~(60), the following important properties of the NOFRFs of MDOF systems with multiple nonlinear components can be obtained.

i) If ) 1 ( L J ≤
, then for the masses (

1 1 - ≤ ≤ J i and n i L L < ≤ ) (
), the following relationships hold.

) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω j G j G j G j G N i N i i i + + = = ( 1 1 - ≤ ≤ J i and n i L L < ≤ ) ( ) (61)
for the masses (

1 ) 1 ( - < ≤ L i J
), the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 2 , 1 ( ) 2 , ( ) 1 , 1 ( ) 1 ,
( ω ω ω ω ω ω j G j G j G j G j G j G N i N i i i i i + + + = = ≠ ( 1 ) 1 ( - < ≤ L i J ) (62) for the masses ( ) ( ) 1 ( L L i L < ≤
), the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 2 , 1 ( ) 2 , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω ω ω j G j G j G j G j G j G N i N i i i i i + + + ≠ ≠ ≠ ( ) ( 1 ) 1 ( L L i L < ≤ - ) (63) ii) If ) ( ) 1 ( L L J L ≤ ≤
, then for the masses ( 1

) 1 ( 1 - < ≤ L i and n i L L < ≤ ) (
) , t h e following relationships hold.

) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω j G j G j G j G N i N i i i + + = = ( 1 ) 1 ( 1 - < ≤ L i and n i L L < ≤ ) ( ) (64)
for the masses ( )

( ) 1 ( L L i L < ≤
), the following relationships hold

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 2 , 1 ( ) 2 , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω ω ω j G j G j G j G j G j G N i N i i i i i + + + ≠ ≠ ≠ ( ) ( 1 ) 1 ( L L i L < ≤ - ) (65) iii) If ) (L L J ≥ , then for the masses ( 1 ) 1 ( 1 - < ≤ L i and n i J < ≤
), the following relationships hold.

) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω j G j G j G j G N i N i i i + + = = ( 1 ) 1 ( 1 - < ≤ L i and n i J < ≤ ) (66)
for the masses (

J i L L < ≤ ) (
), the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 2 , 1 ( ) 2 , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω ω ω j G j G j G j G j G j G N i N i i i i i + + + = = ≠ ( J i L L < ≤ ) ( ) (67)
for the masses (

) ( 1 ) 1 ( L L i L < ≤ -
), the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , 1 ( ) , ( ) 2 , 1 ( ) 2 , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω ω ω j G j G j G j G j G j G N i N i i i i i + + + ≠ ≠ ≠ ( ) ( 1 ) 1 ( L L i L < ≤ - ) (68) iv) For the masses ( 1 ) 1 ) 1 ( , min( 1 - - ≤ ≤ L J i and n i J L L < ≤ ) ), ( max(
), the following relationships of the output frequency responses hold A c c e p t e d m a n u s c r i p t

) ( ) ( ) ( 1 1 , ω ω λ ω j x j x i i i i + + =
15 ( 1 )) 1 ( , min( 1 - ≤ ≤ L J i and n i J L L < ≤ ) ), ( max( ) (69)
where

( ) ( ) [ ] 1 1 , 1 2 1 1 1 , ) ( 1 ) ( + + - + + + + + + - + - + = i i i i i i i i i i i k jc k jc m k jc ω ω ω λ ω ω ω λ (70)
Above fourth properties can be easily extended to a more general case, as the following. v) For any two masses whose positions are

[ ] 1 )) 1 ( , min( , 1 , - ⊆ L J k i , or [ ] 1 )), ( , max( , - ⊆ n L L J k i
, the following relationships hold.

) ( ) ( ) ( ) ( ) ( , ) , ( ) , ( ) 1 , ( ) 1 , ( ω λ ω ω ω ω k i N k N i k i j G j G j G j G = = = ( [ ] 1 )) 1 ( , min( , 1 , - ⊆ L J k i , or [ ] 1 )), ( , max( , - ⊆ n L L J k i ) (71) and ∏ - - = + + + = 1 0 1 , , ) ( ) ( i k d d i d i k i ω λ ω λ (72)
Moreover, the following relationships of their output frequency responses hold

) ( ) ( ) ( , ω ω λ ω j x j x k k i i = (73) vi) If ) 1 ( L J ≤
, for any two masses whose positions are

[ ] 1 ) 1 ( , , - ⊆ L J k i
, the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) ( , ) , ( ) , ( ) 2 , ( ) 2 , ( ) 1 , ( ) 1 , ( ω λ ω ω ω ω ω ω k i N k N i k i k i j G j G j G j G j G j G = = = ≠ ( [ ] 1 ) 1 ( , , - ⊆ L J k i ) (74) vii) If ) (L L J ≥
, then for any masses whose positions are [ ]

J L L k i ), ( , ⊆
, the following relationships hold. ( ,⊆ ) (75) viii) For any two masses whose positions are [ ]

) ( ) ( ) ( ) ( ) ( ) ( ) ( , ) , ( ) , ( ) 2 , ( ) 2 , ( ) 1 , ( ) 1 , ( ω λ ω ω ω ω ω ω k i N k N i k i k i j G j G j G j G j G j G = = = ≠ ( [ ] J L L k i ),
) ( , 1 ) 1 ( , L L L k i - ⊆
, the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , ( ) , ( ) 2 , ( ) 2 , ( ) 1 , ( ) 1 , ( ω ω ω ω ω ω j G j G j G j G j G j G N k N i k i k i ≠ ≠ ≠ (76)
ix) For any two masses whose positions are

[ ] 1 ) 1 ( , 1 - ⊆ L i and [ ] n L k ), 1 ( ⊆ or [ ] 1 ) ( , 1 - ⊆ L L i and [ ] n L L k
), ( ⊆ , the following relationships hold.

) ( ) ( ) ( ) ( ) ( ) ( ) , ( ) , ( ) 2 , ( ) 2 , ( ) 1 , ( ) 1 , ( ω ω
ω ω ω ω j G j G j G j G j G j G N k N i k i k i ≠ ≠ ≠ (77)
The nine relationships between NOFRFs of nonlinear MDOF systems reveal, for the first time, very important properties of the nonlinear systems. They reveal how the linear

A c c e p t e d m a n u s c r i p t

16 system parameters govern the propagation of the nonlinear effect induced by the nonlinear components in the whole system. These properties can be applied to detect and locate the nonlinear components in engineering structures which often is related to faults.

Numerical Study and Discussions

Case Studies

To verify above analysis results, a damped 10-DOF oscillator was adopted, in which the fourth and sixth springs were nonlinear. The damping was assumed to be proportional damping, e.g.,

K C µ = .
The values of the system parameters are 1

10 1 = = = m m , , 10 6 . 3
× = = = = k k k , 8 . 0 1 8 7 6 k k k k × = = = , 9 . 0 1 9 k k × = µ =0.01, 0 ) 3 , 6 ( ) 2 , 6 ( ) 3 , 4 ( ) 2 , 4 ( =
= = = w w w w 2 1 ) 2 , 4 ( 8 . 0 k r × = , 3 1 ) 3 , 4 ( 4
. 0 k r × = , ) 2 , 4 ( ) 2 , 6 ( 5
. 0 r r × = , ) 3 , 4 ( ) 3 , 6 ( 1
. 0 r r × = and the input is a harmonic force, ) 20 2 sin( ) ( t A t u × = π .
When a nonlinear system is subject to a harmonic input, )

= t A t u F (78) it can be derived that the output spectrum ) ( ω j Y cos( ) ( β ω +
of nonlinear systems can be simply expressed as [START_REF] Lang | Evaluation of output frequency responses of nonlinear systems under multiple inputs[END_REF] [ ]

∑ + - = - + - + = 2 / ) 1 ( 1 ) 1 ( 2 ) 1 ( 2 ) ( ) ( ) ( k N n F n k F n k F jk A jk G jk Y ω ω ω ( N k , , 1 , 0 = ) (79)
where [] . means to take the integer part, and

β ω ) 2 ( | | )! ( ! ! 2 1 ) ) 2 ( ( k n j n n F n e A k n k n k n j A + - - = + - (80) ) , , , , , ( ) ) 2 ( ( k
n F F k F F n F n j j j j H k n j G - - - = + - ω ω ω ω ω (81)
According to equations (79)~(81), if only the NOFRFs up to the 4 th order is considered, the frequency components of the outputs of the 10 masses 82), it can be seen that, using the method in [START_REF] Lang | Energy transfer properties of nonlinear systems in the frequency domain[END_REF], two different inputs with the same waveform but different strengths are sufficient to estimate the NOFRFs up to 4 th order. Therefore, in the numerical studies, two different inputs were A=0.8 and A=1.0 respectively. The simulation studies were conducted using a fourth-order Runge-Kutta method to obtain the forced response of the system. In this case, the position of the input force is on the right of the two nonlinear components. The evaluated results of

) ( ) ( ) ( ) ( ) ( 3 ) 3 , ( 1 ) 1 , ( F F i F F i F i j A j G j A j G j x ω ω ω ω ω + = ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 4 ) 4 , ( 2 ) 2 ,
( F F i F F i F i j A j G j A j G j x ω ω ω ω ω + = ) 3 ( ) 3 ( ) 3 ( 3 ) 3 , ( F F i F i j A j G j x ω ω ω = ) 4 ( ) 4 ( ) 4 ( 4 ) 4 , ( F F i F i j A j G j x ω ω ω = ) 10 , , 1 ( = i (82) From equation (
) ( 1 F j G ω , ) ( 3 F j G ω , ) 2 ( 2 F j G ω and ) 2 ( 4 F j G
ω for all masses are given in Table 1. According to Property iii) in the previous section, it can be known that the following relationships should be tenable. ) In this case, the input force is located between the two nonlinear components. The evaluated results of

i i F H i F i F H i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 9 , 8 ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F
i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = ≠ = for 7 , 6 , 2 , 1 = i ) ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F
i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 9 , 8 , 5 , 4 , 3 = i ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , 4 ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( 1 , 2 F
i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = ≠ = for 5 , 4 , 3 = i (83) Table 1, the evaluated results of ) ( 1 F j G ω , ) , 7 , 6 , 2 , 1 = i ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , 4 ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( 1 , 2 F
G ω (×10 -9 ) ) 2 ( 2 F j G ω (×10 -8 ) ) 2 ( 4 F j G ω (×10 - ( 3 F j G ω , ) 2 ( 2 F j G ω and ) 2 ( 4 F j G ω ) ( 1 F j G ω (×10 -6 ) ) ( 3 F j
) ( 1 F j G ω , ) ( 3 F j G ω , ) 2 ( 2 F j G ω and ) 2 ( 4 F j G
ω for all masses are given in Table 3. According to Property ii) in the previous section, it can be shown that the following relationships should be tenable. From the NOFRFs in Table 3,

) ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 9 , 8 , 7 , 6 , 2 , 1
= i ) ( ) ( ) ( ) ( ) ( ) ( 1 , 3 ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( 1 , 1 F i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = ≠ = for 5 , 4 , 3 = i ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , 4 ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( 1 , 2 F
i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = = = for 9 , 8 , 7 , 6 , 2 , 1
= i ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( 1 , 4 ) 4 , 1 ( ) 4 , ( ) 2 , 1 ( ) 2 , ( 1 , 2 F i i F i F i F i F i F i i j j G j G j G j G j ω λ ω ω ω ω ω λ + + + + = ≠ = for 5 , 4 , 3 = i (84) Table 3, the evaluated results of ) ( 1 F j G ω , ) ( 3 F j G ω , ) 2 ( 2 F j G ω and ) 2 ( 4 F j G ω ) ( 1 F j G ω (×10 -6 ) ) ( 3 F j G ω (×10 -8 ) ) 2 ( 2 F j G ω (×10 -8 ) ) 2 ( 4 F j G ω (×10 -
) ( 1 , 1 F i i jω λ + , ) ( 1 , 3 F i i jω λ + , ) 2 ( 1 , 2 F i i j ω λ + and ) 2 ( 1 , 4 F i i j ω λ + ( 9 , , 1 = i
) can be calculated. The results are given in Table 4. It can be seen that the results shown in Tables 4 have a strict accordance with the relationships in (84). The two numerical case studies verify the properties of the NOFRFs of MDOF systems with multiple nonlinear components derived in the present study. These properties can provide a convenient method to detect the positions of the nonlinear components in a MDOF system by analyzing the relationships between the NOFRFs.

F i i jω λ + , ) ( 1 , 3 F i i jω λ + , ) 2 ( 1 , 2 F i i j ω λ + and ) 2 ( 1 , 4 F i i j ω λ + ) ( 1 , 1 F i i jω λ + ) ( 1 , 3 F i i jω λ + ) 2 ( 1 , 2 F i i j ω λ + ) 2 ( 1 , 4 F i i j ω λ + i=1 0.

Discussion

In engineering practice, a wide class of real life structures can be modeled as periodicstructures-like MDOF systems, which are defined as structures consisting of identical substructures connected to each other in identical manner, such as periodically supported beams [START_REF] Duhame | Finite element analysis of the vibrations of waveguides and periodic structures[END_REF]~ [START_REF] Yuan | On the determination of phase constants for the study of the free vibration of periodic structures[END_REF] and plates [32][33]. If one or more components are of nonlinear properties, then the systems can behave nonlinearly. Efforts have been made to study the [START_REF] Marathe | Wave attenuation in nonlinear periodic structures using harmonic balance and multiple scales[END_REF] using one-dimensional and multidimensional MDOF models. In addition, the detection and location of faults and defects in periodic structures and machines are also interesting to many researchers. Zhu and Wu [START_REF] Zhu | The Characteristic receptance method for damage detection in large Mono-coupled Periodic structures[END_REF] have studied the detection of damages in large periodic structures. In their studies, the periodic structure with damage is still considered to be a linear system, and the location and magnitude of damage in large mono-coupled periodic systems were estimated using measured changes in the natural frequencies. However, in engineering practices, the local faults and defects can often make the structures and machines behave nonlinearly. In the latter cases, obviously the properties discovered in this study can provide a convenient way to detect the positions of the faults. For example, consider the 10-element periodic mass-spring system shown in Figure 2, which is a specific form of the mass-spring system investigated in [START_REF] Zhu | The Characteristic receptance method for damage detection in large Mono-coupled Periodic structures[END_REF], and assume the sixth spring is damaged and of nonlinear property. To detect the position of the damaged spring, two excitations can be used to excite the system, and then the NOFRFs up to the fourth order of all masses can be estimated from the responses. Obviously, if the excitation force is acting on the 10 th mass, then according to the property iii), there exist the following relationships:

) (

ω j G j G j G j G i i i i + + = (i =1, 2, 3, 4,) ) ( ) ( ) ( ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( ω ω ω
( ) ( ) ( ) ( ) 3 , 1 ( ) 3 , ( ) 1 , 1 ( ) 1 , ( ω ω ω ω j G j G j G j G i i i i + + ≠ (i =5, 6, 7, 8, 9) (86) (85) )
Clearly, the relationship (86) provides a direct way to detect the position of the damaged spring.

Figure 2, The 10-element periodic mass-spring system whose 6 th spring is fault

The locally nonlinear MDOF (6) can also be used to describe the transversal motion of the tall apartment block with damages shown in Figure 3. Based the same model, Sakellariou and Fassois [START_REF] Fassois | Stochastic output error vibration-based damage detection and assessment in structures under earthquake excitation[END_REF][41] have used a stochastic output error vibration-based methodology to detect the damage in structures where the damage elements were modeled as components of cubic stiffness. No doubt the properties obtained in this study can also provide a convenient way to detect the position of the damage in this kind of structures. In addition, for simplicity, in this study the one-dimensional nonlinear MDOF system, each mass of which has only one freedom degree, is adopted. Nevertheless, following the same procedure used in this study, the nine properties can be extended to the multidimensional nonlinear MDOF systems by simply replacing the scalar forms of m i , c i and k i as matrix forms, consequently, the

i i N + ( 1 , , 0 - = n i ; N N , , ) ( 1 , ω λ j
) also have a matrix form. The multi-dimensional cases could be more complicated than the onedimensional cases, but the NOFRF properties of multi-dimensional nonlinear MDOF systems can be used to analyze a more wide class of structures. For example, the beam shown in Figure 4(a) can be represented as a series of rigid blocks connected by rotational and transverse springs shown in Figure 4(b). The rotational spring approximates the bending of the beam and the transverse spring approximates the shear. Neild, Mcfadden and Williams [START_REF] Neild | A discrete model of vibrating beam using time-stepping approach[END_REF] have used this model to analyze the beam with one breathing crack. When cracks present in this beam, the whole beam can behave nonlinearly. Clearly, this beam is a two-dimensional nonlinear MDOF system, and the detection of the crack position in this beam could be easily achieved using the properties of the NOFRFs for multi-dimensional MDOF systems.

m n J n m 1 J 1 m i J i (b) y 1 θ 1 y i θ i y n θ n (a)
x 1

x n-1 Moreover, as the response spectrum of each mass of nonlinear MDOF systems is totally determined by the associated NOFRFs and the input, based on the relationships between the NOFRFs, new spectra analysis methods can be developed to detect the position of nonlinear elements in the systems, which may only involve a very simple procedure. We are currently working on this and the results will be presented in a future publication.

From equation (70) it can be seen that the ratios between the NOFRFs of two consecutive masses are mainly determined by the linear parameters. Therefore equation (70) reveals how the linear system parameters govern the relationships between the NOFRFs of two consecutive masses. This fact provides a convenient way to estimate these linear parameters, which has been elaborated in [START_REF] Peng | Linear Parameter Estimation for Multi-Degree-of-Freedom Nonlinear Systems Using Nonlinear Output Frequency Response Functions[END_REF].

Conclusions

In this paper, significant relationships between the NOFRFs of MDOF systems with multiple nonlinear components have been derived and verified by numerical studies. The results reveal, for the first time, important properties of this general class of MDOF nonlinear systems. The potential of using these properties to detect and locate faults in engineering structures is also discussed.

It is worth noting here that, theoretically, the obtained relationships about NOFRFs are valid for nonlinear systems whose responses can be described using the Volterra series, which covers a considerably wide range of operating conditions of nonlinear systems [START_REF] Schetzen | The Volterra and Wiener theories of nonlinear systems[END_REF]~ [START_REF] Bartos | Sun, Characterization of Abrupt Nonlinearity by the Volterra-Fourier Method[END_REF][44] [START_REF] Chatterjee | Convergence analysis of Volterra series response of nonlinear systems subjected to harmonic excitation[END_REF]. In practice, because the validity of the Volterra series is dependent on the amplitude of the external force input which is normally controllable during fault detection oriented structural tests, provided that the amplitude of testing input is properly selected, the important relationships between the NOFRFs will hold and can therefore be used for structural fault diagnosis. Moreover, for convenience, in this study the onedimensional MDOF system is adopted, that is each mass has only one degree of freedom, but the obtained results can be easily extended to the multi-dimensional MDOF cases where each mass in the MDOF system has more than one degree of freedom. In the analysis, if m i , c i and k i are taken as matrix forms and ) , , 1 ( , n i x i = are taken as vectors, the same results can be achieved. ), property (63) can be proved.

Proof of Property (ii) and (iii)

Following the same procedure used for proof of Property (i), the second and third properties can be proved. The details are omitted here.

Proof of Property (iv)

The fourth property is also straightforward since, according to equation [START_REF] Stegmayer | Volterra series and neural networks to model an electronic device nonlinear behavior[END_REF], the output frequency response of the ith mass can be expressed as , then the fourth property is proved.

∑ = + + = N k k k i i j U j G j x 1 ) , 1 ( 1 )

Proof of Properties (v)~(ix)

The proof of the properties (v)~(ix) only needs some simple calculations. The details are therefore omitted here.

+

  …,P)the nonlinearity related parameters of the L(i) th nonlinear component displacement and the output frequency response of the i th mass ) the ith row and the lth column of the inverse matrix j th order Volterra kernel associated to the i th mass ) order NOFRF associated to the i th mass ) the n th GFRFs of the i th and (i+1) th masses ) the ratio between the n th NOFRFs of the i th and (i+1) th masses

A c c e p t e d m a n u s c r i p t 5 Equations ( 3 )

 53 ~(

A c c e p t e d m a n u s c r i p t 6

 6

A c c e p t e d m a n u s c r i p t 7 depends

 7 not only on n H (n=1,…,N) but also on the input ) ( ω j U

  …,n; l=1,…,n) is the element of the ith row and the lth column of the inverse matrix of )

A c c e p t e d m a n u s c r i p t 11

 11

A c c e p t e d m a n u s c r i p t 14 4. 3

 143 The Properties of NOFRFs of Locally Nonlinear SystemsWithout loss of generality, assume )

A c c e p t e d m a n u s c r i p t 17 Case 1 .

 171 Input Force Acting on the Eighth Mass (J = 8)

A c c e p t e d m a n u s c r i p t 18

 18

A c c e p t e d m a n u s c r i p t 20 dynamics

 20 of the nonlinear structures[START_REF] Chi | On damped Non-linear Dynamic systems with many degrees of freedom[END_REF]-[

9 A c c e p t e d m a n u s c r i p t 21 Figure 3 ,

 9213 Figure 3, Schematic diagram of the n-storey building block

1 A c c e p t e d m a n u s c r i p t 22 Figure 4 ,

 1224 Figure 4, Modelling of a beam as discrete blocks

  the first three properties, it can be known that,

A c c e p t e d m a n u s c r i p t 25

 25

Table 2 .

 2 It can be seen that the results shown in Table2have a strict accordance with the relationships in (83).

	10 )

Table 2 ,

 2 the evaluated values of ) (

							1 λ + , 1 i i	jω	F	,	1 λ + , 3 i i	(	jω	F	)	,	1 λ + , 2 i i	(	2 j ω	F	)	and	1 λ + , 4 i i	(	2 j ω	F	)
		1 λ + , 1 i i	(	jω	F	)	1 λ + , 3 i i	(	jω	F	)			1 λ + , 2 i i	(	2 j ω	F	)	1 λ + , 4 i i	(	2 j ω	F	)
	i=1	0.5391 -0.0630i	0.5391-0.0630i		0.5085 -0.1741i	0.5085-0.1741i
	i=2
	Case 2.

Input Force Acting on the Fifth Mass (J = 5)

Table 4

 4

	, the evaluated values of	, 1	1	(	)

Acknowledgements

The authors gratefully acknowledge the support of the Engineering and Physical Science Research Council, UK, for this work.

0.7407-0.1588i 0.7407-0.1588i 0.5766-0.3580i 0.5766-0.3580i i=3 0.8219-0.2811i 0.9888-0.4095i 0.2378-1.1151i 0.1388-0.5139i i=4 0.7994-0.3924i 0.7415-0.3778i 0.4796-0.1604i 0.6244-0.3135i i=5 0.6285-0.4587i -0.8354+0.1958i -1.1257-0.1490i -0.9941+0.0731i i=6 0.6197-0.3882i 1.3389+0.8002i 1.3194+1.0981i 1.3194+1.0981i i=7 0.6513-0.3656i 0.9242+1.0986i 1.4455+1.3074i 1.4455+1.3074i i=8 0.6666+0.5608i 0.6666+0.5608i 0.9473+1.3813i 0.9473+1.3813i i=9 0.8299+0.2137i 0.8299+0.2137i 0.7602+0.6028i 0.7602+0.6028i

Appendix 1: Proof of Property (i)

The first property is straightforward. For the masses on the left of the Jth mass, substituting 0 ) (

) into equation (54), it is obtained that ( ) ( )

Subsequently, substituting (A-1) into equation (55) yields

Iteratively using above procedure until i=(J-1), for the masses (

), property (61) can be proved.

Similarly, substituting 1 ) (

Subsequently, substituting (A-3) into equation ( 58), it can be deduced that

Iteratively using above procedure until ) (L L i = , for the masses (

), property (61) can be proved.

Obviously, from equations ( 40) and (57), it is known that

and equation (A-5) into (54), it can be deduced that

)

A c c e p t e d m a n u s c r i p t 24

Iteratively using this procedure until

), property (62) can be proved.