Jun Jiang
email: jun.jiang@mail.xjtu.edu.cn

Determination of the global responses characteristics of a piecewise smooth dynamical system with contact

Keywords: Piecewise smooth systems, Global response characteristics, Stability analysis

In this paper the global response characteristics of a piecewise smooth dynamical system with contact, which is specifically used to describe the rotor/stator rubbing systems, is studied analytically. A method to derive the global response characteristics of the model is proposed by studying each piece of the equations corresponding to different phases of the rotor motion, i.e., the phase without rubbing, the phase with rubbing and the phase of self-excited backward whirl. After solving the typical responses in each phase and deriving the corresponding existence boundaries in the parameter space, an overall picture of the global response characteristics of the model is obtained. As is shown, five types of the coexistences of the different rotor responses and deep insights into the interactive effect of parameters on the dynamic behavior of the model are gained.

Backward whirl "natural" frequency of nonlinear coupled rotor/stator system ω 2

Natural frequency of the rotor system with zero clearance

1 ω b
Normalized whirl frequency of the rotor, ω b = ω w /ω 2 ω w Whirl frequency of the rotor Ω Normalized rotating speed of the rotor, ω/ω 2 Ω e Critical rotating speed for the start of dry friction backward whirl ζ Damping ratio of the rotor system

Introduction

Rotor-to-stator rubbing is a malfunction in rotating machinery that degrades the machine performance and may lead to a catastrophic failure of a whole machine through dry whip in the worst case. A large number of research papers deal with the dynamics of rub-related phenomena in rotor/stator contact systems. Some papers studied the possible rubbing responses induced through the rotor/stator rubbing and tried to derive the relationship between the rubbing responses and the system parameters, as summarized in [START_REF] Muszynska | Rotor-to-stationary element rub-related vibration phenomena in rotating machinery-literature survey[END_REF]. Through the numerical and experimental studies, many rubbing behavior in rotor/stator contact systems have now become well known, for instance, the jump phenomenon and the synchronous full annular rubs [START_REF] Black | Interaction of a whirling rotor with a vibrating stator across a clearance annulus[END_REF][START_REF] Yu | Rotor/seal experimental and analytical study on full annular rub[END_REF][START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[END_REF], the partial rubs in sub-and super-synchronous whirl [START_REF] Bently | Forced subrotative speed dynamic action of rotating machinery[END_REF][START_REF] Childs | Rub-induced parametric excitation in rotors[END_REF][START_REF] Ehrich | High order subharmonic response of high speed rotors in bearing clearance[END_REF], the partial rubs in quasi-periodic whirl [START_REF] Day | Asymptotic expansions in nonlinear rotor dynamics[END_REF][START_REF] Kim | Quasi-periodic response and stability analysis for a nonlinear Jeffcott rotor[END_REF], the chaotic motion [START_REF] Goldman | Chaotic behavior of rotor/stator system with rubs[END_REF][START_REF] Karpenko | Regular and chaotic dynamics of a discontinuously nonlinear rotor system[END_REF] as well as the dry friction backward whirl which includes two regimes: of dry whirl and dry whip [START_REF] Zhang | Dynamic instability of multi-degree-of-freedom flexible rotor systems due to full annular rub[END_REF][START_REF] Crandall | From whirl to whip in rotordynamics[END_REF][START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF][START_REF] Jiang | The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems[END_REF]. The numerical and experimental studies are useful in detecting the possible rubbing responses in rotor/stator contact systems and in capturing some trends of dependence of the responses on system parameters. But a huge amount of work is inevitable in these ways in order to determine the boundaries of different rubbing responses of the rotor/stator contact systems. Furthermore, it was usually found that different researchers might present quite different response cascades and the inherent reason was not well explained. For instance, it was observed in tests of [START_REF] Choi | Investigation on the whirling motion of full annular rotor rub[END_REF][START_REF] Lingener | Experimental investigation of reverse whirl of a flexible rotor[END_REF] that the rotor under the unbalanced excitation would follow a normal behavior with the increase or decrease of the rotating speed, that is, the rotor follows the response cascade: no-rub motion → a synchronous full annular rub → no-rub motion, and vise versa. Jump phenomenon was detected in both run-up and rundown process. However, it was reported that when an additional outside disturbance (hitting on the shaft) was applied on the rotor during the process, even at a very low rotor speed, the dry friction backward whirl with super-synchronous whirl frequency would be triggered from a synchronous response. A quite different response scenario with the increase of the rotor speed was observed in the tests of [START_REF] Yu | Rotor/seal experimental and analytical study on full annular rub[END_REF][START_REF] Crandall | From whirl to whip in rotordynamics[END_REF] as: no-rub motion → a synchronous full annular rub → a partial rub → dry whip (heavy rub). So dry whip was triggered solely by the mass unbalance for the rotor/stator systems in this case, which is different from the case mentioned above where an outside disturbance was needed. This indicates the necessity to further explore the existence boundaries of different responses of rotor/stator contact systems so as to get an overall picture on the interactive effect of system parameters on the appearance of different response cascades. The recent analytical works [START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[END_REF][START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF][START_REF] Jiang | The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems[END_REF] have provided some satisfactory explanation on the response characteristics of rotor/stator contact systems mentioned in above tests.

This paper aims to depict the global response characteristics of a specific rotor/stator contact system by determining the existence conditions of the most typical rotor responses of the piecewise smooth nonlinear system. As known, the piecewise smooth nonlinear model consists of three pieces of equations of motion governing three phases of rotor motion, namely, the phase without rub, the phase with rubs and the phase in self-excited dry friction backward whirl. As will be shown in this paper, the existence condition for the most often observed responses, i.e., synchronous no-rub motion, synchronous full annular rub, partial rub and dry friction backward whirl, can be derived solely from a single piece of the equations even though some responses, i.e., the partial rubs, run across two phases of the rotor motion. After incorporating the existence boundaries of all of the responses together into the parameter space, an overall picture of the response characteristics of the rotor/stator contact system is achieved. Deep insights into the interactive effect of different parameters on the system responses as well as the coexistence behavior are gained. Moreover, the analytical methods used in this paper can be also applied to the more complex rotor/stator systems, i.e., the model considering the stator dynamics and the contact stiffness at contact surfaces.

The paper is organized in the following ways. In Sect. 2, the equations of motion for each phase of motion of the rotor/stator contact system are first introduced. Then the equations for different phases of motion are studied individually in Sect. 3 in order to derive the solutions and the existence boundaries of the corresponding responses. Section 4 will illustrate the global response characteristics of the present rotor/stator contact system in the parameter space through examples. Finally in Sect. 5, the conclusions from this work are drawn.

Mathematical model

The rotor/stator contact system consists of a Jeffcott rotor and a stator that is modeled as an added stiffness (see Fig. 1). This model physically describes the rotor in contact with a non-rotating susceptible circumferential stator or a mechanical seal. The rotor has a weightless shaft carrying a disk with mass m at the middle of the span. The mass center of the rotor is located at a distance e from its geometrical center. The stiffness of the shaft is k s . The clearance between the rotor and the stator is denoted by r 0 . The stator, which is rigidly fixed, has an elastic contact surface modeled as a symmetrical set of radial springs with isotropic The whirl frequency of the shaft and the rotor radius at the axial position of contact are denoted by ω w and r disk respectively. In the present analysis the gravity is neglected.

The governing equations of the rotor/stator contact system can be written as

m ẍ + c ẋ + k s x + Θk b 1 - r 0 r x -μsign(v rel )y = meω 2 cos ωt m ÿ + c ẏ + k s y + Θk b 1 - r 0 r μsign(v rel )x + y = meω 2 sin ωt v rel = ωr disk + ω w r ( 1
)
where r = x 2 + y 2 and v rel represents the relative velocity between the rotor and the stator at the contact point. Θ is Heaviside function with Θ = 0 when r < r 0 and Θ = 1 when r ≥ r 0 . Equation (1) implies that only the sliding contact is assumed to occur between the rotor and the stator.

As can be seen from ( 1) the equations of motion, in fact, consist of three pieces that govern three phases of rotor motion: (i) the phase without no-rub when Θ = 0; (ii) the phase with rubs Θ = 1 (but v rel > 0 all the time); (iii) the phase of self-excited backward whirl (Θ = 1 and v rel fluctuates around zero where the dry friction effect must be taken into account). The first two phases of motion are quite familiar and widely used in the previous works to investigate the dynamics of the rotor/stator rubbing. The last phase of motion exhibits the well-known rubbing response: dry friction backward whirl, which contains two regimes-of dry whirl and dry whipas defined in [START_REF] Crandall | From whirl to whip in rotordynamics[END_REF] and which is not really well explored.

Equation (1) will be first rewritten in the nondimensional form:

X + 2ζ X + βX + Θ 1 - R 0 R X -μsign(V rel )Y = Ω 2 cos Ωτ Y + 2ζ Y + βY + Θ 1 - R 0 R μsign(V rel )X + Y = Ω 2 sin Ωτ V rel = R disk Ω + Rω b ( 2
)
where the non-dimensional variables are defined as:

X = x e , Y = y e , R 0 = r 0 e R disk = r disk e , V rel = v rel e , 2ζ = c √ k b m β = k s k b , ω 2 = k b m , Ω= ω ω 2 , ω b = ω w ω 2 with R = √ X 2 + Y 2 and Θ = 0 when R < R 0 and Θ = 1 when R ≥ R 0 .
The prime indicates the differentiation with respect to the dimensionless time variable τ = ω 2 t .

By defining

W = X + jY and |W | = √ X 2 + Y 2 , ( 2
) can be written in a complex form. Rather than writing a unified complex equation, we will write the equations into three pieces of equations according to the phases of motion.

The piece of equation governing the phase without rubbing

(|W | < R 0 ) is given by W + 2ζ W + βW = Ω 2 e jΩτ ( 3
)
The piece of equation describing the phase with rubbing (|W | ≥ R 0 ) can be written as

W + 2ζ W + βW + (1 + jμ) 1 - R 0 |W | W = Ω 2 e jΩτ ( 4
)
The piece of equations describing the phase of selfexcited backward whirl (|W | ≥ R c with the sign of V rel alternating) can be written as

W + 2ζ W + βW + 1 + jμsign(V rel ) × 1 - R 0 |W | W = Ω 2 e jΩτ V rel = R disk Ω + |W |ω b (5)
As is seen from the above discussion, the dynamics of the rotor/stator system in different phases of motion are governed by different equations of motion depending upon the amplitudes of the rotor deflection.

The equation in each phase of motion is smooth. So we can solve each piece of equation(s) to get the corresponding solution and then carry out the stability analysis to determine its existence boundary. As will be shown, the most often observed responses in the rotor/stator systems, such as the synchronous no-rub motion, the synchronous full annular rub motion and the dry friction backward whirl, can all be determined in this way. Although the quasi-periodic partial rub motion alternates between the phase without rubbing and the phase with rubbing, its existence boundaries can still be determined by investigating the piece of equation governing the phase with rubbing. After the existence boundaries of the solutions from all the pieces of equations are determined, the global response behavior of the system in the parameter space can be depicted.

3 The solutions and stability analysis

The solution in the phase without rubbing

In this phase, the rotor is free from the contact with the stator, the synchronous no-rub motion is the only response and can be derived from ( 3). The periodic solution can be written as

W 0 = Ae j (Ωτ +ϕ) with A = Ω 2 (β -Ω 2 ) 2 + (2ζ Ω) 2 , tgϕ = -2ζ Ω β -Ω 2 (6)
The periodic solution is unconditionally stable. However, the solution is physically valid only when its amplitude is smaller than the clearance, namely A < R 0 , due to the existence of stator. To find the parameter regions where the synchronous no-rub motion exists, let A = R 0 in (6) and we get the parameter equation defining the boundary of no-rub motion:

R 2 0 -1 Ω 4 + R 2 0 4ζ 2 -2β Ω 2 + R 2 0 β 2 = 0 (7)
Equation ( 7) yields two pairs of real roots, each with an equal magnitude and opposite sign. Since the rotating speed is supposed to be positive and the synchronous no-rub response always undertakes forward whirl, two positive roots, Ω l and Ω u , are chosen to define the boundaries of no-rub motion. That is, the synchronous no-rub motion exists when the rotating speed Ω is less than Ω l or greater than Ω u .

The solutions in the phase with rubbing

Synchronous full annular rub motion is a response by which the rotor is in the full annular contact with the stator and whirls forward with constant amplitude at a frequency equal to the rotating speed of the rotor. Equation ( 4) can be used to determine this response by supposing a proper form of the solution according to the characteristics of the synchronous full annular rub.

Let us denote the periodic solution as W 1 = Be j (Ωτ +ϕ) , two nonlinear equations of B and ϕ are yielded after substituting W 1 into (4) and separating the real and imaginary parts of the obtained equation,

β -Ω 2 B + (B -R 0 ) = Ω 2 cos ϕ 2ζ ΩB + μ(B -R 0 ) = -Ω 2 sin ϕ (8)
To eliminate ϕ from the above equations, we square both equations of ( 8) and then sum them up to get a second order polynomial equation of B in the form

(2ζ Ω + μ) 2 + 1 + β -Ω 2 2 B 2 -2R 0 1 + β -Ω 2 + μ(2ζ Ω + μ) B + R 2 0 1 + μ 2 -Ω 4 = 0 ( 9 )
One can solve for the positive real roots of B from ( 9) and then substitute them back into either of the two equations of (8) to get ϕ for the given set of parameters. From the condition with rubbing, B ought to be greater than R 0 .

The response of the synchronous full annular rub corresponds to the stable periodic solution of (4). So the stability of the solution must be analyzed. To carry out the stability analysis on the periodic solution of the nonlinear system, the physical understanding of the response is of good use. It is known that a periodic solution in the inertial coordinate system will correspond to a fixed point in the rotating coordinate system at the same rotating speed as the period of the solution. To transform (4) into the rotating coordinate system, we use the transformation W 1 (τ ) = ρ r (τ )e jΩτ . Substituting it into (4) and making some manipulations, yields

ρ r + 2(ζ + jΩ)ρ r + 1 + β -Ω 2 + j (2Ωζ + μ) ρ r -(1 + jμ)R 0 ρ r |ρ r | = Ω 2 (10)
Equation ( 10) has a fixed point in the form of ρ r0 = Be jϕ with B and ϕ being defined by [START_REF] Day | Asymptotic expansions in nonlinear rotor dynamics[END_REF].

Further, the stability analysis on the fixed point, and thus on the stability of the synchronous full annular rub solution in the inertial coordinate system, will be carried out. To accomplish this task, the complex amplitude of the response is written in the component form ρ r = X r + jY r , where the subscript "r" indicates the variable in the rotating coordinate system.

By defining the state vector in the rotating coordinate system as U = {X r Y r X r Y r } T , (10) can be written in a time-invariant differential equation

U = AU + g( U ) = G( U ) ( 11
)
where A and g( U ) are respectively the coefficient matrix and the nonlinear vector term:

A = ⎡ ⎢ ⎢ ⎣ 0 0 1 0 0 0 0 1 -(1 + β -Ω 2 ) 2Ωζ + μ -2ζ 2Ω -(2Ωζ + μ) -(1 + β -Ω 2 ) -2Ω -2ζ ⎤ ⎥ ⎥ ⎦ and g( U ) = {0 0 H X + Ω 2 H Y } T with H X = R 0 R (X r -μY r ), H Y = R 0 R (μX r + Y r )
. Linearizing [START_REF] Karpenko | Regular and chaotic dynamics of a discontinuously nonlinear rotor system[END_REF] about the solution of the fixed point U 0 = {B cos ϕ B sin ϕ 0 0}, yields a linear, timeinvariant differential equation δU = J U 0 δU [START_REF] Choi | Investigation on the whirling motion of full annular rotor rub[END_REF] where δU = U -U 0 is the perturbation solution to the fixed point. J( U 0 ) is the so-called Jacobian matrix which is given in the following form:

J U 0 = ⎡ ⎢ ⎣ 0 0 1 0 0 0 0 1 -1 -β + Ω 2 (1 -a)μ + 2ζ Ω -2ζ 2Ω -μ -2ζ Ω -1 -β + Ω 2 + a -2Ω -2ζ ⎤ ⎥ ⎦ with a = R 0 /B ( 13
)
From Jacobian matrix given in [START_REF] Lingener | Experimental investigation of reverse whirl of a flexible rotor[END_REF], one can set up the characteristic equation of the solution,

s 4 + b 3 s 3 + b 2 s 2 + b 1 s + b 0 = 0 ( 14
) with b 3 = 4ζ, b 2 = 2Ω 2 + 4ζ 2 -a + 2(1 + β), b 1 = 2μ(2 -a)Ω -2ζ a + 4ζ(1 + β + Ω 2 ), b 0 = Ω 4 + (4ζ 2 + a -2)Ω 2 + 2μζ (2 -a)Ω + (1 + β) 2 - a(1 + β) -μa + μ 2 .
The real parts of the characteristic roots decide the stability of the fixed point, and thus the stability of the synchronous full annular rub solution. Since the stability boundaries are of interest, bifurcation theory is being used as is done in [START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[END_REF]. The saddle-node bifurcation condition, by which a zero characteristic root exists, defines the existence boundary of the solution. The Hopf bifurcation condition, whereby a pair of pure imaginary characteristic roots appears, gives the boundary between the quasi-periodic partial rub and the synchronous full annular rub.

Quasi-periodic partial rubbing is a response of the rotor/stator system, by which the rotor rubs the stator intermittently. There are generally more than one frequency component in the system responses during partial rubs. The equations governing the partial rubs should be a combination of (3) for the phase without rubbing and (4) for the phase with rubbing. Generally, we have difficulty in obtaining the analytical solutions of the piecewise smooth equations. Without the solution of partial rub, it seems impossible to determine their existence boundaries.

Fortunately, we can still sketch the existence boundaries of quasi-periodic partial rub based on the study of other responses and on the understanding of the physical reason for the loss of stability of the partial rub. For instance, it is easy to determine the boundary between no-rub motion and partial rub. Since the synchronous full annular rub bifurcates into the partial rub through Hopf bifurcation, the Hopf bifurcation condition of synchronous full annular rub solution will serve as the boundary between the synchronous full annular rub and the partial rub. As observed in experiment [START_REF] Yu | Rotor/seal experimental and analytical study on full annular rub[END_REF][START_REF] Wu | An experimental study of the influence of disk flexibility and rubbing on rotor dynamics[END_REF] and simulations [START_REF] Jiang | Stability analysis of sliding whirl in a nonlinear Jeffcott rotor with cross-coupling stiffness coefficients[END_REF], the partial rub ceases to exist after the dry friction backward whirl was triggered by the imbalance with the increase rotating speed. In [START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF] the physical reason for the onset of dry friction backward whirl was exploited. It is known that the rotor in resonance at a backward whirl "natural" frequency of the coupled nonlinear rotor/stator system is the reason for the onset of the dry friction backward whirl from the partial rub. From this understanding, a method can be devised to determine the approximate onset condition-another existence boundary of the partial rub.

First, the backward whirl "natural" frequency of the coupled nonlinear rotor/stator system is sought. This frequency corresponds to the mode that the rotor whirls backwards. After setting the right-hand side of (4) to zero and then substituting supposed solution, W = H exp(j ω 0 τ ), into it, we solve the frequency with a negative value and the amplitude (H should be greater than R 0 ) as

ω 0 = - ζ μ + β + ζ 2 μ 2 H = R 0 1 -2 ζ μ ζ μ + β + ζ 2 μ 2 when ζ μ < 1 2 √ 1 + β (15)
Secondly, when the rotor runs at a rotating speed near the magnitude of the backward whirl "natural" frequency, ( 4) is transformed into a form that describes the dynamics of the rotor/stator system near the backward whirl "natural" frequency,

W + 2ζ W + βW + (1 + jμ) 1 - R 0 |H | W = R 0 (1 + jμ) 1 |W | - 1 |H | W + Ω 2 e jΩτ ( 16
)
where the left-hand side of the complex system has the natural frequency equal to the above backward whirl "natural" frequency.

Thirdly, an approximated linear equation to ( 16) is set up after defining

ε = D f = 1 |W | - 1 |H | = O(ε) < 0 ( 17
)
as a small parameter. The reason is that |H | is the amplitude of the rotor in the backward whirl rubbing without imbalance excitation and |W | is the amplitude of the rotor in the backward whirl rubbing with imbalance excitation. Since Ω < 1 in the vicinity of the critical rotating speed of the onset of dry whip, the magnitude of the excitation Ω 2 in (4) is thus a small quantity. So |W | should be only a small modification to |H | and the difference between the corresponding inverse relations should be also small. So ( 16) can be written approximately as:

W + 2ζ W + βW + (1 + jμ) 1 - R 0 |H | W = εR 0 (1 + jμ)W + Ω 2 e jΩτ ( 18
)
The approximated equation can be solved by, e.g., the Multiple Scale Method [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], since W at the rotating speed near the magnitude of the backward whirl "natural" frequency, namely Ω = |ω 0 | + εσ where σ = O(1), is a small detuning.

Assume that the problem depends on many time scales:

T 0 = τ, T 1 = ετ, T 2 = ε 2 τ, . . .
The solution of ( 18) can be represented by an expansion having the form

W (τ ; ε) = W 0 (T 0 , T 1 , . . .) + εW 1 (T 0 , T 1 , . . .) + • • • (19)
The expansion to O(ε 2 ) is carried out in this paper and the corresponding equations of ε 0 and ε 1 are

ε 0 : D 2 0 W 0 + 2ζ D 0 W 0 + κW 0 = Ω 2 exp(j ΩT 0 ) (20) ε 1 : D 2 0 W 1 + 2ζ D 0 W 1 + κW 1 = -2D 0 D 1 W 0 -2ζ D 1 W 0 + R 0 (1 + jμ)W 0 ( 21
)
where

κ = β + (1 + jμ)(1 -R 0 /|H |)
. By this moment, the range of the excitation frequency (the rotating speed) Ω is not specified as done in the analysis procedure of the Multiple Scale Method at the Secondary Resonance [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. The solution of (20) yields

W 0 = a(T 1 ) exp j (ω 0 T 0 ) + θ(T 1 ) + A exp(j ΩT 0 ) ( 22
)
where a(T 1 ) and θ(T 1 ) are the real functions of T 1 and

A = Ω 2 /(κ -Ω 2 + j 2ζ Ω).
Substituting ( 22) into (21) leads to

D 2 0 W 1 + 2ζ D 0 W 1 + κW 1 = -2 jω 0 a -ω 0 aβ -2ζ a + jaβ × exp j (ω 0 T 0 + θ) + R 0 (1 + jμ)a exp j (ω 0 T 0 + θ) + AR 0 (1 + jμ) exp(j ΩT 0 ) (23)
where the prime denotes the derivative with respect to T 1 . Now, substituting the above expression of Ω into (23) and setting the coefficient of terms with exp(j ω 0 T 0 ) on the right-hand side of (23) to zero in order to eliminate the secular term, leads to a complex equation

-2 jω 0 a -ω 0 aθ -2ζ a + jaθ + R 0 (1 + jμ)a × exp(j θ ) + AR 0 (1 + jμ) exp(j σ T 1 ) = 0 (24)
Separating the real and the imaginary parts and then solving a and aθ , gives

a = C 1 (μω 0 + ζ )a + C 1 C 2 sin(σ T 1 -θ + α) aθ = C 1 (μζ -ω 0 )a -C 1 C 2 cos(σ T 1 -θ + α) ( 25
)
where

C 1 = R 0 2(ω 2 0 +ζ 2 ) , C 2 = R 0 (1 + μ 2 )(ω 2 0 + ζ 2 ) × |A|, α = arctan ζ ω 0 + arctan μA r +A i A r -μA i .
A r and A i are the real and imaginary parts of A.

Next, ( 25) is transformed into an autonomous system with T 1 not explicitly appearing by defining

γ = σ T 1 -θ + α (26)
and substituting into (25) to get

a = C 1 (μω 0 + ζ )a + C 1 C 2 sin(γ ) aγ = σ a -C 1 (μζ -ω 0 )a + C 1 C 2 cos(γ ) ( 27
)
The steady-state motions occur when a = γ = 0, which corresponds to the singular points of (27). They correspond to the solution of

C 1 (μω 0 + ζ )a 0 = -C 1 C 2 sin(γ 0 ) σ a 0 -C 1 (μζ -ω 0 )a 0 = -C 1 C 2 cos(γ 0 ) ( 28
)
|W 0 | with a and θ by ( 28) and ( 26) can serve as an approximation of |W | in [START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF]. By using |H | = H in ( 15), D f in ( 17) can be calculated for the given parameters and ε.

Finally, the critical rotating speed for the onset of the dry friction backward whirl or the loss of stability of the partial rub can be solved after the deviation εσ to |ω 0 | can be determined. To do this, we can first set σ a reasonable constant number, i.e., σ = 1 as used in this paper or a value that depends on the system parameters as given in [START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF], then use a trial and error process on ε in order to find the most proper ε that makes D f (as given above) a largest negative value.

The solution in the phase of self-excited backward whirl

Dry friction backward whirl is a motion state of the rotor/stator system, by which the rotor is in continuous contact with the stator, rolling or slipping continuously on the contact surface and whirling backward at a super-synchronous frequency. Dry friction backward whirl is a self-excited vibration of the rotor/stator contact system due to the presence of dry friction effect at the contact surface. Although the mechanism of dry friction backward whirl is well understood, most of the research works used the governing equations, which do not include the dry friction effect to explore this phenomenon. A little progress has been made in analytically determining the characteristics and the parameter region of dry friction backward whirl until now. Below, we will only sketch the solution procedure.

To make the reader to better understand the meaning of ( 5) and the methodology we adopted to approximate the amplitude of dry friction backward whirl, the mechanism for a rotor/stator system that performs dry friction backward whirl is first illustrated.

When rotor-to-stator contact occurs, the friction force at the contact point tends to drive the rotor into backward whirl. In case that the rotor whirls backward, the friction force will apply in the same direction as the rotor whirl and continuously turns the input rotational energy into the energy for the lateral vibration of the rotor system. The vibration amplitude of the rotor is driven higher and higher. The relative velocity at the contact point will change its sign when the rotor deflection is sufficiently large to make the whirl speed exceed the circumferential speed of the rotor at the contact point, i.e., r = -ωr disk /ω w (refer to Fig. 1 [START_REF] Bently | Forced subrotative speed dynamic action of rotating machinery[END_REF]). At this time, the friction force also inverts its direction and is opposite to the rotor whirl to dissipate the energy. The rotor deflection will be driven down. When the rotor deflection decreases to make the relative velocity at the contact point to change sign again, the rotor deflection will be driven up again. In this way, the amplitude of rotor deflection at dry friction backward whirl is kept oscillating around the value with a zero relative velocity, that is, the rotor undertakes a dry friction induced self-excited vibration-dry friction backward whirl.

(b)) or ω b |W | ≥ ΩR disk (see
Now we can use the knowledge of dry friction backward whirl to derive the solution from (3). Since the fluctuation of the response amplitude around the zero relative velocity is much smaller than the total magnitude of the rotor deflection during dry friction backward whirl, the amplitude of dry friction backward whirl can be well approximated by the one at zero relative velocity. From the second equation in [START_REF] Bently | Forced subrotative speed dynamic action of rotating machinery[END_REF], it yields

|W | = -ΩR disk /ω b ( 29
)
By substituting (29) into first equation of ( 5), we obtain the equation governing the motion of dry friction backward whirl in the following form:

W + 2ζ W + βW + (1 + jμ) 1 + ω b R 0 ΩR disk W = Ω 2 e jΩτ (30)
To derive solution for dry friction backward whirl from (30), two important features should be taken into account. (a) Since the contribution of the forced vibration to the total response at dry friction backward whirl is generally too small to be considered, the supposed solution needs only to consider the frequency component with super-synchronous frequency. (b) Most importantly, when the governing equation changes into the form of (30), the solution corresponding to dry friction backward whirl should have a time increasing amplitude, that is, a positive real part in the exponent, because the dry friction effect does not explicitly take effect (its effect is taken into account in (29) to give a constant amplitude). So the solution can be assumed in the form

W 2 = Ce (α+jω b )τ (31)
After substituting (31) into (29), one can get two equations of α and ω b , where ω b indicates the whirl frequency at dry friction backward whirl and α the capability of the system to undertake self-excited vibration. When α is positive, self-excited vibration-dry friction backward whirl-can take place. When α is negative, the responses other than self-excited vibration occur. Therefore, α = 0 is the critical condition for the existence of dry friction backward whirl [START_REF] Jiang | The analytical solution and the existence condition of dry friction backward whirl in rotor-to-stator contact systems[END_REF]. In this way the boundary for dry friction backward whirl in the parameter space can be determined.

Global response characteristics

In this section, two examples are presented to illustrate the overall picture of the global dynamic response characteristics of the rotor/stator contact system with ζ = 0.05, R 0 = 2.00 and R disk = 20R 0 . As can be seen from the above analysis the existence boundaries of no-rub motion, the synchronous full annular rub and the dry friction backward whirl are determined rigorously with very high accuracy. The numerical simulations have proved this point. For the onset condition of dry friction backward whirl from partial rub motion, some approximations are adopted. The validity of the results is also checked by the numerical simulation as given in [START_REF] Jiang | The physical reason and the analytical condition for the onset of dry whip in rotor-to-stator contact systems[END_REF]. So the global dynamical characteristics of the rotor-to-stator contact system will be well demonstrated in the following analysis.

In the following figures, curves Ω l and Ω u indicate the rotating speed whereby the deflection of the linear rotor covers the clearance. Thus, the synchronous no-rub motion is in the rotating speed range less than Ω l and greater than Ω u . Curve HP and line SN are respectively the boundaries of Hopf bifurcation and saddle-node bifurcation of the synchronous full annular rub solution. Curve DW is the boundary by which dry whip is triggered by imbalance and the partial rub ceases to exist. Curve DF gives the existence boundary of the dry friction backward whirl. Also, the parameter regions marked by numbers represent different regions of coexistence, i.e., "0"-no-rub motion + dry friction backward whirl; "1"-synchronous full annular rub + dry friction backward whirl; "2"-partial rub + dry friction backward whirl; "3"-synchronous full annular rub + no-rub motion; "4"-no-rub motion + synchronous full annular rub + dry friction backward whirl.

Figure 2 depicts the global responses characteristics on the parameter planes of Ω-μ with β = 0.04. It is interesting to note that the dry friction backward whirl exists already from a very low rotating speed when the friction coefficient is larger than a certain value which satisfies the condition given in (15) as μ > 0.102. When the friction coefficient is below this value, the only rubbing response is the synchronous full annular rub solution. Due to the coexistence of the no-rub motion and the synchronous full annular rub (see region "3"), the "jump" phenomenon will occur in this region. When the friction coefficient is above this value, one can find the regions that the dry friction backward whirl coexists with other responses. So it is quite possible to induce the more destructive dry friction backward whirl from a relatively mild rubbing response or even from no-rub motion when the rotor is subjected to a sudden outside disturbance. It is even worse that there is a region (enclosed by DW and Ω u ) where dry whip is the sole stable response. Thus, it is important to avoid a large contact friction, especially that which exceeds the one defined by the criterion in [START_REF] Zhang | Dynamic instability of multi-degree-of-freedom flexible rotor systems due to full annular rub[END_REF]. The coexistence of three responses can be also detected in a very small region-"4".

Figure 3 shows the global response characteristics on the parameter planes of Ω-β with μ = 0.15. It is easily seen from Fig. 3 that with the increase of the stiffness ratio (it can be taken as the reduction of the stator stiffness but keeping the rotor stiffness unchanged, namely, using a soft stator), both Ω l and Ω u increase but the latter increases more rapidly. This makes the rotating speed range with the synchronous full annular rub larger and the speed range with the no-rub motion smaller. This seems to be unprofitable to the rotor response. However, it is found from Fig. 3 that when the stiffness ratio is larger than a certain value, i.e., greater than about 0.3 in this case, the Fig. 3 Rotor response characteristics on parameter plane of Ω-β with μ = 0.15, where ζ = 0.05, R 0 = 2.0 and R disk = 20R 0 . The regions of only three of the five typical types of coexistence are found in this case. Ω l and Ω u are the critical rotating speeds of no-rub motion. SN and HP are the saddle-node and Hopf bifurcation boundaries of the synchronous full annular rub. DF and DW are the existence and onset boundaries of dry friction backward whirl synchronous full annular rub will no longer bifurcate into the more severe partial rub. More importantly, the induction of the destructive dry whip through the imbalance will not occur. In fact, when the stiffness ratio is greater than about 2.4, the imbalance can no longer trigger the dry whip although the rotor undertakes the partial rub. From this sense, a larger stiffness ratio benefits the rotor response. Since the criterion given in [START_REF] Zhang | Dynamic instability of multi-degree-of-freedom flexible rotor systems due to full annular rub[END_REF] is met in the range of the stiffness ratio examined in Fig. 3, the dry friction backward whirl exists in the whole region right to line DF. This means that the dry friction backward whirl coexists with all other rotor responses we have examined. From the point of view of the global responses, there are only three types of coexistence regions in this case, except for region "3" and region "4" mentioned above.

Conclusions

In this paper, the global response characteristics of a piecewise smooth rotor/stator contact system are investigated. An effective method to get the global response characteristics of this specific model has demonstrated that the piecewise smooth equations of motion can be analyzed piece by piece so that the existence boundaries of the typical responses in each phases of motion, i.e., in the phase without rub, with rub and of self-excited backward whirl, have been determined. Then, they are incorporated into the same parameter space to get the global response characteristics of the rotor/stator contact system. The results show deep insights into the interactive effect of different parameters on the system responses as well as the coexistence behavior. The main conclusions can be summarized as follows.

From the point view of global responses, five typical types of coexistence are found for the piecewise smooth rotor/stator contact system in this paper. They are the coexistence of no-rub motion and synchronous full annular rub, of no-rub motion and dry friction backward whirl, of synchronous full annular rub and dry friction backward whirl, of partial rub and dry friction backward whirl and of no-rub motion, synchronous full annular rub and dry friction backward whirl. Except for the first type of coexistence, the other types of coexistence are not quite familiar and occur, in fact, when the criterion given in ( 15) is met. From this criterion and also from diagrams of the global response characteristics in the parameter planes, it is known that small friction on the contact surfaces and large stiffness ratio of the rotor-to-stator stiffness will benefit the rotor rubbing behavior by avoiding the occurrence of the severe partial rub or the destructive dry whip.

C

  Real constant in the solution part of backward whirl motion H Complex constant in the supposed solution R 0 Non-dimensional clearance R Non-dimensional deflection of the shaft center R disk Non-dimensional radius of the disk W Complex deflection of the rotor, X + jY X, Y Non-dimension displacements of the shaft center c Damping of the rotor e Rotor mass eccentricity k s , k b Stiffness of the rotor shaft and the stator m Mass of the rotor r 0 Clearance between rotor and stator r Displacement of the shaft geometric center r disk Radius of the disk at the contact point t Time x, y Horizontal and vertical displacements of the shaft center α Real part of exponent of the solution part of backward whirl motion β Stiffness ratio of rotor-to-stator, k s /k b ϕ Phase angle μ Coefficient of friction τ Non-dimensional time ω Rotating speed of the rotor ω 0

Fig. 1 a

 1 Fig. 1 a Schematic diagram of Jeffcott rotor with the stator clearance. b The section of the rotor at the position of disk to show the forces applied on the rotor during the backward whirl of the rotor

Fig. 2

 2 Fig. 2 Rotor response characteristics on parameter plane of Ω-μ with β = 0.04, where ζ = 0.05, R 0 = 2.0 and R disk = 20R 0 . The regions of the five typical types of coexistence are identified. Ω l and Ω u are the critical rotating speeds of no-rub motion. SN and HP are the saddle-node and Hopf bifurcation boundaries of the synchronous full annular rub. DF and DW are the existence and onset boundaries of dry friction backward whirl

Acknowledgements

This work is supported by the National Natural Science Foundation of China (NSFC) under Grant 10472086 and by the Project for Excellent Young Teachers of MOE, China.


Mehran Naghizaderokni

Mohamadreza Hassanlo

Mohamadreza Ramzani

A developed procedure for predicting the risk of liquefaction: A case study of Rasht City

Keywords: Probability of liquefaction, SPT, Safety factor, Empirical equation

Soil liquefaction is one of the most important and complex seismic geotechnical and engineering topics. Most experimental methods in the study of liquefaction have been based on deterministic analysis, and parameters, such as soil resistance and earthquake loads, have been considered without dispersion and fault. Statistical analysis, particularly reliability analysis, is a new and comprehensive approach used to solve and evaluate problems, as well as entering of uncertainties in calculations. In this study, information on 50 boreholes were obtained using a database of Standard Penetration Test (STP), and diverse experimental procedures, such as NCEER2001, Idriss and Boulanger, Highway bridge of Japan, OCDI and reliability approach, were used in the evaluation of soil liquefaction in Rasht and finally, the results were compared with each other. By comparing the different methods, it was observed that these methods do not correspond with clay soils and their results were different from each other. Also, safety factor of greater or lesser than one is the mean of safety and occurrence of liquefaction, therefore, reliability analysis was used to ensure the probability of liquefaction. Finally, two empirical relationships based on the probability of liquefaction (PL) and safety factor (FS) were proposed from which the liquefaction potential can be calculated directly. Soil liquefaction risk can be assessed by the proposed relationship between (PL) and (FS) based on deterministic approaches.

INTRODUCTION

Iran is among the countries located in the earthquake belt, and the nature of the soil in many areas indicates high liquefaction potential under dynamic wind effects of earthquake. Since all the soils do not have the potential to be liquefied, the first step in liquefaction risk assessment is usually to evaluate the liquefaction potential [START_REF] Seed | Simplified procedure for evaluating soil liquefaction potential[END_REF]. Most experimental approaches, such as the Idriss approach, NCEER2001 procedure, Highway bridge of Japan and the Overseas Coastal Area Development Institute of Japan (OCDI), are based on a definitive analysis in evaluating liquefaction, and parameters including soil resistance and earthquake loads are considered precisely, without dispersion and faults [START_REF] Cetin | Standard penetration test-based probabilistic and deterministic assessment of seismic soil liquefaction potential[END_REF]. In addition, there is no exact correlation between the probability of liquefaction and safety factor in the use of these methods. Statistical analysis, particularly reliability analysis is a new and comprehensive approach to studying and solving problems, as well as entering of uncertainties in calculations [START_REF] Juang | Assessing probability-based methods for liquefaction potential evaluation[END_REF]. Different methods of reliability analysis, such as the second moment of the first order (FORM), Estimation of point (PEM), Monte Carlo simulation (MCS) and new method combined FROM-PEM, can be used in the evaluation of soil liquefaction. In this study, information on 50 boreholes were obtained from the database of standard penetration test (SPT) as shown in Figure 1. Different experimental procedures such as NCEER2001, Idriss and Boulanger, Highway bridge of Japan, OCDI method and reliability procedure were used in the study of soil liquefaction in Rasht, and the results were compared with each other [START_REF] Juang | Reliability-based method for assessing liquefaction potential of soils[END_REF]. It is very important to assess the liquefaction potential of Rasht city because of the sandy soils present in some areas and earthquake prone area as well as high levels of ground water that can lead to saturated soil. Thus, one of the main objectives of this study is to identify the liquefaction potential segments and determine the most appropriate method of evaluating liquefaction potential [START_REF] Bahrainy | Urban planning and design in a seismic-prone region (the case of Rasht in Northern Iran)[END_REF].

Fig.1 Location of boreholes in Rasht City

International Journal of GEOMATE, Jan., 2017, Vol. 12, Issue 29, pp. 2878-2884 Geotec., Const. Mat. & Env, ISSN: 2186-2982(Print), 2186-2990(Online), Japan International Journal of GEOMATE, 000., 0000, Vol. 00, Issue 00, pp. 0000-0000

GENERAL GEOLOGY

Gilan plain includes an alluvial area located on the northern slopes of the Alborz Mountains and the southern coast of the Caspian Sea. This plain is divided into three narrow plains: eastern, northwestern and central plain. Rasht city is located in the central part of Gilan plain at a distance of approximately 30 km from the Caspian Sea.

Fig.2 Location of Rasht City in Gilan plain

Rasht city is located at 37° and 1 min to 37° and 27 min of the north latitude from the equator and 48° and 35 min to 49° and 3 min of longitude from the meridian of the Caspian Sea and its area is about 136 km. The existing sediments in Rasht city is related to the fourth age of geology. Rasht city is located on the confluence of 5 types of sediment including coastal sediments, new marine sediments, wetland sediments, the new alluvium delta and marine sediments in terms of geology [START_REF] Dezhkam | Simulating the urban growth dimensions and scenario prediction through sleuth model: A case study of Rasht County, Guilan, Iran[END_REF].

The density of existing clay soils in Rasht city is high (SPT average number equal to 15) and the soils are relatively pre-consolidated (with OCR Average equal 2.75). The natural moisture content of the clay varies due to changes in groundwater level, and its average is 24.4%. The plasticity limit, plasticity index and shrinkage limit of clay in Rasht city are 54, 29 and 13.4%, respectively. Also, the existing grain soil varies in terms of density. The cohesion of these soils is in the range of 0 to 2.0 kg/cm2 (mean 1.0), the angle of internal friction is in the range of 17 to 44° (mean 26) and the standard penetration number is in the range of 8 to 50 (mean 18) [START_REF] Samadi | Shallow subsurface geology and Vs characteristics of sedimentary units throughout Rasht City, Iran[END_REF]. In a research conducted by the Geotechnical Group International Institute of Earthquake Engineering, the liquefaction zoning maps of Iran were plotted, where Rasht city was shown to have high liquefaction capacity.

THE METHODS USED TO ASSESS LIQUEFACTION BASED ON NSPT

The safety of factor was calculated at different depths and the soil of Rasht was analyzed to identify the liquefiable and non-liquefiable areas [START_REF]The Technical Committee for earthquake geotechnical engineering, ISSMGE[END_REF]. All records were then analyzed to identify areas of liquefaction and non-liquefaction by the Overseas Coastal Area Development Institute of Japan (OCDI) procedure. Unlike previous approaches, the main criterion of this method is the possibility of liquefaction and a definitive amount cannot be calculated. Through this method. (N1)60 should be calculated using the proposed equation to predict liquefaction for records in Rasht. σ ) is the effective overburden pressure of the subsoil calculated in the approach at different depths in Rasht. In the following step, the equivalent acceleration was calculated using Eq. ( 4) [START_REF] Iwasaki | Soil liquefaction studies in Japan: state-of-the-art[END_REF]. max ' 0.7

eq v A g τ σ = (2)
This Eq. is used in the OCDI method where ( max τ ) is the maximum shear stress (KN/m 2 ) and g is the gravitational acceleration. Liquefaction can be predicted utilizing the equivalent N-value and equivalent acceleration. The soil layer was categorized according to a diagram that was divided into four segments. A total of 50 diagrams were drawn for each bore log and the liquefaction potential was evaluated based on (Aeq) and (N65) for different records in Rasht [START_REF] Tokimatsu | Evaluation of settlements in sands due to earthquake shaking[END_REF].

Performance function is recognized as one of the major elements in reliability analysis. Some of the uncertainty in determining the (CRR) and (CSR) are also considered in the order provided by the (R) and (S) and can act as random variables which are normally distributed. Performance function can be defined as (Z=R-S) [START_REF] Juang | Firstorder reliability method for probabilistic liquefaction triggering analysis using CPT[END_REF]:

If ( 0 Z < ) the performance functions is designated as 'failed'. If ( 0 Z > ) the performance functions is designated as the safe zone and liquefaction does not occur in that area, and if (Z=0) the performance state is designated as a 'limit state'. After obtaining some of the uncertainties inherent in the estimation of (CSR) and (CRR), the (R) and (S) can be accepted as random variables. The liquefaction probability is defined as the probability that ( 0 Z 〈 ) [15]. However, an exact evaluation of this probability is not an easy task. Indeed, it is so tough to accurately determine the (PDFs) of random variables, such as (R) and (S). A simplified calculation approach, the first order and secondmoment method, has been improved to meet this need. The method utilizes the statistics of the basic independent random variables, such as (R) and (S). If the probability density functions (PDF) and the cumulative probability function (CPF) of (Z) are defined as (fz (Z)) and (Fz (z)), respectively, the liquefaction probability (PL) then equal the probability of ( 0 Z 〈 ). If the mean values and standard deviations of (R) and (S) are ( R µ ) , ( S µ ), ( R σ ), ( S σ ) according to the first order and second moment method, the mean, standard deviation and coefficient of variance of the (Z), by using the following Eq [START_REF] Wotherspoon | Assessment of liquefaction evaluation procedures and severity index frameworks at Christchurch strong motion stations[END_REF]:

2 2 S R S R Z Z Z σ σ µ µ σ µ δ + - = = (3)
The reliability index ( β ) is defined as the inverse of the coefficient of variation ( z δ ), and is used in calculating the reliability of the results of liquefaction assessment. It is assumed that (R) and (S) are independent variables with a normal distribution to display the process of the reliability analysis. According to this assumption, the performance function (Z) can be in a normal distribution of ( ( )

2 , z z Z µ δ 
). By placing the (PDF) of (Z), for (PL):

) ( 1 ) ( β β Φ - = - Φ = L P (4)
In the aforementioned Eq., (PL) represents failure probability, ( z σ ) represents the standard deviation, ( z µ ) represents the mean value, ( β ) represents the reliability index and ( ( ) Φ β ) represents the cumulative probability. A computer program was written in the MATLAB environment to assess the liquefaction potential based on reliability for approximately 500 records in the study area. The seismic information, mean, and coefficient of variation associated with effective parameters were introduced to the program as fixed input parameters to assess liquefaction and the parameters relevant to the genetic algorithm. The geotechnical information required for assessment of liquefaction involves groundwater level, soil bulk density, depth of soil layers, soil type, (SPT) number, the percentage of fines (sieve 200) as variable input parameters.

The application reliability index values were calculated from boreholes to different depths based on the flowchart designed and the probability of liquefaction was then evaluated in these depths.

RESULTS AND DISCUSSION

Due to the variety and heterogeneity of existing deposits in Rasht city on the surface and at different depths, the best way to evaluate the characteristics of liquefaction in the city is by data collection. For this purpose, conducted reports and geotechnical studies in Rasht city on construction projects, such as private buildings, administrative and bridges, were collected and the required data were extracted. The present information resources were developed using a database of the Standard Penetration Test (SPT) related to 500 data with 50 boreholes, which were conducted by the Consultant Engineers Company of Iran soil, Moshaver Asas Gostar Khak and Pay, and Advisor of Gill Design and Control in different areas of Rasht City. In this study, plan acceleration was selected (3.0 g) based on Iran 2800 Regulation (Third Edition), and the Richter magnitude of the earthquake was determined as 7. After analyzing the geotechnical data of Rasht city, 2 boreholes were selected for discussion from the different segments of Rasht city based on engineering judgment. To compare the results of the liquefaction potential assessment with different approaches, the maximum liquefaction depth obtained from the different methods were compared with each other. Based on the above description and considering the safety factor of less than or equal to 1 for NCEER2001, High bridge of Japan, and Idriss and Boulanger approaches as a measure of liquefaction, the results of the field operations, the charts and the comparison of the results of the liquefaction potential assessment in the selected boreholes are shown in Figures 3 to 6 At a depth of 2 m, in the first boreholes having sand with bad granularity, a safety factor of NCEER2001, Idriss and Boulanger, and Highway Bridge of Japan methods were identified as close to 1, and this indicates that liquefaction will occur. However, according to reliability procedure, it has no liquefaction potential but based on the OCDI method, liquefaction is possible.

At 10 m depth of silty sand, the results of analyzing different methods are similar and the soil is considered as non-liquefiable. At depths of 14, 16 and 18-m in which the layers of soil are made of clay and silt, the results were completely different from each other and as such, an accurate and definitive idea could not be determined. Fig. 6 The results of comparison of different ratios of cyclic resistance based on liquefaction potential evaluation with different methods in the fourth boreholes

In boreholes with depths of 2.5 and 5.5 m, having sand with good granularity and silty sand, respectively, the results were similar and the soils were recognized as non-liquefaction and liquefaction, respectively. For boreholes with a depth of 8 m, having sand with bad granularity, the probability of liquefaction indicated that liquefaction incidence and liquefaction nonoccurrence are equally probable, thus NCEER2001, Idriss and Boulanger, and the OCDI methods indicated that liquefaction will occur, whereas a comparison of the calculated highway bridge of Japan procedure suggested that there will be no occurrence of liquefaction. At a depth of 10 m, the results could not be trusted due to the presence of clay soil with high plastic properties. For boreholes with a depth of 14 m having sand with bad granularity, a comparison of the calculated OCDI approach showed that the soil is in the 3 area and the possibility of liquefaction is low. On the other hand, the Idriss and Boulanger method suggested that liquefaction will not occur. However, the results obtained from the reliability approach showed that the possibility of liquefaction and nonliquefaction is equal.

EXTRACTED EQUATION FROM PROPOSED MODEL

The final step was performed with the help of the GP.NET software (Genetic programming). Program Settings were carried out in the setting menu and it was possible to adjust the number of equation constants. The type of relations considered was adjusted in the Function menu. The program became set to begin the analysis by entering the Run menu program. Finally, the output equation became available after completion of the analysis time in the Model menu. The output equation derived from liquefaction assessment based on the result of the Highway bridge of Japan and the GP.NET Software is shown in Relation [START_REF] Boulanger | Probabilistic Standard Penetration Test-Based Liquefaction-Triggering Procedure[END_REF] and Figure [START_REF] Samadi | Shallow subsurface geology and Vs characteristics of sedimentary units throughout Rasht City, Iran[END_REF] In this study, the potential of liquefaction in Rasht city was evaluated using NCEER2001, Idriss and Boulanger, the Highway bridge of Japan, OCDI and reliability procedures.

According to seismic zone, due to the diversity of alluvial sediments, high level of groundwater, and saturated soil, this area has a high potential for liquefaction. Therefore, the liquefaction potential assessment of this area is important to avoid probable dangers and hazards, especially in times of severe earthquake.

CONCLUSION

By collecting and testing 500 reports of geotechnical boreholes in different depths of Rasht city, about the public profile of Rasht city soil and the existing geotechnical properties, the liquefaction potential of different layers can be explained as follows: the values were first obtained from different areas in Rasht with no absolute ideal distribution; and secondly, the values were obtained by relying on the results of other studies.

1. By comparing the liquefaction potential assessments based on the different methods, it was observed that the climatic conditions of Iran have some similarities with the United States. Moreover, the main criteria of Idriss and Boulanger's method is based on (N1)60 and the obtained results have high accuracy. Therefore, it can be concluded that Idriss and Boulanger's method is a better option for liquefaction potential assessment in Iran.

2. In the Highway bridge of Japan approach and OCDI method due to the climatic conditions of Japan and the soil type of this area, the soil liquefaction ability is assumed to be highly conservative. This is not affordable in terms of economic issues. In the OCDI procedure, which is based on the classification of areas to assess liquefaction possibility, the importance of structure can play a significant role in liquefaction diagnosis in areas II and III.

3. Analysis of the reliability approach is recognized as the best means of providing a logical framework for considering the uncertainty in assessing liquefaction potential. The main result of the liquefaction calculation, based on the reliability method, is that the calculation of probability of liquefaction incidence is possible qualitatively using a confidence coefficient.

4. Therefore, based on the results obtained in the present study, an appropriate method should be used based on the geotechnical area and climatic conditions. Thus, the utilization of both probabilistic and deterministic procedures in International Journal of GEOMATE, 000., 0000, Vol. 00, Issue 00, pp. 0000-0000 assessing liquefaction potential seems to be more reasonable and conservative in achieving more definite liquefaction potential.

5. In this research, it was determined that a confidence coefficient greater or lesser than 1 does not mean safety and/ or liquefaction in cadence for liquefaction, hence, in order to ensure liquefaction probability, reliability-based method analysis should be used.

6. Regarding the proposed relationship between (PL) and (FS), the liquefaction probability of soil layers can be obtained using the deterministic methods. This is a big advantage for geotechnical engineers who make use of common methods based on confidence coefficient for liquefaction potential assessment.

  this equation, (N65) is the equivalent N-value, (N) is the N-value of the subsoil, and ( ' v

Fig. 3 Fig. 4

 34 Fig.3 Specifications of second boreholes (Rasht -Golbagh Namaz -next to the Administration Laboratory)

Fig. 5

 5 Fig.5 Specifications of the second boreholes (Rasht -Seyghalan Square).

2 [ 2 (

 22 Fig.7 Graph of safety factor showing changes in the behavior of Highway Bridge of Japan method with probability using GP software

Fig. 8

 8 Fig.8 Graph of safety factor showing changes in the behavior of Idriss and Boulanger's method with probability using the GP software

Table 1

 1 . Summary of data and calculation of liquefaction in the studied area (boreholes 1)

	P	l	β	Aeq	N65	FS 3	FS 2	FS 1	CR R3	CR R2	CR R1	CSR 3	CSR 2	CS R1	FC	N	Soil	De pth	R o w
	28	0.23	437	39.3	0.8 1	0.9 5	0.9 5	0.32	0.2	0.21	0.39	0.38	0.4	3	15	SP	2	1
	12	1.17	431	42.1	0.8 3	1.3 5	0.1 9	0.32	0.26	0.26	0.38	0.34	0.4	3	18	SP	4	2
	45	-0.07	425	31.8	0.7 7	1.1 4	1.1 2	0.28	0.19	0.24	0.37	0.3	0.3 9	9	15	SP	6	3
	3	1.42	401	33.6	1.1	2.4 0	2	0.38	0.3	0.42	0.34	0.22	0.3 7	37	19	SM	10	4
	65	0.11	352	33.3	1.6	1.1 3	0.6 1	0.52	0.29	0.11	0.32	0.16	0.3 2	80	18	CH	14	5
	48	0.14	323	35.4	1.4 9	4.0 3	0.7 9	0.46	0.32	0.13	0.31	0.13	0.2 9	65	22	ML	16	6
	88	-0.56	295	20.2	1.4 6	2.6 9	2.2 4	0.43	0.18	0.34	0.3	0.11	0.2 7	85	16	CH	18	7

Table 2

 2 Summary of information and calculations of liquefaction in the studied area (boreholes 2)

	P	l	β
	58	0.01	435	28.6	0.68	0.73	0.69	0.26	0.15	0.15	0.40	0.37	0.40	3	11	SP	2.5	1
	19	0.84	426	36.6	1.03	1.59	1.85	0.38	0.28	0.41	0.39	0.31	0.39	42	15	SP	5.5	2
	35	0.43	415	36.6	0.89	1.62	1.14	0.32	0.32	0.24	0.38	0.26	0.38	3	19	SP	8	3
	27	0.46	401	44.5	1.49	2.99	1.27	0.52	0.62	0.22	0.37	0.25	0.37	63	25	CH	10	4
	39	0.32	352	26.4	0.82	1.78	0.93	0.26	0.16	0.16	0.32	0.16	0.32	3	18	SP	14	5
	19	0.84	396	21.9	0.90	2.69	1.39	0.27	0.18	0.21	0.27	0.11	0.27	22	14	SP	18	6


Mohamdreza Hossanlourad

Mehran Naghizadeh Rokni

Mohamadreza Hassanlo

Akbar Badrlou

DISPERSIVE CLAY STABILISED BY ALUM AND LIME

Keywords: Dispersive Clay, Alum, Lime, Unconfined Compression Test, Direct Shear Test

Dispersion phenomenon can cause erosion problems during hydraulic construction projects such as earth dams, channels, hydraulic installations and embankments. However, aluminum sulfate can be used to stabilize dispersive clay. Tests were done to consider the effects of alum and lime additions on shear strength parameters of clay as secondary effects. Clay for samples was taken from the Mirzakhanluo Dam, located in Tarom County, northwestern Zanjan Province of Iran. In this research, the clay sample was subjected to direct shear testing; the clay was tested with optimum moisture under varying percentages of alum (0, 0.3, 0.6, 1, 3, and 5 as percentage of dry soil weight) for 28 days. Also, unconfined tests were conducted on samples with equal amounts of moisture, alum and lime under the different treatment durations of 7, 14, and 21 and 28 days, separately. Results for the direct shear test indicated that the maximum cohesion was recorded in the sample with 1% alum. Results also showed that the internal friction angle decreased in samples with alum percentage up to the level of 0.6% and at levels higher than that it showed an increase. Generally, results showed that soil shear strength was higher in soil stabilized by alum than in unstabilized soil. Results of the uniaxial test, determined that alum and lime additions to dispersive soil and its treatment, compression strength was not limited and elastic modulus showed an increase. The highest evaluations were determined for amounts of elastic modulus and uniaxial strength in soil with 3% lime or 5% alum on the 28 th day of treatment.

INTRODUCTION

Type of soil is an important consideration in selection for resources. A lack of understanding of soil type has caused many problems in hydraulic construction projects and in some cases destruction like shrink-swell can occur. Clay soil in many parts of the world can contain plate-like minerals that contribute to shrink-swell characteristics. Dispersion is another common phenomenon in some types of soil [START_REF] Yukselen | Suitability of the methylene blue test for surface area, cation exchange capacity and swell potential determination of clayey soils[END_REF].

Soil in which its mineral content is dispersed when in contact with water is called dispersive clay. This type of clay soil is easily eroded and will flow under conditions of low-level stress and low hydraulic gradient [START_REF] Foster | The statistics of embankment dam failures and accidents[END_REF]. Clay usually contains calcium, potassium and manganese; dispersive clay contains sodium [START_REF] Sherard | Piping in Earth Dams of Dispersive Clay (Paper introduced by Norman L. Ryker)[END_REF].

Dispersion in clay occurs through complicated chemical and physical mechanisms and is related to the absorption quality of its particular mineral structure and ion exchange capacity [START_REF] Sherard | Piping in Earth Dams of Dispersive Clay (Paper introduced by Norman L. Ryker)[END_REF]. Contact with water activates these mechanisms that cause the minerals in dispersive clay to become suspended [START_REF] Ouhadi | Assessment of the stability of a dispersive soil treated by alum[END_REF]. This dispersion phenomenon can cause piping in earth dams, erosion and destruction of roads and water channels and destruction of foundations [START_REF] Sherard | Piping in Earth Dams of Dispersive Clay (Paper introduced by Norman L. Ryker)[END_REF].

Aluminum sulfate, known as alum or white alum, is an important aluminum compound. Aluminum sulfate is used in water purification and as a mordant in dyeing and textile printing. In water purification, it causes impurities to coagulate into larger particles and then to settle at the bottom of a container (or be filtered out) more easily.

According to previous studies, white alum has higher solubility in water than lime and lower cost while being mixed with soil in order to reduce dispersion and is less harmful than lime on contact with human skin. However, using this material has certain disadvantages including high cost, high acidity and CO2 that causes problems associated with plant growth and corrosion [START_REF] Yukselen | Suitability of the methylene blue test for surface area, cation exchange capacity and swell potential determination of clayey soils[END_REF].

Studies have reported use of aluminum sulfate in a dam project for stabilizing part of an embankment in a 152m channel near Sabradinho City in Brazil [START_REF] De Rezende | The use of quarry waste in pavement construction[END_REF]. In Iran, the effect of aluminum sulfate was evaluated in terms of dispersive rate in the Gotvand Dam project [START_REF] Nateghi | Control negative effects of blasting waves on concrete of the structures by analyzing of parameters of ground vibration[END_REF]. According to these above-mentioned investigations, addition of 0.6 % aluminum sulfate changes soil to a non-dispersive state [START_REF] Turkoz | The effect of magnesium chloride solution on the engineering properties of clay soil with expansive and dispersive characteristics[END_REF]. Also, Jafari reports that addition of 3% International Journal of GEOMATE, Jan., 2017, Vol. 12, Issue 29, pp. 156-162 Geotec., Const. Mat. & Env., ISSN: 2186-2990, Japan; DOI: http://dx.doi.org/10.21660/2017.29.93287 alum to dispersive soil, which was tested in this research, made the soil completely non-dispersive [START_REF] Asgari | Swelling and dispersion potential of soils from geotechnical view[END_REF]. Mohammadi et al. reports that by adding about 1% alum to dispersive clay, dispersion is decreased and consolidation and swelling is improved [START_REF] Asgari | Swelling and dispersion potential of soils from geotechnical view[END_REF].

In this research, white natural alum (AL2 (SO4)3.18H2O) and lime (or calcium hydroxide (Ca (OH) 2)) were mixed with dispersive clay. Effectiveness of the mentioned materials was determined as follows: firstly, by replacing aluminum trivalent cation (AL3+) and calcium bivalent cation (Ca2+) with sodium monovalent cation (Na+), the soil structure changed from disperse to clod, the ejection force between the granules decreased and double layer thickness decreased; finally, with increased cohesive force, there was a decrease in dispersive potential in samples, the granules got closer together and the swelling showed a remarkable decrease [START_REF] Locat | Laboratory investigations on the lime stabilization of sensitive clays: shear strength development[END_REF].

MATERIALS USED AND THE TESTING PROGRAM

The dispersive clay used in this study was from clay deposits of the 2 nd Mirzakhanlou Dam, located 16 km distance from Darram area in Zanjan. According to test results and different standards of evaluating dispersive soil, this soil was ranked as highly dispersive and erosive (Tables 1 and2) [START_REF] Asgari | Swelling and dispersion potential of soils from geotechnical view[END_REF]. The results of hydrometric and aggregation tests (single and double) were obtained. In Figure 1 circles with voids represent a double hydrometric and filled circles represent a single hydrometric.

The double hydrometer test is a good example of misinterpretation due to ambiguity. Despite the test being an ASTM standard, many laboratories simply duplicate the standard hydrometer analysis procedure (ASTM), which invariably produces incorrect results. Variability of the results obtained from the double hydrometer test appears to be the cause of much ambiguity and discrepancy in the classification systems used during this research. Problems relating to the double hydrometer tests have the potential to mislead results since the test is associated with a number of different parameters in the rating systems. Inaccurate results from the double hydrometer test can significantly affect the correlation of the final rating, particularly when this test method is used as the reference method for preliminary classification of soil dispersion. In order to resolve these problems, a mixture of single and double hydrometric results were used. Table 3 shows percentages of granule passing and ASTM standard soil ranking. Figure 2 indicates that clay in the unified ranking system was the plastic type (CH) [START_REF] Jha | Volume change behavior of lime treated gypseous soil-influence of mineralogy and microstructure[END_REF]. The step to produce and treat samples was to apply the density test through application of the proctor compaction test on natural soil (Table 4). Tests were done of percentages of alum mixed with soil and duration of the treatment period (28 days). Dispersive tests showed that after 21 days, the dispersive soil mixed with aluminum sulfate was completely stabilized and that an increased treatment duration had no further effect on stabilization.

RESULTS OF DIRECT SHEAR TEST

Direct shear testing was done according to ASTM D3080 standard with 1mm/min velocity by the CU method (Consolidated Undrained) [START_REF] Locat | Laboratory investigations on the lime stabilization of sensitive clays: shear strength development[END_REF]. Direct shear testing was first performed on dispersive clay without additives and on samples with different percentages of aluminum sulfate (0.3, 0.6, 1, 3, and 5). All samples had maximum density and optimum moisture and were treated for 28 days. After treatment, direct shear tests were done with different loads (0.5, 1, 1.5 kg/cm 2 ).

Results shown in Figure 3 indicate that adding aluminum sulfate increased shear strength of the soil. Variations of apparent cohesion and internal friction angle with alum content are presented in Figures 4 and5, respectively. Figure 4 shows that cohesion increased from 0.3 to 0.5 kg/cm 2 ; i.e. about 67% by adding 1% alum. Further addition of alum to the dispersive clay led to a reduction in apparent cohesion. However, apparent cohesion at alum content above 1% still showed improvement compared to results for unstabilized dispersive clay. Figure 5 demonstrates that addition of alum reduced the internal friction angle from 32.5 to 27.5 degrees; however, at alum content of 5% the internal friction angle increased to 37.5 degrees. In other words, addition of alum beyond 6% increased the internal friction in clay. The uniaxial compression test is a special kind of triaxial shear test. This test can be described as an unconsolidated and undrained triaxial test without comprehensive stress it is used to define Cu cohesion in clay according to Eq. ( 1) below [START_REF] Bro | Estimating undrained strength of clays from direct shear testing at fast displacement rates[END_REF][START_REF] Goodarzi | Stabilization treatment of a dispersive clayey soil using granulated blast furnace slag and basic oxygen furnace slag[END_REF]:

2 u u q C = (1)
Where 𝑞𝑞 𝑢𝑢 the uniaxial compression strength and Cu is is the undrained cohesion of soil.

The results and procedure of the uniaxial test were determined according to ASTM 2166.

The first step was to test dispersive clay without additives, then further tests were applied with additions of, 0.3, 0.6, 1, 3 and 5% of alum and lime was added to the soil to the soil dry weight, separately. Then, samples were produced in optimized moisture and maximum density (Table 4) and loaded on days 7, 14, 21 and 28.

Results of Unconfined Compression Tests

According to Figure 7, on the 7 th day of treatment, adding alum increased the uniaxial strength, which was at its highest level at 5% alum. In soil samples with lime, the uniaxial strength increased until 1% and, after that showed a decrease until 5%. Figure 8 shows that soil with lime had 183% increase in strength and soil with alum had an increase of 87%. On the 14 th day of treatment, the increasing rate of uniaxial strength still showed growth and increasing alum; 5% caused an 80% increase in strength. In soil samples with lime, the uniaxial strength increased until 1% and then it started to decrease until 5%. The maximum strength reached 190%. According to Figure 9, on the 21 st day of treatment, the uniaxial strength still increased and addition of alum resulted in an increase; at 5%, the strength increase was about 103%. In soil with lime there was also an increase in strength that reached 231% [START_REF] Anggraini | Reinforcement Benefits of Nanomodified Coir Fiber in Lime-Treated Marine Clay[END_REF].

According to Figure 10, the increasing amount of uniaxial strength after 28 days of treatment was twice that of the 5% alum sample. In soil with lime, the increasing rate continued until 3% and reached about 5 times as much and after that it showed a decrease until 5%.

In Figures 11 and12, variations of the uniaxial strength versus treatment time are shown with different alum and lime percentages; results indicate increased uniaxial strength in soil stabilized with alum and lime. Aluminum has good ability for transmitting ions on the first day, so the strength of growth of alum-stabilized soil was faster at the beginning of the treatment and the growth rate decreased after the 28 th day.

Comparison of the Elastic Modulus among Soil Samples

Figure 13 indicates that adding alum increased evaluations of elastic modulus, showing a %50𝑞𝑞 𝑢𝑢 increase. In all diagrams, an increasing rate of elastic modulus was higher in the lime-stabilized soil than in the alum-stabilized soil [START_REF] Sivapullaiah | Gypsum induced strength behaviour of fly ash-lime stabilized expansive soil[END_REF][START_REF] Kavak | Long-term behavior of lime-stabilized kaolinite clay[END_REF]. The highest value of elastic modulus for the alumstablized specimens was 5% alum content with a 37% increase over the untreated clay, while the highest elastic modulus for the lime-stablized samples 409% increase over evaluations for untreated clay.

CONCLUSION

According to studies on dispersive soil stabilized by alum and lime and using uniaxial and direct shear tests, the following conclusions were made:

1-Adding alum at different percentages both controlled soil dispersion and increased its shear strength. Adding alum caused increased shear strength of up to 50%.

2-Increasing alum content to 1% led to an increase of up to 61% in apparent cohesion of the dispersive clay. Any further increase in alum content led to a reduction in apparent cohesion.

3-Increasing alum content caused a reduction of the internal friction angle of the dispersive clay. The minimum friction angle of the dispersive clay was reduced. The minimum friction angle was observed at 0.6% alum content, after which any further increase of alum content resulted in an increased internal friction angle.

4-Under unconfined compression tests, both strength and stiffness evaluations of the dispersive clay were improved through addition of lime or alum.

5-Aluminum has high ion transmission ability on the first day so the increasing rate of the axial strength was faster in the lime modified soil and this growth rate decreased after the 28 th day. For soil modified by lime, the highest amount of uniaxial compression strength was with 5% lime, which was about twice that recorded for unmodified soil.

6-The largest increase in elastic modules (of up to 5 times that of the untreated clay) was observed for the specimen with 3% lime content.

Fig. 1 Fig. 2

 12 Fig. 1 Particle size distribution by simple and double hydrometric methods [6, 11] Table 3The results of used soils granular [7, 12] AA SH TO Clas s Gs Pas 0.005 doubl e hydro metri c Pas 0.005 simpl e hydr o metri c D60 mm D3 0 m m D1 0 m m

Fig. 3 Fig. 5 Fig. 6

 356 Fig. 3 Mohr -Coulomb shear failure envelope for different amounts of alum

Fig. 7 Fig. 8 Fig. 9

 789 Fig. 7 Unconfined compression strength versus amount of additives (7 days' curing time)

Fig. 13

 13 Fig. 13 Elastic modulus versus amount of additive (28 days' curing time).

Table 1

 1 Properties of the used dispersive clay[START_REF] Turkoz | The effect of magnesium chloride solution on the engineering properties of clay soil with expansive and dispersive characteristics[END_REF]

	Classification	Classification	Percentage of
	based on	based on	dispersiveniss
	pinhole test	double	by double
		hydrometric	hydrometric
		test	test
	D1: Totally	Totally	59.20
	dispersive	dispersive

Table 2

 2 Chemical specifications of the used dispersive clay[START_REF] Turkoz | The effect of magnesium chloride solution on the engineering properties of clay soil with expansive and dispersive characteristics[END_REF]

	P H	Na	+	K	+	Ca	2	+	Mg	2	+ E C	T D	C L	S O
		meq		meq	meq			meq			S	+	4
		lit		lit		lit			lit			%		2 -
	7.	174		9 34.			13.9 9	0.	10.	38
	2	.3				8					6	6	22	.9

Table 4

 4 Variations of the optimum moisture content and maximum dry density[START_REF] Turkoz | The effect of magnesium chloride solution on the engineering properties of clay soil with expansive and dispersive characteristics[END_REF][START_REF] Anggraini | Reinforcement Benefits of Nanomodified Coir Fiber in Lime-Treated Marine Clay[END_REF]

	Percentage
	of added
	alum as percentage	0	0.3	0.6	1	3	5
	of soil dry
	weight
	Percentage
	of the	30	28.	27.	27.	28.	27.
	Optimum	.4	1	7	3	4	9
	Moisture
	Maximum
	Dry
	Density
	Of the
	Soil
	Gr/cm 3

ACKNOWLEDGEMENTS

I am using this opportunity to express my gratitude to everyone who supported me to complete this paper. I am thankful for their aspiring guidance, invaluably constructive criticism and friendly advice during the project work.