M Nicolas

A comprehensive study on the behavior of a rigid block on an oscillating ground with friction, elastic and viscous forces

Keywords: friction, oscillator, seismic response

This paper investigates the behavior of a non-linear mechanical model where a block is driven by an oscillating ground through Coulomb friction, a linear viscous damper and a linear spring. The governing equation is solved analytically for different partial configurations: friction only, friction with viscous damping, friction with a linear restoring force, and for the complete model. Using dimensionless groups, the analysis of the block motion provides a comprehensive set of information on the motion regime (stick, stick-slip or permanent sliding), on the dominant energies or forces, on the resonance and on the amplification of the ground oscillation by the system. The limit between the stick-slip regime and the permanent slipping regime is found either analytically or numerically. It is also shown that there exists a set of parameters for which the friction force, the viscous dissipative force and the elastic restoring force are equal.

Introduction

While the mass-spring-dashpot system is extensively presented in physics textbooks, the addition of a strongly non linear force such as the Coulomb friction force leads to a more complex problem. Indeed, when an external energy is brought to the system not only by the spring but also through the dry friction between two solid surfaces, the friction force is both a damping force and a driving force.

The coupling of friction and an external vibration is interesting in many physical and mechanical systems. A vertical pile of rigid blocks is a paradigm for many old built structures such as greek columns [START_REF] Papaloizou | Komodromos Planar investigation of the seismic response of ancient columns and colonnades with epistyles using a custom-made software[END_REF], medieval churches [START_REF] Lancioni | Dynamics and failure mechanisms of ancient masonry churches subjected to seismic actions by using the NSCD method: The case of the medieval church of S. Maria in Portuno[END_REF] or five-stories wooden pagodas in Japan [START_REF] Nakahara | Earthquake response of ancient five-story pagoda structure of Horyu-Ji temple in Japan[END_REF]. When the ground is shaken because of seismic activity, the friction force between the blocks may not be sufficient to prevent the collapse of the structure. Much more recently in the building history, the development of anti seismic buildings [START_REF] Coladant | Seismic isolation of nuclear power plants -EDF's philosophy[END_REF][START_REF] Syed | The challenging requirements of the ITER anti seismic bearings[END_REF] is reached by the use of soft and dissipative links between beams and columns, as sketched in Fig. 1(a): a reinforced elastomer supplies elasticity while a dashpot supplies viscous dissipation between the pillar and the beam. Another recent system for seismic-response reduction is the friction pendulum system (FPS) where a contact element coated with teflon slides on a spherical concave surface (as sketched on Fig . 1b). A simple non-linear mechanical model of the horizontal restoring force implies a Coulomb friction force as well as an elastic restoring force [START_REF] Landi | Comparison of different models for friction pendulum isolators in structures subjected to horizontal and vertical ground motions[END_REF]. Despite these solutions are already implemented on actual buildings, the understanding of the relative importance of the different mechanical elements is still to establish. This paper aims to bring a fundamental insight on the energy transfer through a friction-elastic-viscous link between an oscillating ground and a rigid block. We focus here on the role of the non-linear friction force on the resonance frequency and on the relative motion between the block and the ground. The transition between the stick regime, the stick-slip regime and the permanent slip regime is investigated both theoretically and numerically with the use of dimensionless groups.

The governing equation and analysis method

The model

The mechanical model of the pillar-beam link is presented on Fig. 1(c).

We consider a single block of mass M linked to a horizontally shaking ground through a friction-viscous-elastic set of forces. We note X the position of the block relatively to the ground, Ẋ its velocity, K the spring stiffness (in A 0 , a frequency Ω/2π, and an initial phase ϕ. This system is slightly different than the model studied by Hundal [START_REF] Hundal | Response of a base excited system with Coulomb and viscous friction[END_REF], López [START_REF] López | Energy dissipation of a friction damper[END_REF] or Chatterjee [START_REF] Chatterjee | Resonant locking in viscous and dry friction damper kinematically driving mechanical oscillators[END_REF]: in this present model, the spring and dashpot is a link between the mass and the oscillating ground and not between the mass and a non-moving wall. Taking the oscillating ground as a reference frame, the relative motion equation of the block is

M Ẍ = -µ( Ẋ)M g Ẋ | Ẋ| -KX -γ Ẋ -M A 0 Ω 2 cos(ΩT + ϕ). (1)
In this study, two Coulomb friction models µ( Ẋ) are considered. The first and simplest model is the Coulomb friction model with a unique friction coefficient µ s . The second model introduces a kinematic friction coefficient µ k when the block slides on the ground, usually lower than the static friction coefficient µ s . The table 1 summarizes the values for models 1 and 2. Many other friction models exist (for example see the review [START_REF] Hinrichs | On the modelling of friction oscillators[END_REF]) but since they need more than two parameters, they we prefer simple Coulomb models to highlight the role of the friction force.

Introducing the following dimensionless variables and groups

x = X A , ẋ = Ẋ √ A 0 g , k = K M A 0 g , η = γ M A 0 g , α = A 0 Ω 2 µ s g , ω = √ αµ s the motion equation (1) becomes ẍ + µ( ẋ) ẋ | ẋ| + kx + η ẋ + ω 2 cos (ωt + ϕ) = 0. (2)
In this non-dimensional model, the governing parameter is the reduced acceleration α, and the non-dimensional frequency is ω = √ αµ s , which depends on the static friction coefficient. In addition, we assume the following initial conditions : x(0) = 0, ẋ(0) = 0 and we set ϕ = π/2 to ensure an initial zero acceleration of the ground.

Analysis method and motion characteristics

Despite the non-linear behavior of the friction force, the motion equation ( 2) is piecewise linear, and the block either sticks or slides on the ground.

With the given initial conditions (see above), the block starts to slide at a time t 0 > 0 when the inertia force overcomes the static friction force.

An analytical solution of Eq. ( 2) may be found, and since two integration constants arise, these constants are solved using the initial conditions x(t 0 ) and ẋ(t 0 ) = 0. The block then stops at a time t 1 > t 0 when its velocity is zero again ẋ(t 1 ) = 0. This time t 1 is numerically solved using a Newton method.

At t = t 1 the block either sticks to the ground if the sum of the inertia and elastic forces is lower than the static friction force, or slides again with a velocity sign change if the sum of the inertia and elastic forces is greater than the static friction force. The first case leads to a stick-slip motion and the block sticks until time t 2 > t 1 where the start condition is reached again i.e. when the sum of the inertia force and the elastic force is greater than the static friction force. The second case leads to a permanent slipping motion and the next start time is simply t 2 = t 1 , with t 2 -t 1 = π/ω, the half-period of oscillation of the ground. The sequel of the motion is calculated with the same method as previously with a new start time t 2 instead of t 0 for the computation of the integration constants, and subsequently the analytical solution of Eq. ( 2) is given through a list of {t i } times, where even i are start indexes, and odd i are stop indexes. An illustration of the t i times is given on Fig. 2(a) for a stick-slip motion. We checked our analytical method with a numerical solving of Eq. (2) using an implicite numerical scheme with a10 -3 time-step. Both methods were implemented using Python packages and codes.

The block behavior is characterized with different measurable quantities.

First, we measure the the block amplitude A which is half of the peak-topeak block displacement x. We also introduce the time-ratio ζ which is the time when the block slips on the ground compared to the total time:

ζ = ∆t slip ∆t slip + ∆t stick = (t i+3 -t i+2 ) + (t i+1 -t i ) t i+4 -t i ( 3
)
where i is the index of a start time. A permanent stick regime is thus observed for ζ = 0 (ST), a stick and slip motion for 0 < ζ < 1 (SS) and a permanent slip regime is observed for ζ = 1 (SL). A last relevant quantity is the root mean squared velocity defined by

ẋrms = 1 τ t+τ t ẋ2 (t )dt . (4)

Forces and energies

The relative importance of the three forces in the system can be assessed through non dimensional groups. Assuming that a steady state is established with a bounded oscillation amplitude, a generic solution is x(t) = A(µ s , µ k , k, η)F (ωt) where F is a periodic bounded function with |max(F )| = 1. We therefore introduce

ρ ev = k 2 η 2 αµ s , ρ vf = η 2 A 2 αµ s µ 2 k , ρ ef = kA µ k (5)
which are the ratio of the elastic force to the viscous force, the ratio of the viscous force to the friction force, and the ratio of the elastic force to the friction force respectively.

By multiplying Eq. ( 2) with the relative velocity ẋ and integrating over time, the energy conservation equation is

t 0 ẍ ẋdt + t 0 ẋxdt +µ k t 0 | ẋ|dt +η t 0 ẋ2 dt +αµ s t 0 ẋ cos (ωt + ϕ) dt = 0 (6)
which may be written as

e k + e p + e f + e v = -w i (7)
where e c is the kinetic energy, e p is the potential energy, e f is the friction dissipated energy, e v is the viscous dissipation energy, and w i is the work of the inertia force.

If the time integration is taken over one cycle (between t and t + τ , τ = 2π/ω), there is no variation of the kinetic and potential energies because of the periodicity of the x(t) and ẋ(t) functions, and the energy conservation equation is reduced to

∆e f + ∆e v + ∆w i = 0 (8)
The first term ∆e f is the energy dissipated by Coulomb friction during sliding. Over one cycle of oscillation, it is equivalent to

∆e f = µ k t+τ t | ẋ(t )|dt = µ k t+τ t |dx| = 4µ k A (9)
configuration friction viscous damping elasticity damping

configuration 1 • configuration 2 • • configuration 3 • • configuration 4a • • • weak configuration 4b • • • critical configuration 4c • • • strong Table 2: Configurations studied in this paper.
where A is the amplitude of the block motion. The second term is the energy dissipated by viscous dissipation over one cycle:

∆e v = η t+τ t ẋ2 (t )dt = ητ ( ẋrms ) 2 . (10)
Finally, the the third term of Eq. ( 8) is the work of the inertia force

∆w i = ω 2 t+τ t ẋ(t ) cos (ωt + ϕ) dt (11)
which is the energy input per cycle of ground oscillation.

In this paper, we investigate separately the role of the friction force, of the elastic restoring force and of the viscous dissipation force, and we seek to highlight the effects of these forces on the amplitude A, and on the time-ratio ζ depicting the sliding regimes. The different configurations are summarized in table 2. All the results presented in this paper correspond to a permanent periodic motion, i.e. when the transient motion has vanished.

Results for configurations 1-3

In this section we study three partial models where only the friction force applies (configuration 1), the friction force and the viscous dissipative force are present (configuration 2) and where the friction force is coupled with the elastic restoring force (configuration 3).

Configuration 1: friction only

When the friction is the only applied force (k = 0, η = 0), it plays the ambiguous role of the driving force and the damping force. The Eq. ( 2) then

reduces to ẍ = -µ( ẋ) ẋ | ẋ| -ω 2 cos (ωt + ϕ) . (12)
Such a mechanical model is known to present two acceleration thresholds.

The non-dimensional acceleration threshold between the sticking regime and the stick-slip regime is simply α SS = 1. The criterion α α SL for the permanent slipping regime may be calculated as follows. Naming

t 0 = ω -1 [arccos(-µ s /ω 2 ) -ϕ] (13)
the time when the block starts to slide and t f the time when the block stops sliding, the sliding duration t f -t 0 is governed by

µ k (t f -t 0 ) + ω [sin(ωt f + ϕ) -sin(ωt 0 + ϕ)] = 0. ( 14
)
Since the block is sliding at any time when t f -t 0 = π/ω (half of a period), this gives a criteria on the vibration parameter

α SL = 1 + π 2 4 µ k µ s 2 . ( 15
)
With a unique friction coefficient (µ s = µ k ), this model has been proposed as early as the 1930's by Den Haartog [START_REF] Den | Forced vibrations with combined Coulomb and viscous friction[END_REF], then by Alspaugh [START_REF] Alspaugh | Analysis od Coulomb friction vibration dampers[END_REF], and extended by Westermo & Udwadia [START_REF] Westermo | Periodic response of a sliding oscillator system to harmonic excitation[END_REF] with α SL = 1 + π 2 /4 ≈ 1.86 (see López [START_REF] López | Energy dissipation of a friction damper[END_REF] or Benedetti et al. [START_REF] Benedetti | An angular particle sliding down a transversally vibrated smooth plane[END_REF] for the analytical calculation of α SL ).

During the sliding, the solution of Eq. ( 12) is

x(t) = c 1 + c 2 t - 1 2 µ k t 2 + cos(ωt + ϕ) = 0 ( 16
)
where the constants c 1 and c 2 are the solutions of respectively. With this friction-only model, the block amplitude A increases with the vibration parameter α from the unity threshold and tends towards unity for a large acceleration parameter α (Fig. 2f), whatever the friction model. There is no obvious discontinuity on the amplitude curves when the motion regime changes from stick-slip (SS) to permanent slip (SL). For a larger acceleration α > α SL , the block slides continuously (Fig. 2c,d), the direction change is indicated by a slope discontinuity on the velocity curve.

  1 t 0 0 1     c 1 c 2   =   x 0 + 1 2 µ k t 2 0 -cos(ωt 0 + ϕ) ẋ0 + µ k t 0 + ω sin(ωt 0 + ϕ)   (
0 1 2 3 t/τ 1 0 x (a) α =1.2 ζ =0.57 2A t 0 t 1 t 2 t 3 t 4 0 1 2 3 t/τ 0.5 0.0 0.5 ẋ (b) 0 1 2 3 t/τ 3 0 x (c) α =2.0 ζ =1 0 1 2 3 t/τ

Configuration 2: friction and viscous dissipation

When the viscous dissipation term is non-zero, without the elastic restoring force (k = 0), the governing equation is

ẍ + µ( ẋ) ẋ | ẋ| + η ẋ + ω 2 cos (ωt + ϕ) = 0. ( 18
)
with an explicit solution for t t 0 (when sliding)

x(t) = c 1 - c 2 η e -ηt - µ k ηω (ωt + ϕ) + 1 η 2 + ω 2 ω 2 cos(ωt + ϕ) -ηω sin(ωt + ϕ) (19)
where the constants are the solutions of

  1 -1 η e -ηt 0 0 e -ηt 0     c 1 c 2   =   x 0 + µ k ηω (ωt 0 + ϕ) -1 η 2 +ω 2 [ω 2 cos(ωt 0 + ϕ) -ηω sin(ωt 0 + ϕ)] ẋ0 + µ k η + ω 2 η 2 +ω 2 [ω sin(ωt 0 + ϕ) + η cos(ωt 0 + ϕ)]   (20)
We show on Figs. 3(a,b) the α SL (η) curves for both friction models. For a weak dissipation parameter (η 2), these curves may be well approximated by

α SL (η) = 1 + (-0.48 µ s µ k + 1.18)η 1 + π 2 4 µ k µ s 2 (21)
with a linear increase of the threshold with the viscous parameter η. This expression is plotted as blue dotted lines in Figs. 3(a,b).

Compared to the configuration 1, adding a viscous dissipative force helps to reduce the block amplitude, as shown on Fig. 3(c), computed for α = 2.5. This threshold value depends on α and the limit curve separating these two dissipation regimes is plotted as a dotted line on Figs 3(a,b). The viscous parameter η has a strong effect on the characteristic velocity ẋrms (Fig. 3f), very similar to the amplitude.

Configuration 3: friction + spring

When an elastic restoring force is coupled with the friction force, the motion equation is

ẍ + µ( ẋ) ẋ | ẋ| + kx + ω 2 cos (ωt + ϕ) = 0. ( 22
)
and its solution is

x(t) = c 1 cos( √ kt) + c 2 sin( √ kt) - µ k k - ω 2 k -ω 2 cos(ωt + ϕ) (23)
where the constants c 1 and c 2 are the solutions of

  cos( √ kt 0 ) sin( √ kt 0 ) - √ k sin( √ kt 0 ) √ k cos( √ kt 0 )     c 1 c 2   =   x 0 + µ k k + ω 2 k 2 -ω 2 cos(ωt 0 + ϕ) ẋ0 -ω 3 k 2 -ω 2 sin(ωt 0 + ϕ)   (24)
In this configuration, the friction force is the only dissipation force when sliding, while the static friction force and the elastic force are driving the

in the SL regime, where the oscillation amplitude increases indefinitely (Fig. 4e).

Despite an analytical solution for this configuration, it would be cumbersome to derive an analytical expression for the permanent slip vibration parameter α SL as a function of k. A numerical solution is thus proposed here.

The figure 5(a,b) summarize the motion regimes, showing the α SL (k) curves for the two friction models in the (k, α) space. The curve separating the SS and the SL regimes has a non trivial shape. For 4.1 < k < 5.3 with friction model 1 or 3.7 < k < 5.4 with friction model 2, an increase of α from 1 causes the system to have a SS/SL/SS/SL sequence. Moreover, for a low value of k, with the friction model 1 (Fig. 5a), the permanent slip parameter α SL decreases to a minimum then increases with a nearly constant slope. With the friction model 2, the α SL (k) curves also decreases but intercepts the α = 1 permanent stick boundary. This shows that for 0.18 < k/µ s < 0.41 the block can not show a stick-slip motion. The limit ρ ef = 1 is drawn as a magenta curve (Fig. 5a) and shows that the elastic force is dominant except for a very low value of the elastic parameter (k < 0.2) or for a low acceleration. For the friction model 2, the same trend is observed for the elastic parameter (the friction force is dominant for k < 0.15), whatever the acceleration parameter.

Without a strong dissipation force, the amplitude diverges at the reso- with β as a growth coefficient depending mainly on k and weakly on the friction model (Fig. 5f).

This configuration is relevant for the friction pendulum system, where the stiffness is K F P S = M g/R (M is here the mass of the pillar and R is the curvature radius of the FPS, see Fig. 1b). In our model, the dimensionless stiffness is then k = A 0 /R, the ratio of the horizontal ground displacement on the curvature radius. In an actual FPS, the curvature is of the order of one meter, and a horizontal ground displacement is typically 1 to 10 cm during an earthquake. This leads to a 10 -2 to 10 -1 value for k.

Results for the complete model

We now consider the solutions of Eq. ( 2) without neglecting any force.

For a simple forced and damped oscillator, it is known that three motion behaviors exist, depending on the respective values of the damping parameter η and the spring stiffness k. With our symbols, the critical damping is expressed by

η 2 -4k = 0 (27)
and we provide in the appendix the explicit solutions of Eq. ( 2) for the underdamped behavior (η 2 -4k < 0), the critically damped oscillator (η 2 -4k = 0) and the overdamped behavior (η 2 -4k > 0).

A first result is a mapping of the block behavior in the (k, α) plane for η = 0.5 (Fig. 6). For α 1, the relative motion is non zero and the boundary between the stick-slip and the permanent slip regime is drawn as a red curve. We also present as a light gray area the space where the block amplitude is larger than unity, i.e. when the block motion is amplified compared to the ground motion. This area is roughly centered around the undamped resonance line (dashed black line). The limit between the underdamped and the over-damped behavior is drawn as a vertical blue dashed line. By comparing the dimensionless forces (see Eqs. 5), we plot ρ ef = 1 (magenta), ρ vf = 1 (orange) and ρ ev = 1 (green). Thus, there exists an area where the friction force is dominant (light blue area), an area where the viscous dissipation force is dominant (light red area) and finally an area where the elastic restoring force is dominant (dashed area). It is interesting to remark that these three areas connect as a triple point (black circle) where the three forces are of equal magnitude. This occurs for an acceleration α and a stiffness k = η √ αµ s .

As usual, the resonance is reached when then block amplitude reaches a maximum at a resonance acceleration α r . The blue curve represents the resonance line (k r , α r ) where the block oscillation is maximum. For a large k value, this line is asymptotic to α = k/µ s , drawn as a dashed black line, the resonance condition without viscous dissipation.

Since this mapping depends on η, we present on Fig. 7 resonance line

α =k/µ s ρ vf =1 ρ ef =1 ρ ev =1 η 2 -4k =0
triple point frame.

This mapping helps to characterize the block motion for any set of mechanical parameters {α, k, η, µ s , µ k }. In particular, a good estimate of these parameters for an actual system may help to predict the occurence of a stickslip behavior or of a permanent slip behavior. Indeed, the stick-slip motion presents periodic singularities for the acceleration and thus produces strong force fluctuations. And this model also helps to predict the range of parameters for which the amplitude of the block is A > 1, thus when the system amplifies the ground motion.

While the results presented in Figures 6 and7 were computed for the friction model 1, the block amplitude (Fig. 8a

Conclusions

We proposed an analytical and numerical study of the motion of a block linked by a spring, a dashpot and a friction contact to a rigid oscillating ground. For each configuration, we provided the solution of the motion equation, then analyzed these solutions by computing the slipping time ratio, resonance line

α =k/µ s ρ vf =1 ρ ef =1 ρ ev =1 η 2 -4k =0
triple point and comparing the energies and forces for a continuous range of parameters.

From this work, the main conclusions are:

• the slipping limit has a non-trivial expression as soon as k is non zero The solution x(t) for Eq. ( 2) can be analytically derived during the slipping motion for η 2 -4k > 0:

•
x(t) = c 1 e -s 1 2 t + c 2 e -s 2 2 t - µ k k - ω 2 D (k -

Figure 1 :

 1 Figure 1: (a) Sketch of a damping link between a column and a beam in a anti seismic building. (b) Sketch of a simple friction pendulum system. (c) Corresponding mechanical model studied in this paper.

  17)An example of the block trajectory and its corresponding velocity is shown on Fig.2a, showing a typical stick and slip motion for α = 1.2. It must be noted that even with a single friction coefficient (red curves), the stick-slip regime exists. This model is indeed different from the traditional masson-belt model where the stick-slip phenomenon occurs only for a kinetic friction coefficient different (and lower) than µ s . With the friction model 2 (black curves), the block shows a larger displacement amplitude and a larger peak velocity (Fig.2b) due to the lower friction coefficient during the sliding motion. We show the influence of α on ζ and A on Fig.2(e) and (f)

Figure 2 :

 2 Figure 2: Configuration 1. (a) position and (b) velocity for α = 1.2, (c) position and (d) velocity for α = 2.0. The red curves are computed for the friction model 1, and the black curves are computed for the friction model 2. Full symbols indicates a stop of the motion whereas empty symbols indicates a start of the motion. For this example the static friction coefficient is µ s = 0.25. (e) solid line: sliding time-ratio ζ increasing with the acceleration parameter α, for friction model 1 (red curve) and model 2 (black curve). Dashed line: characteristic velocity ẋrms . (f) the block amplitude A varying with the α vibration parameter for the friction model 1 (red curve) and the friction model 2 (black curve). The stars indicates the threshold between the stick-slip regime and the permanent slip regime.

Figure 3 : 5 .

 35 Figure 3: Configuration 2. (a) the permanent slip threshold α SL (bold red line) and the ρ vf = 1 line (dashed red line) as a function of the viscous parameter η for friction model 1. (b) the permanent slip threshold α SL (bold black line) and the ρ vf = 1 line (dashed black line) as a function of the viscous parameter η for friction model 2. The blue dotted lines in (a,b) represent the expression (21) for both friction models. (c) Block amplitude vs η for friction model 1 (red) and friction model 2 (black). (d) The variation of the energy dissipated by the dashpot during one vibration cycle. (e) The ratio ∆e v /∆e f as a function of η. The stars indicate ∆e v = ∆e f . (f) The mean rms velocity as a function of η. Plots (c)-(f) are computed with α = 2.5.

  block and thus are injecting energy into the system. Since the friction dissipation force is finite, the block may experience a long transient before a steady oscillating state is reached. Some typical behaviors are shown on figure4for two different accelerations, α = 1.2 (a,b,c) and α = 3 (d,e,f), and for the friction model 1. A first obvious effect of the elastic restoring force is to put the block back to a zero mean position, even if the block is displaced at the start of the motion. While in the SS regime the transient is short and a steady-state is reached after a few cycles, the transient depends on the stiffness k value in the permanent slip regime. For k < k r , the transient is short (Fig.4d), while for k > k r the transient may last more than 10 oscillation cycles (Fig.4f). Here the threshold k r value is a resonant stiffness simply expressed by k r, SL = αµ s ,

Figure 4 :Figure 5 : 1 .

 451 Figure 4: Configuration 3. examples of the motion of the block x(t) in the stick-slip regime (a, b, c: α = 1.2) and in the permanent slip regime (d,e,f: α = 3). Curves are computed for the friction model 1. The full symbols indicates a stop of the motion whereas empty symbols indicates a start of the motion.

  the behavior of the system for friction model 1 for different values of η. Above a weak value of η, the α SL (k) curve becomes single-valued (see Fig.7bfor η = 0.05). With an increase of η, the α SL threshold increases, following the same trend as for the configuration 2.An increase of the viscous parameter implies obviously an extension of the viscous-dominated area (light red area), as shown on Fig.7(e) for η = 1. In this configuration, the triple point lies in the stick-slip domain and the A > 1 area as well as the resonance line does not lie in the (0 < k < 6, 0 < α < 6)

Figure 6 :

 6 Figure 6: Configuration 4. Map of the block behavior in the (k, α) plane for η = 0.5. Friction model 1. The light gray area sketches the condition where A > 1. The friction force is dominant in the light blue area, the viscous dissipation force is dominant in the light red area, and the the hashed area represents the dominance of the elastic force. Also plotted is the triple point where the three forces are of equal magnitude.

  , b) is weakly affected by the choice of the friction model, and for a weak dissipation (η = 0.01), the A(α) curve shows the same fluctuations as in configuration 3, above and below the resonant acceleration α r . The magnitude of these sub-harmonic and superharmonic resonances increases with the stiffness value (Fig. 8a for k = 1 and Fig. 8b for k = 2). The influence of the viscous dissipation is clearly shown on Fig. 8c (k = 1) or Fig. 8d (k = 2), where the η parameter is varied from 0.01 to 0.75. When the dissipation increases, the A(α) fluctuations are damped, and for η > 0.1 they even are undistinguishable.

Figure 7 :

 7 Figure 7: Configuration 4. Effect of η on the motion regime in the (k, α) space for different values of η. For all plots, the friction model 1 is chosen. The light gray area sketches the condition where A > 1. The friction force is dominant in the light blue area, the viscous dissipation force is dominant in the light red area, and the the hashed area represents the dominance of the elastic force. Also plotted is the triple point where the three forces are of equal magnitude.

Figure 8 :

 8 Figure 8: Configuration 4. The block amplitude vs α for the friction model 1 (red curves) and the friction model 2 (black curves) and for = 1.0, η = 0.01 (a) and k = 2.0, η = 0.01 (b). Effect of η for k = 1.0 (c) and k = 2.0 (d) for friction model 1. All curves correspond to an underdamped behavior.

4 6. 1 .k k - ω 2 D 2 √µ k k + ω 2 (k+ω 2 ) 2 2 √ 2 √ 6 . 3 .

 41222222263 the friction model µ( ẋ) has a weak effect on the behavior and on the block amplitude• sub-harmonic and super-harmonic resonances due to the coupling of the friction and the stiffness • for a given set of parameters {α, k, η, µ s , µ k }, this model and its solving method provide a fast answer on the block behavior regarding the stickslip or permanent slipping regime, its amplitude and the dominant force in play.With this comprehensive knowledge of the dominant energies and forces, and of the motion regimes (stick-slip or slipping regimes), this simple model of a single block may be extended to a vertical pile of identical or different blocks.6. Appendix: analytical solutions for configuration Configuration 4a: underdampingWhen η 2 -4k < 0, we introduceω k = 1 2 4k -η 2 , D = (k -ω 2 ) 2 + η 2 ω 2(28)so that the solution of the motion equation (2) writesx(t) = e -η 2 t [c 1 cos(ω k t) + c 2 sin(ω k t)] -µ (k -ω 2 ) cos(ωt + ϕ) + ηω sin(ωt + ϕ) (29)where c 1 and c 2 are the solutions ofe -µ k k + ω 2 D [(k -ω 2 ) cos(ωt 0 + ϕ) + ηω sin(ωt 0 + ϕ)] ẋ0 -ω 3 D [(k -ω 2 ) sin(ωt 0 + ϕ) -ηω cos(ωt 0 + ϕ)]  (30)witha 11 = ω k sin(ω k t 0 ) + η 2 cos(ω k t 0 )(31)a 12 = ω k cos(ω k t 0 ) -η 2 sin(ω k t 0 )(32)a 21 = cos(ω k t 0 )(33)a 22 = sin(ω k t 0 ) (34)6.2. Configuration 4b: critical dampingWhen η 2 -4k = 0, the solution of Eq. (2) isx(t) = (c 1 + c 2 t)e - ω 2 ) 2 (k -ω 2 ) cos(ωt + ϕ) + kω sin(ωt + ϕ) (35)when the block slips. The constants c 1 and c 2 are to be determined with the (k -ω 2 ) cos(ωt 0 + ϕ) + kω sin(ωt 0 + ϕ) ẋ0 + ω 3 (k+ω 2 ) 2 -(k -ω 2 ) sin(ωt 0 + ϕ) + kω cos(ωt 0 + ϕ) Configuration 4c: overdamping

ω 2 ) 1 =+ µ k k + ω 2 D

 212 cos(ωt + ϕ) + ηω sin(ωt + ϕ) η + η 2 -4k, s 2 = η -η 2 -4k. (38)The constants c 1 and c 2 are the solutions of [(k -ω 2 ) cos(ωt 0 + ϕ) + ηω sin(ωt0 + ϕ)] ẋ0 -ω 3 D [(k -ω 2 ) sin(ωt 0 + ϕ) -ηω cos(ωt 0 + ϕ)]  (39)

Table 1 :

 1 Description of the two friction models used in this study.

	N•m -1 ), γ the viscous dissipation constant (in kg•s -1 ), and µ( Ẋ) the friction
	law. The ground has a monochromatic horizontal motion with an amplitude


Braden Snook

Elisabeth Guazzelli

Jason E Butler

Vorticity alignment of rigid fibers in an oscillatory shear flow: Role of confinement

Rigid fibers suspended in a viscous, Newtonian fluid at high concentrations can be aligned in the direction perpendicular to the flow-gradient plane (vorticity direction) by applying an oscillatory shear flow. A simple model, which considers only excluded volume and self-mobilities, can accurately predict the orientation distributions measured in experiments by Franceschini et al. ["Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter," Phys. Rev. Lett. 107, 250603 (2011)]. Furthermore, simulations reveal that the alignment of the fibers in the vorticity direction depends strongly on the presence of the bounding walls. C 2012 American Institute of Physics. [http://dx.doi.org/10.1063/1.4770141] The rheological and mechanical properties of suspensions composed of rod-like particles are sensitive to the spatial and orientational distribution of the particles within the suspension. Understanding and controlling the flow-induced alignment of these elongated particles is of fundamental importance and has implications for numerous material processing applications, such as the manufacture of paper 1 and fiber composites. [START_REF] Papathanasiou | Flow-Induced Alignment in Composite Materials[END_REF] Recent measurements [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] indicated that the orientation distribution in semi-concentrated suspensions of non-colloidal fibers can be controlled by altering the strain amplitude within an oscillatory shear flow. The orientation distribution changes little from the initial state at sufficiently small amplitudes of strain. At large strain amplitudes, the fibers align in the direction of the flow, similarly to the steady shearing flow of a suspension of rigid fibers. [START_REF] Stover | Observations of fiber orientation in simple shear flow of semi-dilute suspensions[END_REF][START_REF] Sundararajakumar | Structure and properties of sheared fiber suspensions with mechanical contacts[END_REF] While these results might be expected, the experiments also measured a strongly preferred alignment of the fibers in the vorticity (transverse to both the flow and gradient planes) direction for intermediate strain amplitudes. This surprising result remains to be understood. Alignment with the vorticity direction is observed when oscillating fibers suspended in a weakly elastic fluid. [START_REF] Petrich | Oscillatory shear alignment of a non-Brownian fiber in a weakly elastic fluid[END_REF] However, the experiments 3 of interest here were performed with a Newtonian fluid. For an elastic fluid, the orientation of even a single fiber spins toward the vorticity direction, [START_REF] Leal | The slow motion of slender rod-like particles in second-order fluid[END_REF][START_REF] Iso | Orientation in simple shear flow of semi-dilute fiber suspensions. 1. Weakly elastic fluids[END_REF] whereas the orientation of a single fiber oscillated in a Newtonian fluid traces a reversible path defined by the Jeffery orbit corresponding to the initial orientation. [START_REF] Okagawa | Memory impairment in flowing suspensions. I. Some theoretical considerations[END_REF] The irreversible dynamics of suspensions of non-colloidal particles suspended and sheared in Newtonian fluids has been studied extensively; the vast majority of work has focused on spherical particles, but early work used fibers. [START_REF] Okagawa | Memory impairment in flowing suspensions. I. Some theoretical considerations[END_REF][START_REF] Okagawa | Suspensions: Fluids with fading memories[END_REF][START_REF] Ennis | Memory impairment in flowing suspensions. II. Experimental results[END_REF] The reversibility of the Stokes equations implies that the particles should retrace their paths upon reversal of the flow. However, any small source of irreversibility in a physical system (e.g., Brownian motion, particle surface roughness, particle deformability, or inertia) causes a loss of memory. Oscillating the suspensions provides a particularly convenient format for studying irreversibility in suspension systems. [START_REF] Breedveld | Measuring shear-induced diffusion and rheology of noncolloidal suspensions: Time scales and particle displacements[END_REF][START_REF] Pine | Chaos and threshold for irreversibility in sheared suspensions[END_REF] For sufficiently large strains and concentrations, particles do not return to their starting configurations after one or more cycles and their displacements exhibit a chaotic, random motion. The need for relatively large concentrations of particles implies that short-range interactions are the primary origin of the chaotic behavior.

This work addresses the origin of the unexpected and irreversible alignment of the fibers in the direction of vorticity, which cannot be attributed to the elasticity of the suspending fluid. Phys. Fluids 24, 121702 (2012) FIG. 1. (a) Each fiber of aspect ratio A = L/d is described by its center of mass x i and orientation p i . A repulsive force acts between two fibers i and j at the points of closest approach separated by the distance h ij . (b) The x-direction is that of the flow; the flow varies linearly in the y-direction, or gradient direction, and the vorticity axis is in the z-direction. The fiber's projection on the flow-vorticity plane and the flow direction forms the angle α:

cos 2 (α) = p 2 x /( p 2 x + p 2 z ).
Rather the vorticity alignment, and other behaviors of the orientation distribution, was attributed by Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] to short-range interactions between fibers in a process that can drive the system into a reversible, or "absorbing," state. A simple algorithm demonstrates that short-range interactions between particles can account for alignment in the vorticity direction, depending upon the concentration and strain amplitude. The model reveals that the extent of alignment depends strongly on whether or not the suspension of fibers is confined in the gradient direction by closely-spaced bounding walls. Only those simulations with a wall-spacing matching that used in the experiments accurately predict the mean orientation reported by Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] The model is similar to that used by Sundararajakumar and Koch 5 to simulate the steady shear of semi-concentrated suspensions of rigid fibers. The motion of each fiber is determined by a balance between the hydrodynamic drag forces and short-range forces used to maintain the excluded volume. Inertia is ignored in these calculations, since both the maximum Reynolds number and frequency of oscillation were small in the experiments. Long range fluid disturbances, or hydrodynamic interactions, are neglected within the model, as they have little impact on dynamics in the semiconcentrated regime; [START_REF] Sundararajakumar | Structure and properties of sheared fiber suspensions with mechanical contacts[END_REF] comparing results of steady shear from this model with simulations with hydrodynamic interactions 14 confirms the assertion.

The motion of each rigid fiber, i, is fully described by the center-of-mass velocity, ẋi , and the rotational velocity, ṗi , where x i and p i are the center-of-mass position and unit vector parallel to the fiber's major axis as shown in Fig. 1. The center-of-mass motion is

ẋi = u (x i ) + ξ -1 I + p i p i • F i (1)
and rotational motion is

ṗi = • p i + B I -p i p i • E • p i + 12ξ -1 L 2 T i × p i . ( 2
)
The mobility of the center-of-mass is ξ -1 I + p i p i , where ξ -1 = ln (2A)/4πμL, μ is the fluid viscosity, L is the fiber length, A = L/d is the aspect ratio, and d is the fiber diameter. The mobility is derived from slender-body theory. [START_REF] Batchelor | Slender-body theory for particles of arbitrary cross-section in Stokes flow[END_REF][START_REF] Cox | The motion of long slender bodies in a viscous fluid. Part 1. General theory[END_REF] A simple oscillatory shear, u (x, t) = γ (t) ye x , is imposed, with flow in the x-direction and the gradient in the y-direction as shown in Fig. 1. The shear rate, γ (t), is a square wave of magnitude one and the period is determined by the strain amplitude, γ .

In the absence of collisions, each fiber translates with the velocity evaluated at the center-ofmass, u (x i ), and rotates in a Jeffery orbit: [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] the orientation rotates fully with the rate of rotation of the flow, = 1 2 (∇u) -(∇u) T , and only a fraction B = A 2 e -1 / A 2 e + 1 of the rate of extension of the flow, E = 1 2 (∇u) + (∇u) T . The aspect ratio A has been replaced by an effective one, A e , in the parameter B. The effective ratio, as developed through experiment and theory, [START_REF] Bretherton | The motion of rigid particles in a shear flow at low Reynolds number[END_REF][START_REF] Trevelyan | Particle motions in sheared suspensions. I. Rotations[END_REF][START_REF] Mason | Particle motions in sheared suspensions: Orientations and interactions of rigid rods[END_REF] accounts for the rotation of a spherocylinder as opposed to an ellipsoidal particle, for which Jeffery [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] originally derived the equation. Phys. Fluids 24, 121702 (2012) Collisions between the fibers alter the motion of both the center-of-mass and orientation. To prevent the overlap of a fiber i with a fiber j, a short range repulsive force is applied,

f i j = 0 if h i j > , f 0 n i j if h i j ≤ , (3)
where h ij is the minimum separation distance between two fibers and n i j is their common normal, ± p i × p j / p i × p j . The sign of the normal is chosen so that the force repels the fibers and the force on fiber j is equal and opposite to that on i, hence there is no net force. The simulation results are insensitive to the exact form of the repulsive force used, so long as the range is small compared to the fiber diameter and the displacements caused by the force are resolved at a fine level by holding the time step small. Here, the range is set to 0.1d and the force f o = | γ Lξ/2|. On rare occasions when fibers overlap, the force f o is increased relative to the overlap distance, ensuring that the particles separate. All interactions on fiber i sum to give the total, non-hydrodynamic force F i ,

F i = N f j=1 f i j , (4)
and torque T i ,

T i = N f j=1 s i j p i × f i j , (5)
where s ij is the point along fiber i at which the collision occurs with fiber j and the sums exclude the occurrences of i = j. The total number of fibers simulated, N f , is set by the concentration and size of the simulation box. The concentration of fibers and strain amplitude are varied within the simulations. The aspect ratio is set to A = 11 to correspond with the experiments and A e = 8.8. For fully periodic calculations that utilize Lee-Edward's boundary conditions, the simulation box is cubic with a length of 3.25L. For simulations with boundaries, fibers are confined in the gradient direction by applying the same force between a fiber and the wall as used between contacting fibers. The wall separation of 1.45L used in the simulations corresponds to that used in the experiments of Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] Periodic boundaries in the flow and vorticity directions are separated by 5.25L. The number of particles simulated ranges from 265 to 1232, corresponding to volume fractions φ of 0.05 to 0.20. Equations ( 1) and (2) are integrated using an Euler method with an adaptive time step selected to maintain a maximum displacement less than 0.01d. Twelve initial configurations were used to calculate each mean value and reported errors are those between the mean values of each configuration.

Primarily, two measures are calculated from the data generated by the simulations and presented here: one for the average orientations and another that tracks the irreversible activity of the suspension of fibers. The order parameter S α = 1 -2 cos 2 (α) is calculated from the simulation results to facilitate comparison with the experimental measurements 3 of alignment that were reported solely in terms of this parameter. As shown in Fig. 1, the angle α is that between the flow direction (x) and the projection of a fiber's orientation onto the flow-vorticity (x-z) plane; the brackets, • , indicate an average over all runs and fibers. Note that S α is 0 for a suspension having a random orientation distribution and is either -1 or 1 for a suspension of fibers aligned perfectly with the flow or vorticity direction, respectively. The number of contacts between fibers are counted over each oscillation and reported as N c . In this model, irreversible motions of the particles are generated only by contact forces, thus, N c is a measure of irreversible activity within the suspension.

Figure 2(a) compares experimentally measured and calculated values of S α as a function of strain amplitude. Simulations with boundaries separated by a distance equivalent to the experiments captures the trend of the alignment behavior, including the surprising alignment of the fibers with the vorticity direction in the vicinity of γ ≈ 3. Using fully periodic boundaries does not reproduce the experimental results. the initial condition at N = 0 and progressively aligns in the vorticity direction through N = 400, corresponding to the point at which the comparison is made with experiments in Fig. 2(a). The simulation predicts an increasing alignment of the fibers with the vorticity direction until reaching a steady state at N ≈ 4000. Figure 3(a) provides more details about the evolution of S α with N and shows that the value of S α for γ = 3 transitions from 0.60 at N = 400 to 0.70 at N = 4000; the significant difference in the orientation distribution represented by this difference of 0.1 in the value of S α is demonstrated in Fig. 2(b).

Figure 3(a) shows that significant changes in orientation occur when oscillating with γ > 1. The curves exhibit an initial decrease in S α , indicating a movement toward flow alignment for initial oscillations, followed by an increase that is more pronounced for simulations with boundaries. In this latter case, S α attains positive values that reflect a vorticity alignment of the fibers. Note that the initial distribution, both spatial and orientational, is a random one in the simulations, whereas steady shear was applied in the experiments before beginning the oscillations. This may explain the lack of an initial decrease in S α as reported in the third figure of Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] Figure 3(b) shows that an order of magnitude reduction in the number of collisions between fibers occurs for γ = 3 and 5 for simulations with bounding walls. A similar trend is observed for simulations without confinement for the same strain amplitudes, though N c is an order of magnitude larger. In these cases, N c reaches steady state by the same value of N as S α . The number of collisions also decreases with N for γ = 1, but a steady state is not obtained within the same number of cycles as S α . Though S α does not change significantly, individual fibers continue to change orientation due to contacts.

Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] interpreted the lack of a change in the order parameter S α after the first few oscillations for γ < 2.2 to the suspension attaining an absorbing state in which individual particles return to their starting positions and orientations after each cycle. They also associated the vorticity alignment with the attainment of an absorbing state, whereas Fig. 3 indicates activity across all strain amplitudes for the duration of the simulations. The discrepancy between the numerical predictions and experimental measurements of activity could be explained by resolution limits in the measurements. The simulations predict orientation changes and displacements of only a few microns for individual particles at γ = 1 and N > 10; such small changes would be difficult to detect even with a microscope. Also, lubrication interactions that exist between fibers in the viscous fluid could prevent, or at least reduce, the occurrence of contacts. The question of the existence of contacts in the presence of lubricating fluids is a question of active investigation for suspensions of spheres, [START_REF] Arp | The kinetics of flowing dispersions. IX. Doublets of rigid spheres (experimental)[END_REF][START_REF] Metzger | Clouds of particles in a periodic shear flow[END_REF] Phys. Fluids 24, 121702 (2012) though arguments [START_REF] Petrich | Interactions between contacting fibers[END_REF] that contacts are more likely between fibers than spheres have been made. The current simulations do not include lubrication with the goal of generating a minimal model that can explain the unexpected alignment phenomena. All data presented up to this point have been for a volume fraction of 0.20. Figure 4 shows the steady values of S α at φ = 0.20 and two additional volume fractions. At the lower volume fraction of 0.05, the orientation of the fibers tend to align more with the flow, rather than the vorticity, direction for both the bounded and unbounded cases.

Confining the fibers has the most impact on the orientation distribution at φ = 0.20 and γ = 3, where the results of simulations with full periodicity predict a mean alignment in the direction of flow. The probability distribution of the angle α given in Fig. 5(a) shows the small preference for flow alignment when using periodic boundaries and also quantifies, with more detail, the strong alignment in the vorticity direction when using bounding surfaces.

The maximum extent of alignment is predicted at φ = 0.15 and γ = 3, where S α = 0.88 ± 0.02 for the simulations of confined fibers. Using fully periodic boundaries for the same concentration and strain amplitude also generates a preferred alignment in the vorticity direction, where S α = 0.60 ± 0.05. Comparing the details of the distributions, as done in Fig. 5(b), demonstrates the significant orientation difference between the two types of boundary conditions. This latter result shows that bounding walls are not required to achieve some alignment of the fibers in the vorticity direction, but confining the fibers clearly enhances the extent of organization and alignment in the vorticity direction. The alignment phenomena, whether in the flow or vorticity direction, are a complex function of the interactions between the imposed straining flows and their amplitude, the concentration of the fibers, and the geometry of the shear cell.

Confining elongated particles between bounding walls generally promotes organization of the orientations, as for liquid crystals, which transition from an isotropic to nematic phase when squeezed between planar walls. [START_REF] De Gennes | The Physics of Liquid Crystals[END_REF] Likewise, the boundaries impose some order on the particles in the immediate vicinity of the wall (see Fig. 2(b)), even for the initial condition; fluctuations in orientation generated by the oscillatory shear, rather than thermal fluctuations as for liquid crystals, provide the mechanism for the particles to search for a favorable arrangement. Perhaps, the more interesting point here is the preference for vorticity, versus flow, alignment.

The lack of a boundary that organizes the orientations results in a significantly lower extent of, or even no, alignment in the direction of vorticity as compared to the bounded simulations. A simple test was performed to verify that the solid boundaries themselves, and neither the initial orientation distribution created by the walls, nor different geometry of the simulation box, are responsible for the vorticity alignment predicted by the simulations. The initial distributions were generated within the bounded cell, but then the simulation was performed with periodic boundaries. These simulations returned only a slightly lower tendency for the fibers to align in the flow direction. For φ = 0.20 and γ = 3, the simulation predicts S α = -0.08 ± 0.03 at N = 400, which is close to the value of S α = -0.17 ± 0.03 for periodic boundaries when starting from a random particle configuration. The sensitivity of the fiber alignment on the initial distribution was explored further by aligning the suspension in steady shear prior to oscillating, as was done in the experiments. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] For φ = 20 with γ = 3 and bounding walls, the simulations with the pre-shear predict that any differences in the mean alignment disappear by N = 400, where the pre-sheared cases give S α = 0.63 ± 0.05 and the initially random, but confined, distributions give S α = 0.6 ± 0.1.

In conclusion, the minimal model employed here, which considers only self-mobilities and excluded volume of the fibers, reproduces the experimental observation 3 of vorticity alignment of the orientation distribution for oscillatory shear when boundaries are included. Franceschini et al. [START_REF] Franceschini | Transverse alignment of fibers in a periodically sheared suspension: An absorbing phase transition with a slowly varying control parameter[END_REF] attributed the vorticity alignment to particle-particle collisions solely, omitting the wall contributions. Though some alignment can occur without confining the suspension, the confinement clearly plays a strong role in the experimental observations, which utilized a wall spacing of 1.45 fiber lengths. Accordingly, a simple and complete principle of organization is still lacking as the alignment depends subtly upon three factors: the amplitude of the straining flow, the concentration of the fibers, and the boundary conditions. For the latter factor, only the two cases of fully periodic boundaries and a wall separation identical to that in the experiments were considered. The wall separation at which the orientational dynamics of the suspension transitions from the highly confined to the unbounded behavior remains to be explored.

FIG. 2 .

 2 FIG. 2. (a) Comparison of simulation results for S α with the experimental values for φ = 0.20 and at the the last oscillation reported by Franceschini et al. 3 The inset shows the oscillation number (N R ) at which the data were reported. (b) Images of the microstructure from a simulation with boundaries, φ = 0.20, and γ = 3 prior to oscillating (N = 0), at the last point recorded in the experiments (N = 400), and at steady state (N = 4000) (enhanced online) [URL: http://dx.doi.org/10.1063/1.4770141.1] [URL: http://dx.doi.org/10.1063/1.4770141.2] [URL: http://dx.doi.org/10.1063/1.4770141.3].

FIG. 3 .

 3 FIG. 3. Results for a volume fraction of φ = 0.20. (a) The evolution of S α with the number of oscillations for both the periodic and bounded simulations. (b) The total number of collisions, N c , experienced between the fibers during each oscillation.

FIG. 5 .

 5 FIG.5. The probability distribution of the angle α, where 0 o is the flow direction and 90 o is the vorticity direction (see Fig.1), at steady state for bounded and fully periodic conditions. Results for a volume fraction of (a) φ = 0.20 and (b) φ = 0.15. Note that the plots utilize different scales.

  Steady values of the order parameter S α as a function of strain amplitude for simulations performed with fully periodic boundaries (Periodic) and with bounding walls (Bounded).

							Phys. Fluids 24, 121702 (2012)
		1
			Bounded φ=0.05
		0.8	Bounded φ=0.15
			Bounded φ=0.20
		0.6	Periodic φ=0.05
			Periodic φ=0.15
		0.4	Periodic φ=0.20
	α	0.2
	S
		0
		-0.2
		-0.4
		-0.6
		0	1	2	3	4	5
			γ (Strain Amplitude)
	FIG. 4.

The assistance of Mr. Rusty Stiles with running the simulation code and organizing the results is greatly appreciated. We also thank Dr. Franceschini and Professor Pine for their comments and discussion. This work was supported by the Partner University Fund for particulate flows and the National Science Foundation (Grant No. 0968313).