Fluid dynamics of dilatant fluid
aa r X i v : . [ c ond - m a t . s o f t ] D ec Fluid dynamics of dilatant fluid
Hiizu Nakanishi
Department of Physics, Kyushu University 33, Fukuoka 812-8581, Japan
Shin-ichiro Nagahiro
Department of Mechanical Engineering, Sendai National College of Technology, Natori, Miyagi 981-1239, Japan
Namiko Mitarai
Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, DK-2100, Copenhagen Ø, Denmark (Dated: September 11, 2018)Dense mixture of granules and liquid often shows a sever shear thickening and is called a dilatantfluid. We construct a fluid dynamics model for the dilatant fluid by introducing a phenomenologicalstate variable for a local state of dispersed particles. With simple assumptions for an equation of thestate variable, we demonstrate that the model can describe basic features of the dilatant fluid suchas the stress-shear rate curve that represents discontinuous severe shear thickening, hysteresis uponchanging shear rate, instantaneous hardening upon external impact. Analysis of the model revealsthat the shear thickening fluid shows an instability in a shear flow for some regime and exhibits theshear thickening oscillation , i.e. the oscillatory shear flow alternating between the thickened and therelaxed states. Results of numerical simulations are presented for one and two-dimensional systems.
PACS numbers: 83.80.Hj,83.60.Rs, 83.10.Ff,83.60.Wc
I. INTRODUCTION
One of the most common materials of the dilatant fluidis a dense mixture of cornstarch and water, and it canbe used to demonstrate a number of counter-intuitivebehaviors that the shear thickening medium shows: sud-den solidification upon externally applied stress, quickre-fluidization after removal of the stress, formation ofholes and protrusions under strong vibration[1, 2], etc.These behaviors come from severe shear thickening andhysteresis, that dense colloid or dense mixture of granulesand liquid often show. The shear viscosity increases al-most discontinuously by orders of magnitude at a certaincritical shear rate[3], which makes the fluid almost rigidagainst the sudden application of stress. It is called a “di-latant fluid” by analogy with the behavior of a granularmedium[4]; when a granular medium is densely packedin a bag that is flexible but non-stretchable, it cannotbe deformed because the volume is constant. The gran-ular medium must dilate upon deformation due to the principle of dilatancy by Reynolds[5].There are several peculiar features in the shear thick-ening of the dilatant fluid: (i) the thickening is so severeand instantaneous that it might be used even to make abody armor to stop a bullet[6], (ii) the relaxation after re-moval of the external stress occurs within a few seconds,that is quick but not as instantaneous as in the thickeningprocess, (iii) the medium in the thickened state behaveslike a rigid material allowing little elastic deformation aslong as it is under stress, (iv) the viscosity shows hystere-sis upon changing the shear rate[7], (v) noisy fluctuationshave been observed in the response to an external shearstress in the thickening regime[7, 8].Despite of the apparent analogy between the behav-iors shown by these media, it is not clear if the shear thickening of the dilatant fluid has something to do withthe property of dilatancy of granular media. Origi-nally, the shear thickening in colloid systems were re-garded as a result of the disorder transition of the layerand/or string structure developed in the low shear rateregime[9–12]. The dispersed particles align due theshear flow to give shear thinning in low shear regime,but the turbulent motion in high shear regime destroysthis structure to give shear thickening. Such layerand/or string structures have been observed in numer-ical simulations[13] and experiments[14], and in somecases the shear thickening occurs when the structure isbroken[9]. However, there are some other cases where nosignificant structure change are observed upon discon-tinuous shear thickening[15–18]. Hydrocluster formationhas been proposed as an alternative origin of the shearthickening[16, 19, 20]. Due to hydrodynamic interactionamong particles in the fluid, there is a certain condi-tion that clusters of particles grow and they can givelarge viscosity. Such a cluster structure of particles hasbeen first identified in numerical simulations[19], thensuggested by SANS experiments[16]. More direct obser-vation has been made using fast confocal microscopy[21].Jamming is another possibility under active debate inrecent years in connection with the glass transition. Indense granular system, the jamming can cause the diver-gence of viscosity[3, 8, 22–26]. In connection with thedilatancy, Brown and Jaeger studied the discontinuousshear thickening and obtained somewhat empirical con-stitutive relations[27].There are no microscopic theories for the dilatant fluidyet in the sense that the shear thickening is derived fromthe elementary interactions among constituents of themedium, i.e. granules and fluid, but there are a cou-ple of semi-empirical theories: the soft-glassy rheology(SGR) model[28] and the schematic mode coupling the-ory (MCT)[29]. The SGR model is the model based onthe stochastic dynamics with the activation energy thatdepends on the stress. This model is extended to de-scribe the shear thickening by introducing the stress de-pendent effective temperature. MCT, which gives reason-able description of the glass transition, has been extendedschematically by introducing shear rate dependent in-tegral kernel. Both of the theories are semi-empiricaland have been demonstrated to show the discontinuousshear thickening, but they have not been incorporated inthe fluid dynamics to study its flowing behavior of themedium.Recently, the present authors constructed a fluid dy-namic model for the dilatant fluid by phenomenologicallyintroducing an internal state variable, which determinesthe viscosity of the medium[30]. The state variable it-self is determined by the local stress[31]. The purposeof this paper is to present detailed study on the flow-ing property of the medium represented by the model.We demonstrate that the model shows the discontinu-ous shear thickening transition and the hysteresis uponchanging the shear rate as has been observed in exper-iments. It is also shown that the steady shear flow be-comes unstable for a certain parameter range against theshear thickening oscillation, where the medium alternatesbetween the thickened and the relaxed states.The paper is organized as follow. The model is in-troduced in Sec.2, and it is examined for a simple uni-form shear flow configuration in Sec.3. Similar analysisis given for the gravitational slope flow and Poiseuilleflow in Sec.4. The response to an impact is simulatedin Sec.5. Effects of inhomogeneity is studied in Sec.6 bytwo-dimensional simulations. Summary and discussionsare given in Sec.7.
II. MODEL
The model is based on the fluid dynamics with an in-ternal state variable that describes the local structureof particles dispersed in the liquid. The viscosity of themedium is determined by the internal state, which inturn changes in response to the local shear stress. Weintroduce each element of the model in the following. a. Fluid dynamics:
The dynamics of the medium asa fluid is represented by the velocity field v ( r ), and isgoverned by the hydrodynamic equation, ρ Dv i Dt = ∂∂x j (cid:16) − P δ i,j + σ i,j (cid:17) + ρg i , (1)where the Lagrange derivative is introduced: DDt ≡ ∂∂t + v j ∂∂x j . (2)The symbols ρ , P , and σ i,j represent the density, thepressure, and the ( i, j ) component of the viscous stress (a) (b) FIG. 1: Schematic pictures for granular configurations: a re-laxed state(a) and a jammed state(b). tensor ˆ σ , respectively. The last term in Eq.(1) representsthe body force on the fluid due to the gravitational accel-eration g i . We employ Einstein’s rule for the summationover repeated suffixes.We consider the incompressible fluid, thus the pressure P is determined by the incompressible condition ∇ · v ( r ) = 0 . (3)The viscous stress tensor is assumed to be expressedthrough the ordinary relation σ i,j = η ( φ ) ˙ γ i,j , (4)with the shear rate tensor˙ γ i,j ≡ ∂v i ∂x j + ∂v j ∂x i − δ i,j ∂v l ∂x l . (5)Note that Eq.(4) does not represent a linear viscositybecause the viscosity η is not constant but depends onthe internal state variable φ of the medium. b. Internal state of the medium: The dilatant fluidcontains dispersed granular particles, which provides thesystem with an internal degree of freedom for a macro-scopic description. Fig.1 shows a schematic illustrationfor a relaxed state(a) and that for a jammed state(b).The internal state may have a vector or even higher or-der symmetry in general, but in this work we study asimple case where the state is represented by a scalarfield φ ( r ). We assign φ = 0 for the relaxed state and φ = 1 for the jammed state.For a given flow field v ( r ), we assume that there existsa stationary value φ ∗ , toward which the state variable φ changes as τ DφDt = − (cid:16) φ − φ ∗ (cid:17) (6)with the time scale τ .We may assume that τ is constant in the case where theinternal state changes due to the thermal fluctuation orsome other mechanism independent of the flowing field.However, we adopt the variable time scale τ that is in-versely proportional to the local shear rate ˙Γ, τ = r ˙Γ − (7)with a dimensionless constant r , because it is more natu-ral to suppose that the state change is driven by the flowdeformation. Note that this form of τ does not introducea new time scale to the system and makes it respondquite peculiarly to an external impact.The stationary value φ ∗ is determined by the local flowand we assume that it is an increasing function of thelocal stress S . We employ a simple form φ ∗ ( S ) = φ M · ( S/S ) S/S ) (8)with the characteristic shear stress S . The parameter φ M represents the value of the state variable in the highstress limit and should depend upon the volume frac-tion of the granules and some other parameters of themedium.For the scalar values of the shear rate ˙Γ and the shearstress S in Eqs.(7) and (8), we adopt the definitions˙Γ ≡ r
12 Tr[ˆ˙ γ · ˆ˙ γ ] , S ≡ r
12 Tr[ˆ σ · ˆ σ ] , (9)which reduce to the ordinary shear rate and shear stressin the case of simple shear flow. c. Viscosity: The shear thickening property of themodel comes from the φ -dependence of the viscosity, forwhich we assume η ( φ ) = η exp (cid:20) A φ − φ (cid:21) (10)with the viscosity in the relaxed state η and a dimension-less parameter A . We have introduced the Vogel-Fulchertype strong divergence at the jamming point φ = 1 in or-der to represent severe thickening observed in the dilatantfluid. In Eq.(10), the state variable φ plays an analogousrole with the inverse temperature in the glass transition.Note that the state variable φ cannot be φ > φ M > ց φ ր d. Unit system: For numerical presentation, we em-ploy the unit system where η = S = ρ = 1 , (11)namely, the time, length, and mass are measured by theunits τ ≡ η S , ℓ ≡ r η ρ τ , m ≡ ρ ℓ , (12)respectively. The rate 1 /τ gives the scale for the shearrate where thickening occurs, and the length scale ℓ isthe corresponding hydrodynamic length scale. For thecornstarch suspension of 41 wt%[3], these parametersmay be estimated as S ≈
50 Pa, η ≈
10 Pa · s, and ρ ≈ kg / m , which give τ ≈ . ℓ ≈ θ g (d) ∆ xL (c) (b)(a) P xz h S e S e xz h−h h z x R FIG. 2: Simple flow configurations and the coordinate system:(a) shear flow, (b) gravitational slope flow, (c) Poiseuille flow,and (d) impact by a bullet.
III. SIMPLE SHEAR FLOW UNDEREXTERNAL SHEAR STRESS
First, we will study behaviors of the dilatant fluid fora simple shear flow under an externally applied shearstress(Fig. 2(a)). The velocity field is assumed to be v = ( u ( z, t ) , ,
0) and the external stress imposes theboundary condition S ( z, t ) (cid:12)(cid:12)(cid:12) z = ± h = S e , (13)where we have introduced the notation for the shearstress S ( z, t ) ≡ η ( φ ) ˙ γ ( z, t ) (14)with the shear rate ˙ γ ( z, t ) ≡ ∂u ( z, t ) ∂z . (15) h is the half width of the flow and S e is the applied stressat the boundaries(Fig. 2(a)). Then, Eqs.(1) and (6)become ρ ∂u ( z, t ) ∂t = ∂∂z S ( z, t ) , (16) r ∂φ ( z, t ) ∂t = −| ˙ γ ( z, t ) | (cid:16) φ ( z, t ) − φ ∗ (cid:0) S ( z, t ) (cid:1)(cid:17) , (17)In the following, we first examine a steady flow solu-tion, then perform the stability analysis for the solutionand the numerical simulation for these equations of mo-tion. S e γ. ∗ φ M = 0.75= 0.85= 1.0 = 2.0 −1 −2 −1 S e γ. ∗ FIG. 3: (Color online) The stress-shear rate relation for theviscosity given by Eq.(18) for various φ M with A = 1. Theinset shows the plots in the logarithmic scale.
1. Steady flow solution
The steady solution for Eqs.(13) ∼ (17) can be readilyobtained as φ = φ ∗ ( S e ) , ˙ γ = S e η (cid:0) φ ∗ ( S e ) (cid:1) ≡ ˙ γ ∗ ( S e ) . (18)From these equations, we can obtain the relationshipbetween the stress and the shear rate, which is plottedin Fig.3 for various values of φ M with A = 1. In thelogarithmic plots, the straight line with the slope 1 cor-respond to the linear stress-shear rate relation with aconstant differential viscosity. One can see that there aretwo regimes: the low viscosity regime in the low shearstress and the high viscosity regime in the high shearstress. Between the two regimes, there is a branch wherethe shear rate decreases for increasing shear stress. Thestate in the middle branch can be unstable against in-finitesimal perturbation.From this stress-shear rate relation, we expect thereshould be hysteresis upon changing the shear rate; if thesystem starts from the low shear rate on the lower branch,the stress increases continuously, but before the systemreaches the end of the lower branch, it should jump to theupper branch by discontinuous increase of the stress. Ifthe system starts from the high shear rate on the upperbranch and the shear rate decreases, the stress shouldjump to the lower branch before the system reaches theend of the upper branch. This sudden increase/decreaseof stress corresponds to the discontinuous change of vis-cosity in the shear thickening.
2. Linear stability of the steady flow
Now, we examine the linear stability of the steadyshear flow given by Eq.(18). For a full analysis, an arbi-trary perturbation should be allowed, but here we exam-ine the linear stability against the restricted perturbationwhere the velocity is in the x -direction and the spatial dependence is only on the z coordinate: v ( r , t ) = ( ˙ γ ∗ z + δu ( z, t ) , , , φ ( r , t ) = φ ∗ ( S e )+ δφ ( z, t ) , (19)then the dynamics is analyzed using Eqs.(16) and (17).Even within this restriction, we will see the steady shearflow in the middle branch may become unstable and theoscillatory flow arises.The linearized equations for the perturbation are nowgiven by ρ ∂∂t δ ˙ γ ( z, t ) = η ∗ ∂ ∂z δ ˙ γ ( z, t ) + η ′∗ ˙ γ ∗ ∂ ∂z δφ ( z, t ) , (20) r ∂∂t δφ ( z, t ) = ˙ γ ∗ (cid:16) φ ′∗ η ∗ δ ˙ γ ( z, t ) + (cid:0) − φ ′∗ η ′∗ ˙ γ ∗ (cid:1) δφ ( z, t ) (cid:17) , (21)where the primes denote the derivative by its argument,and we have introduced the abbreviated notations, η ∗ ≡ η (cid:0) φ ∗ ( S e ) (cid:1) , η ′∗ ≡ dη ( φ ) dφ (cid:12)(cid:12)(cid:12)(cid:12) φ = φ ∗ ( S e ) , φ ′∗ ≡ dφ ∗ ( S e ) dS e . (22)Then, the growth rate Ω k of the perturbation for theFourier component with the wave number k in the z di-rection is determined by (cid:12)(cid:12)(cid:12)(cid:12) ρ Ω k + k η ∗ , k η ′∗ ˙ γ ∗ − ˙ γ ∗ φ ′∗ η ∗ , r Ω k − ˙ γ ∗ (cid:0) − φ ′∗ η ′∗ ˙ γ ∗ (cid:1) (cid:12)(cid:12)(cid:12)(cid:12) = 0 . (23)This gives a positive real part of Ω k for the wave num-ber k that satisfies0 < k < k c ≡ r (cid:18) ρη ∗ (cid:19) S e (cid:18) − d ˙ γ ∗ dS e (cid:19) (24)in the case d ˙ γ ∗ /dS e <
0, i.e. S e is in the unstable branchof the shear stress-shear rate curve. Since the smallestpossible wave number k for the perturbation is π/ (2 h )and η ∗ /ρ is the kinematic viscosity for the steady flow,we can interpret this result in the way that the steadyshear flow in the unstable branch is unstable as long asthe system width is larger than the momentum diffusionlength due to the viscosity. For a given external shear stress S e in the unstablebranch, the flow becomes unstable for the system widerthan 2 h c ≡ π/k c , where the growth rate Ω k has a finiteimaginary part ω c given by ω c ≡ s S e rρ k c = 1 r s S e (cid:18) − d ˙ γ ∗ dS e (cid:19) ˙ γ ∗ . (25)Note that the scales of k c and ω c are typically set by 1 /ℓ and 1 /τ although their actual values depend on r andthe other system parameters, i.e. A and φ M .
3. Shear thickening oscillation in unstable shear flow
The oscillatory behavior of the shear flow in the un-stable regime can be seen by numerically integrating A v e r age S hea r R a t e Time A =1, φ M =1, r =0.1, S e =1.1 h =1.32.03.0 FIG. 4: (Color online) Oscillation of the average shear rate u ( h ) /h in the shear flow for h = 1 .
3, 2, and 3 with A = φ M =1, r = 0 .
1, and S e = 1 .
1. The (green) line, that overlaps theplot for h = 1 .
3, shows the plot for f ( t ) = c + c e − t/τ sin( ωt + θ ) with ω = 4 . τ = 5 . θ = 0 . c = 0 .
33, and c = 0 . Eqs.(16) and (17) with Eqs.(14) and (15). In Fig.4, theaverage shear rates u ( h ) /h for (anti-)symmetric solutionsare plotted as a function of time for various system width h with the constant shear stress S e = 1 . A = φ M = 1 and r = 0 .
1. The initial state isprepared as the steady solution (18) for S e = 1. For thisset of parameters, k c = 1 .
18, which gives h c = 1 .
33 and ω c = 3 . h = 1 .
3, which is smaller than h c , the flow showsoverdumped sinusoidal oscillation with the angular fre-quency 4.0, which is close to ω c . For larger h , the oscilla-tion becomes self-sustained and non-linear; the gradualbuildups of the flow speed are followed by sudden drops.This non-linear oscillation of shear thickening fluid isshown in more detail in Fig.5, where the time develop-ment of the φ ( z ) and u ( z ) are plotted. Only the positivehalf of the solution is plotted for the (anti-)symmetricsolution. In the plots, the oscillation starts from the al-most uniform shear flow in the high viscosity state with alarger value of φ . This flow cannot be completely uniformbecause the external shear stress S e is in the unstablebranch, thus the flow speed builds up gradually as theinternal state φ relaxes to reduce the viscosity, but even-tually, φ starts increasing when the shear stress becomeslarge enough. Then, larger value of φ causes higher vis-cosity, which decelerates the flow speed, but this causeseven higher value of φ because the inertia stress due tothe deceleration is added on the top of the stress by theshear flow, which results in the sudden drop of the flowspeed. IV. GRAVITATIONAL SLOPE FLOW ANDPOISEUILLE FLOW
Similar analyses are performed for a gravitational slopeflow and Poiseuille pipe flow.
Simple Shear Flow t = 6.4
Simple Shear Flow t = 6.4t = 6.6t = 6.6t = 6.8t = 6.8t = 7.0t = 7.0t = 7.2t = 7.2t = 7.4t = 7.4t = 7.6t = 7.6t = 7.8t = 7.8t = 8.0t = 8.0t = 8.2t = 8.2t = 8.4t = 8.4t = 8.6t = 8.6t = 8.8t = 8.8t = 9.0t = 9.0
0 0.5 1 1.5 2 0 0.5 1 φ z t = 9.0
0 0.5 1 1.5 2 0 0.5 1 uz FIG. 5: (Color online) Time development of φ ( z ) (left) and u ( z ) (right) in the shear flow oscillation for A = φ M = 1, r = 0 . S e = 1 .
1, and h = 2. Only the positive parts of theflow ( z >
0) are presented.
A. Gravitational slope flow
For the gravitational slope flow (Fig.2(b)), Eqs.(1) and(6) should be solved with the boundary conditions v (cid:12)(cid:12) z =0 = 0 , ˆ σ · n (cid:12)(cid:12) z = h ( x,y ) = 0 , (26)where we have assumed that the bottom of the flow islocated at z = 0 and the flow depth at ( x, y ) is given by h ( x, y ); the vector n represents the normal vector to theflow surface. The gravitational body force is given by g = ( g sin θ, , − g cos θ ) ≡ ( g k , , − g ⊥ ) (27)with the slope angle θ .For the flow field v = ( u ( z, t ) , , ρ ∂u ( z, t ) ∂t = ∂S ( z, t ) ∂z + g k , (28)0 = − ∂P ( z ) ∂z − g ⊥ , (29) r ∂φ ( z, t ) ∂t = −| ˙ γ ( z, t ) | (cid:16) φ ( z, t ) − φ ∗ (cid:0) S ( z, t ) (cid:1)(cid:17) (30)with the shear stress (14) and the boundary conditions(26) are given by u (0) = 0 , ∂u ( z, t ) ∂z (cid:12)(cid:12)(cid:12)(cid:12) z = h = 0 . (31)Eq.(29) can be solved immediately to give the pressure P ( z ) = g ⊥ ( h − z ) + P (32)with the atmospheric pressure P .The steady solution for Eqs.(28) and (30) under theboundary condition (31) is given by φ ( z ) = φ ∗ (cid:0) S g ( z ) (cid:1) , ˙ γ ( z ) = S g ( z ) η (cid:16) φ ∗ (cid:0) S g ( z ) (cid:1)(cid:17) (33)with the gravitational shear stress S g ( z ) ≡ g k ( h − z ) . (34) z u Gravitational Slope Flow (a) A =1 φ M =1 h =2 g || = S u r f a c e s peed g || (b) A =1 φ M =1h = 321 0 0.5 1 0 5 10 15 20 S u r f a c e s peed Time (c) A= φ M =1, r =0.1, h =2 g || = 0.60.70.81.0 FIG. 6: (Color online) Steady Gravitational flows: (a) theflow speed profiles as a function of z , (b) the surface flow speedvs g k , and (c) the time development of the surface speed. From these, the flow speed u ( z ) and the flux per unitwidth Φ G can be calculated by u ( z ) = Z z ˙ γ ( z ′ ) dz ′ , Φ G = Z h u ( z ) dz. (35)In Fig.6, the flow speed profiles and the surface speedsgiven by Eq.(35) plotted for several sets of parameters.The depth dependences of the flow speed are shown inFig.6(a) for some values of g k ; For small g k , the flow speeddepends upon the depth parabolically as in a Newtonianfluid, while, for larger g k , the flow speed profile develops aconvex part, which corresponds with the unstable branchof Fig.3 in the shear flow. In Fig.6(b), the surface speed u ( h ) are plotted as a function of g k for some values of h .One can see they decreases for large g k , which means thatthe fluid flows slower for larger inclination angle. This isbecause the viscosity of the fluid becomes large in thehigh shear stress caused by large g k . B. Poiseuille Flow
Pressure driven pipe flow with the cylindrical symme-try around the x -axis (Fig.2(c)) is governed by the equa-tion ρ ∂u ( r, t ) ∂t = + ∆ PL + 1 r ∂∂r (cid:16) rS ( r, t ) (cid:17) , (36) r ∂φ ( r, t ) ∂t = −| ˙ γ ( r, t ) | (cid:16) φ ( r, t ) − φ ∗ (cid:0) S ( r, t ) (cid:1)(cid:17) (37)with the shear stress and the shear rate S ( r, t ) = η ( φ ) ˙ γ ( r, t ) , ˙ γ ( r, t ) = ∂u ( r, t ) ∂r . (38) r u Poisseuille Flow (a) A =1 φ M =1 R =3 ∆ P/L = F l u x ∆ P/L (b) A =1 φ M =1 R = 432 0 5 10 15 0 5 10 15 20 F l u x Time (c) A= φ M =1, r =0.1, R =3 ∆ P/L = 0.781.0 1.2 1.3
FIG. 7: (Color online) Poiseuille flow: (a) the flow speedprofiles as a function of r , (b) the flow flux vs pressure gradient∆ P/L , and (c) the time development of the flux.
Here, ∆ P ( >
0) is the pressure drop along the pipe overthe length L , and r is the distance from the central axis: r ≡ p y + z .The steady flow solution for this configuration is givenby φ ( r ) = φ ∗ ( S P ( r )) , ˙ γ ( r ) = S P ( r ) η (cid:16) φ ∗ (cid:0) S P ( r ) (cid:1)(cid:17) (39)with the Poiseuille shear stress S P ( r ) ≡ −
12 ∆
PL r. (40)In Fig.7, the flow speed profiles u ( r ) and the flow fluxΦ defined as Φ ≡ Z R u ( r )2 πr dr (41)are plotted. General features of the flow is analogous tothose of the gravitational flow, and the flow flux decreasesupon increasing the pressure gradient for the large pres-sure gradient because of the shear thickening. C. Shear thickening oscillation in gravitational flowand Poiseuille flow
These steady flows become unstable when the shearstress is in the range of the unstable branch at some re-gion of the flow. The oscillations in the surface flow speedand the flow flux are plotted for the gravitational andPoiseuille flow in Figs.6(c) and 7(c), respectively. Theshear thickening oscillation appears in a large enoughsystem for a certain range of external drive g k or ∆ P/L ;The system length scale should be larger than the viscouslength scale of the flow, and the external drive should bein the range where some part of the flow is in the unsta-ble branch. From the plots, one can see the oscillationdisappears when the external drive is either too small ortoo large. In the former case, the fluid behaves as New-tonian while in the latter case the size of the unstableregion becomes too small. The shape of the oscillationin the non-linear oscillation regime is saw-teeth like, i.e.gradual increases followed by sudden drops, as we havediscussed in the simple shear flow case.The spatial variation of oscillatory flow are shown inFig.8. The general feature is the same with that of theshear flow, but the value of φ is zero at the surface ofthe gravitational flow and at the center of Poiseuille flowbecause the shear is zero. V. RESPONSE TO AN EXTERNAL IMPACT
One of the peculiar features of the dilatant fluid is in-stantaneous hardening by an external impact. It hardensalmost immediately upon application of an external im-pact and allows little deformation like rigid material. Ithas been demonstrated that the hardening is so rapidthat the material can be used for a body armor to stopa bullet[6]. Such instantaneous hardening cannot be ex-plained by the transformation between steady configu-rations of granules, but must be a result of the failureto rearrange the granular configuration due to some ob-struction. Upon sudden impact, the granules are inhib-ited to rearrange their configurations due to to either dis-sipation by the interstitial fluid or the jamming by directcontacts. In the case of slow deformation, the stress islow and the lubrication due to the fluid allows the gran-ules to re-arrange themselves so that they can pass eachother.In the present model, this aspect of the medium is rep-resented by Eq.(7) that the relaxation rate of the inter-nal state is proportional to the shear rate. For a suddendeformation, the state variable changes to a high stressvalue as the medium deforms; When the medium is dense( φ M &
1) and the external impact is strong enough, thestate reaches the jammed state after a certain amount ofdeformation, which is almost independent of the speed ofdeformation.In order to demonstrate this aspect of the model, weperform simple simulations that the layer of fluid of thethickness h is driven by a sudden motion of the upperboundary wall at z = h with the fixed lower boundary at z = 0 (Fig.1d). Let U ( t ) ≡ u ( h, t ) be the velocity of theupper wall. Initially, the fluid is at rest, u ( z, t ) = 0 , φ ( z, t ) = 0 , U ( t ) = 0 for t < , (42)then the upper wall is moved suddenly by the velocity u at t = 0, U (0) = u . For t >
0, the velocity of the upper wall is determined by m dU ( t ) dt = − η (cid:0) φ ( h, t ) (cid:1) ∂u ( z, t ) ∂z (cid:12)(cid:12)(cid:12) z = h , (43)with Eqs.(16) and (17), where m is the mass of the upperwall per unit length.Fig.9 shows the displacement X ( t ) of the upper wall, X ( t ) = Z t U ( t ′ ) dt ′ , (44)for the three cases, φ M = 0 .
8, 1, and 2 for various ini-tial speeds u increases. The wall decelerates rapidly asthe fluid thickens in response to the stress, and eventu-ally stops. For φ M = 0 .
8, the final wall displacementincreases as the initial speed u . On the other hand, for φ M = 1 and 2, the final displacement hardly depend on u when u >
5. This is because the fluid gets jammedat a certain strain as it deforms, and cannot deforms fur-ther. However, when the initial speed is small enough,the upper wall does not stop quickly because the fluiddoes not thicken as shown Fig.9(c).
VI. TWO DIMENSIONAL INHOMOGENEOUSFLOW
Now, we present the results of numerical simulationsfor two dimensional system in the simple shear configu-ration (Fig.2(a)) in order to examine how the inhomo-geneity in the x direction affects the system behavior,especially in the case of shear thickening oscillation.The velocity field is assumed to be in the x − z plane, v = ( u ( x, z, t ) , , w ( x, z, t )), and in the x direction weemploy the periodic boundary condition with the systemlength L . We take L = 10 h in the present simulations.The fluid dynamic equation (1) is integrated using thestandard MAC(Marker-and-Cell) method[32] for the in-compressible fluid, and | ∇ · v | is kept less than 10 − . Eu-ler method is employed for the time integration of Eqs.(1)and (6).The motion of the plates at z = ± h is controlled sothat the average shear stress on the plate is equal to S e ,1 L Z L η (cid:0) φ ( r , t ) (cid:1) ˙ γ xz ( r , t ) (cid:12)(cid:12)(cid:12) z = ± h dx = S e . (45)As for the initial configuration at t = 0, we assumethat the fluid is at rest and the state variable φ is closeto zero with small fluctuations introduced at every com-putational grid point r i , v ( r ,
0) = 0 , φ ( r i ,
0) = ξ i , (46)where ξ i is a random variable uniformly distributed over[0 , ǫ ) with a small parameter ǫ . We take ǫ = 10 − . Notethat, for the case of ǫ = 0, all quantities do not dependon x , thus the simulations reduce to the one dimensionalcase given in Sec. III. (a) Gravitational Slope Flow t = 10.4t = 10.4t = 10.6t = 10.6t = 10.8t = 10.8t = 11.0t = 11.0t = 11.2t = 11.2t = 11.4t = 11.4t = 11.6t = 11.6t = 11.8t = 11.8t = 12.0t = 12.0t = 12.2t = 12.2t = 12.4t = 12.4t = 12.6t = 12.6t = 12.8t = 12.8t = 13.0t = 13.0
0 0.5 1 1.5 2 0 0.5 1 φ z t = 13.0
0 0.5 1 1.5 2 0 0.5 1 uz (b) Poiseuille Flow t = 15.0t = 15.0t = 15.2t = 15.2t = 15.4t = 15.4t = 15.6t = 15.6t = 15.8t = 15.8t = 16.0t = 16.0t = 16.2t = 16.2t = 16.4t = 16.4t = 16.6t = 16.6t = 16.8t = 16.8t = 17.0t = 17.0t = 17.2t = 17.2t = 17.4t = 17.4t = 17.6t = 17.6
0 1 2 3 0 0.5 1 φ r t = 17.6
0 1 2 3 0 0.5 1 ur FIG. 8: (Color online) Time development of φ ( z ) (left) and u ( z ) (right) in the oscillatory flow of the gravitational slope flow(a) and the Poiseuille flow (b). The parameters are A = φ M = 1 and r = 0 . g k = 0 . h = 2 for the gravitational flowand with ∆ P/L = 1 and R = 3 for Poiseuille flow. D i s p l a c e m en t X Time(a) φ M = 0.8 Time(b) φ M = 1.0 Time(c) φ M = 2.0 u = 4020105 D i s p l a c e m en t X Time(d) φ M = 1.0 u = 521 FIG. 9: (Color online) The time dependence of the displace-ment X after the impact for the system of φ M = 0 . u = 40, 20, 10, and 5,and for the system of φ M = 1 . u = 5, 2, and 1 (d).The other parameters are h = 2, r = 0 .
1, and A = m = 1. A. Flow diagrams
Fig.10 shows a flow diagram in the S e − h plane for φ M = 0 .
85 with A = 1 and r = 0 . φ M =1). The diagrams are determined by the simulations atthe points with marks.In the steady shear flow (grey) and the oscillatory flow(purple) regions, the initial fluctuations do not grow, thusthe flows remain homogeneous in the x direction and arethe same with those in the corresponding one dimen-sional cases in Sec.III (Fig.11); the dashed lines show theboundary for the two regimes in the one-dimensional casegiven by k c ( S e ) = π h (47)using the definition of k c in Eq.(24) as a function of S e .In the low S e side, one can see that this coincides withthe corresponding boundary in the two-dimensional casebetween the steady shear flow and the oscillatory flow.The major difference between the one and the twodimensional cases is that these two homogeneous flowregimes are limited to the smaller S e side. In the larger FIG. 10: (Color online) Flow diagram of the shear flow for φ M = 0 .
85. The internal state variable φ does not depend on x in the region colored with gray and blue. The inset is thesame diagram obtained for φ M = 1. The other parametersare A = 1, r = 0 .
1, with L = 10 h . S e region, the initial fluctuations in the state variable φ grows, thus the flow results in the inhomogeneous flowin the case of φ M = 0 . φ M = 1 . B. Inhomogeneous oscillatory flow
First, we examine the flow for φ M = 0 .
85. In thiscase, the viscosity does not diverge and the medium keepsflowing. In Fig.12, the time evolution of the upper platevelocity U p is plotted along with the case without initialfluctuations. The flow shows irregular oscillation withsmaller amplitude compared with the noiseless case.The snapshots of φ for a single cycle of oscillation inFig.13 reveal that the whole system is not thickened andthe oscillation is governed by a few thickening bands. Atthe time when U p reaches its minimum (a), the thickeningbranch with the high value of φ (the red region) is beingextended along the direction of (1 , U p gradually increases FIG. 11: (Color online) Time development of the upperplate velocity U p in the oscillatory flow regime in the two-dimensional simulation. The initial fluctuations decay quicklyand the flows show homogeneous oscillation as in the case ofthe one-dimensional system.FIG. 12: (Color online) The time evolution of the upper platevelocity U p in the inhomogeneous oscillatory flow regime with φ M = 0 .
85. The other parameters are h = 5, L = 50, S e =1 . A = 1, and r = 0 .
1. The initial fluctuation is given by ǫ = 10 − . The uniform oscillation flow with ǫ = 0 (the dashedline) is shown for comparison. (b), and eventually the branch breaks off and U p reachesmaximum (c). Then, high shear rate makes the brokenbranches extend again to the other side to cause suddendeceleration (d). The thickening branches first appearboth in (1 ,
1) and (1 , −
1) direction, but the latter tendto disappear and transforms into the (1 ,
1) direction inthe course of time, and only the thickening branches inthe (1 ,
1) direction remain.
C. Jamming caused by the instability
For the case of φ M = 1, the viscosity can divergesand the instability of the homogeneous flow causes thejamming to stop the flow.Figs.14 and 15 shows the simulation results for φ M = 1. FIG. 13: (Color online) The snapshots of the state variable φ taken during a cycle of oscillation presented in Fig.12. Thearrows indicate flow velocity.FIG. 14: The spatial variation of viscosity η at z = 0 atseveral times in the jamming regime. The parameters are h = 3, L = 30, S e = 1 . φ M = 1 with A = 1 and r = 0 . The time evolution of the viscosity distribution at z = 0is shown in Fig.14. Initially, the viscosity is rather uni-form with some fluctuations, but peak structure appearssoon around t = 3 with a certain characteristic lengthscale. Some of the peaks grow sharply and the thicken-ing regions strongly localize( t & P (a) andthe state variable φ (b) at t = 4. The fluctuation of φ first stands out just below the moving plates, but higher φ regions form band structure, and extend along the prin-cipal axis of the shear deformation, (1 ,
1) and (1 , − FIG. 15: (Color online) (a) The spatial distribution of thepressure P and (b) the internal state variable φ in the systempresented in Fig.14 at t = 4.FIG. 16: (Color online) The time evolution of (a) the up-per plate velocity U p , and (b) the maximum viscosity in thesystem presented in Fig.14. The dashed lines represent theoscillatory flow without fluctuations. plate, and they jam the system.In Fig.16(a), we present the velocity of upper plate U p as a function of time. The solid line shows the time evo-lution starting from the internal state with fluctuations,and the dotted line represents the case without fluctua-tions for comparison. The state variable suddenly loseshomogeneity at t = 2 .
3, then thickening branches ap-pears, and the velocity U p drops to zero. The maximumviscosity is always found inside the thickening branchfor t & .
3, and its value sharply increases as plottedin Fig.16(b). We cannot simulate the system up to thetime when the plates motion actually stops because thenumerical time integration becomes difficult as the vis-cosity becomes large, since it requires smaller time step.In the present case, however, we expect the system isjammed because the decrease of U p and the increase ofmaximum viscosity are rapid and monotonic. VII. SUMMARY AND DISCUSSIONS
The shear thickening shown by a dense mixture ofgranules and fluid has some peculiar features: (i) instan-taneous hardening, (ii) fast relaxation to flowing state,(iii) rigid thickened state, (iv) hysteretic thickening tran-sition, (v) oscillatory flowing behavior. We constructed a fluid dynamics model by introducing a phenomenolog-ical state variable and showed that the model can de-scribe these features. Especially, we demonstrated that the shear thickening oscillation appears in various shearflow configurations. a. Comparison with visco-elasticity:
A visco-elasticfluid such as a polymer melt shows analogous behavior tothe dilatant fluid; It behaves like a solid in a short timescale and like a fluid in a long time scale. This may becompared with that of the dilatant fluid, i.e. instanta-neous hardening in response to an external impact andfluidization after relaxation of the applied stress. How-ever, there are some important differences; the visco-elastic medium changes its behavior according to the ob-servation time scale, and it also allows large elastic de-formation even in a short time solid like behavior. Onthe other hand, the dilatant fluid changes its behavioraccording to the stress, i.e. it stays hardened while it isunder the stress and starts flowing within a few secondafter the removal of the applied stress. The dilatant fluidallow little deformation even under large stress. b. Shear stress thickening:
In constructing themodel, we assume the fluid is shear-stress thickening , i.e.the viscosity depends upon the state variable φ , and thesteady value of the state variable φ ∗ is determined bythe local stress as in Eq.(8). It is instructive to see whatwould happen if we assume φ ∗ as a function of the shearrate φ ∗ ( ˙ γ ). In this case, the viscosity is directly givenas a function of the shear rate, η (cid:0) φ ∗ ( ˙ γ ) (cid:1) , thus we shouldnot have a discontinuous thickening unless we assume adiscontinuity either in φ ∗ ( ˙ γ ) or in η ( φ ).Experimentally, the most direct evidence for the shear-stress thickening should be obtained by the observationthat the pipe flow flux is not monotonically increasing asa function of the applied pressure gradient. c. State variable: Although we introduced the statevariable φ phenomenologically, we suppose that the vari-able represents a certain microscopic property of themedium, such as contact numbers between grains, associ-ated with the restrictions against local rearrangement ofgranular configuration. The variable could be a vector ora tensor, but we examined the scalar case for simplicity.It is notable that even the scalar state variable producesthe anisotropic stress chain like structure in the systemas we have seen in the two-dimensional simulations. d. Jamming and response to impact: A remarkablefeature of the dilatant fluid is that hardening responseis so instantaneous that the medium can be used for abody armor to stop a bullet[6]. We believe that such aninstantaneous severe hardening cannot be explained bya transformation between two steady states, i.e. fromthe fluid state under low stress to the rigid state underhigh stress. Instead, it must be a result that the rapidrearrangement in the granular configuration is inhibited.There are two possible mechanisms that inhibit the gran-ular rearrangement in the densely packed medium: theReynolds dilatancy and the formation of stress chains.The Reynolds dilatancy can inhibit rearrangement by1the fluid friction because the rearrangement should in-duce the strong interstitial fluid flow through granules inorder to compensate the local volume change caused bythe dilatancy. The stress chains through direct contactsbetween granules can be formed by the application exter-nal stress and prevent granules from being rearranged.In the present model, such an aspect of hardening isrepresented by the fact that the relaxational time scalefor the state change is not constant but proportional tothe shear rate (7). Then, the state variable φ reaches asteady value φ ∗ in a time scale where the strain changesby the amount r . Consequently, the medium with φ M & r by ahard impact. e. Shear thickening oscillation: One of interestingresults of the present model is the shear thickening os-cillation. The steady shear flow is unstable when theflow is in the unstable branch of the shear stress-shearrate curve and is wide enough compared to the diffusionlength scale by the viscosity. In this case, the flow showsthe oscillation between the thick state in the high shearregime and the thin state in the low shear regime. Intwo dimensions, uniform oscillation appears only in thesmaller external shear stress region, but some oscillatorybehavior remains even when inhomogeneity develops inthe flow.The shear thickening oscillation in the homogeneousflow looks similar to the stick-slip motion in a frictionalsystem. However, there are some important differences;the stick-slip motion starts with the sudden accelera-tion caused by the slip weakening resistance under theconstant speed driving through a mechanism with finiterigidity, while the shear thickening oscillation starts withthe sudden deceleration caused by the shear thickeningtransition under the constant force driving.Aradian and Cates have also studied the dynam-ics of shear thickening fluids and found oscillatorybehaviors[33]. They focused, however, on the regimewhere the structural relaxation time is much larger thanthe fluid dynamical time scale, which is appropriate forliquid crystal systems. Consequently, the time scale ofthe oscillation is of the order or larger than the struc-tural relaxation time scale, therefore the dynamics iscompletely dissipative. On the other hand, in the presentcase, the structural relaxation time of the state variableis set by the shear rate and always comparable with thefluid dynamical time scale, thus the period of the oscilla-tion is determined by the fluid dynamics. It is also clearthat the stress from the inertia plays an important rolein the thickening phase in the oscillation. f. Experiments:
As for the hysteresis, Deegan ob-served the hysteresis loop of the viscosity under the os-cillatory stress for the cornstarch-fluid system and con-sidered it to be the mechanism for persistent holes[34].His system is thinner and in the regime of mild shearthickening in comparison with the system studied by Fallet al[3]. In the present work, we try to model rather se-vere shear thickening in choosing the functional form of −6−4−2 0 2 4 6−1.5 −1 −0.5 0 0.5 1 1.5 S t r e ss Shear rate τ = r /γ. , r = 2φ M = 0.6ω = 0.2 S e = 54321 FIG. 17: (Color online) Hysteresis loops in the shear stressvs the shear rate in response to oscillatory shear stress forvarious amplitudes. The oscillatory stress S ( t ) = S e sin( ωt )is applied to the upper plate located at z = h with the lowerplate fixed at z = 0. The system parameters are h = 0 . φ M = 0 . r = 2, and A = 1. the viscosity Eq.(10), but the present model qualitativelyreproduces main feature of his data, by adjusting the pa-rameters for φ M and r (Fig.17 compared with Fig.5 of[34]). For this set of the parameters, the initial noise de-cays, thus the one and two dimensional simulations givethe same results.Regarding the shear thickening oscillation, we couldnot find any literature on the experiment which showsclear oscillation. This may partly because the systemneeds to be large enough, typically wider than ℓ or afew centimeters, and the inhomogeneity may develops inthree dimensional system, which can obscure the oscilla-tion. The noisy fluctuation often observed in rheometerexperiments under constant stress may be explained bythis mechanism[7, 8]. Nevertheless, one may easily noticethe oscillation around 10 Hz simply by pouring the densewater-starch mixture out of a container. We believe thatthis oscillation should be explained by the shear thicken-ing oscillation. We are planning experiments that allowquantitative comparison with our results.Although in the different physical context, the clearoscillatory flows[35] along with a discontinuous transi-tion and hysteresis[36, 37] have been observed in the liq-uid crystal system that shows the shear thinning due tothe state dependent viscosity. Such behavior could bealso described using the phenomenological model like thepresent one.Chaotic dynamics has been observed in dilute aqueoussolutions of a surfactant in the experiment under the con-stant shear rate in the shear thickening regime[38, 39]. Itwas interpreted as the stick-slip transition between thetwo states of the fluid structure, thus physical relevanceto the present instability is not clear, but we also foundthe chaotic dynamics in the present model in the caseof the constant relaxation time τ with the large systemwidth.2 Acknowledgments
This work is supported by KAKENHI(21540418)(H.N.) and the Danish Council for Independent Research, Natural Sciences (FNU) (N.M.). [1] F. S. Merkt, R. D. Deegan, D. I. Goldman, E. C. Rericha,and H. L. Swinney, Phys. Rev. Lett. , 184501 (2004).[2] H. Ebata, S. Tatsumi, and M. Sano, Phys. Rev. E ,066308 (2009).[3] A. Fall, N. Huang, F. Bertrand, G. Ovarlez, and D. Bonn,Phys. Rev. Lett. , 018301 (2008).[4] H. Freundlich and F. Juliusburger, Trans. Faraday Soc. , 920 (1935).[5] O. Reynolds, Phil. Mag. S5 , 469 (1885).[6] N. J. Wagner and J. F. Brady, Physics Today , Oct.27 (2009).[7] H. Laun, R. Bung, and F. Schmidt, J. Rheol. , 999(1991).[8] D. Lootens, H. van Damme, Y. H´emar, and P. H´ebraud,Phys. Rev. Lett. , 268302 (2005).[9] R. Hoffman, Trans. Soc. Rheol. , 155 (1972).[10] R. Hoffman, J. Colloid and Interface Sci. , 491 (1974).[11] R. Hoffman, J. Rheol. , 111 (1998).[12] H. Barnes, J. Rheology , 329 (1989).[13] J. J. Erpenbeck, Phys. Rev. Lett. , 1333 (1984).[14] M. Chow and C. Zukoski, J. Rheol. , 33 (1995).[15] H. Laun, R. Bung, S. Hess, W. Loose, O. Hess, K. Hahn,E. H¨adicke, R. Hingmann, and F. Schmidt, J. Rheol. ,743 (1992).[16] J. Bender and N. J. Wagner, J. Rheol. , 899 (1996).[17] B. J. Maranzano and N. J. Wagner, J. Chem. Phys. ,10291 (2002).[18] R. G. Egres and N. J. Wagner, J. Rheol. , 719 (2005).[19] J. Brady and G. Bossis, J. Fluid Mech. , 105 (1985).[20] J. Melrose and R. Ball, J. Rheol. , 937 (2004).[21] X. Cheng, J. H. McCoy, J. N. Israelachvili, and I. Cohen,Science , 1276 (2011).[22] E. Brown and H. M. Jaeger, Phys. Rev. Lett. , 086001(2009).[23] J. Melrose and R. Ball, J. Rheol. , 961 (2004).[24] R. S. Farr, J. R. Melrose, and R. C. Ball, Phys. Rev. E , 7203 (1997). [25] M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, andP. Claudin, Phys. Rev. Lett. , 1841 (1998).[26] E. Bertrand, J. Bibette, and V. Schmitt, Phys. Rev. E , 060401 (2002).[27] E. Brown and H. M. Jaeger, arXiv: [cond-mat.soft] ,1010.4921v1 (2010).[28] D. A. Head, A. Ajdari, and M. E. Cates, Phys. Rev. E , 061509 (2001).[29] C. B. Holmes, M. E. Cates, M. Fuchs, and P. Sollich, J.Rheol. , 237 (2005).[30] H. Nakanishi and M. Namiko, J. Phys. Soc. Jpn. ,033801 (2011).[31] A similar state variable has been discussed for a gran-ular system in connection with the friction coefficientstrengthening[40–42].[32] F. H. Harlow and J. E. Welch, Phys. Fluids , 2182(1965).[33] A. Aradian and M. E. Cates, Phys. Rev. E , 041508(2006).[34] R. D. Deegan, Phys. Rev. E , 036319 (2010).[35] A. S. Wunenburger, A. Colin, J. Leng, A. Arn´eodo, andD. Roux, Phys. Rev. Lett. , 1374 (2001).[36] D. Bonn, J. Meunier, O. Greffier, A. Al-Kahwaji, andH. Kellay, Phys. Rev. E , 2115 (1998).[37] O. Volkova, S. Cutillas, and G. Bossis, Phys. Rev. Lett. , 233 (1999).[38] R. Bandyopadhyay, G. Basappa, and A. K. Sood, Phys.Rev. Lett. , 2022 (2000).[39] R. Bandyopadhyay and A. K. Sood, EPL (EurophysicsLetters) , 447 (2001).[40] W. Losert, J.-C. G´eminard, S. Nasuno, and J. P. Gollub,Phys. Rev. E , 4060 (2000).[41] M. Lubert and A. de Ryck, Phys. Rev. E , 021502(2001).[42] W. P. Vellinga and C. P. Hendriks, Phys. Rev. E63