Despite the recent interest in the discontinuous shear-thickening (DST) behavior, few computational works tackle the rich hydrodynamics of these fluids. In this work, we present the first implementation of a microstructural DST model in smoothed particle hydrodynamic (SPH) simulation. The scalar model was implemented in an SPH scheme and tested in two flow geometries. Three distinct ratios of local to non-local microstructural effects were probed: zero, moderate, and strong non-locality. Strong and moderate cases yielded excellent agreement with flow curves constructed via the Wyart–Cates (WC) model, with the moderate case exhibiting banding patterns. We demonstrate that a local model is prone to a stress-splitting instability, resulting in discontinuous stress fields and poor agreement with the WC model. The mechanism of stress splitting has been explored and contextualized by the interaction of local microstructure evolution and the stress-control scheme. Analytic solutions for a body-force-driven DST channel flow have been derived and used to validate the SPH simulations with excellent agreement in velocity profiles. Simulations carried out at increasing driving forces exhibited a decrease in flow. We showed that even the simple scalar model can capture some of the key properties of DST materials, laying the foundation for further SPH study of instabilities and pattern formation.

Suspensions of both colloidal1 and non-colloidal2,3 particles are ubiquitous in industrially relevant flows,4 and their rich properties have been subject of extensive study. Some dense suspensions of fine non-Brownian particles can exhibit discontinuous shear thickening (DST) behavior,5 which has been a subject of great interest in the last decade. With possible implications in industries spanning from manufacturing, environmental engineering, and military application, understanding the underlying mechanism responsible for the unique behavior has become paramount. Multiple physical explanations have been proposed,6 with the most widely adopted being the mechanism of Seto et al.7 or Wyart–Cates (WC) model,8 where the particles are held apart in a frictionless state by a characteristic repulsive force. Upon application of external stress, the interparticle repulsion is overcome, bringing about a transition to a frictional state where the developed contact network inhibits the ease of motion of the suspension. Recently, significant effort has been dedicated to characterizing and understanding the effect of the force chain networks in DST.9–12 In the work of Sedes et al.,13 the macroscopic microstructure parameter f, which in the WC model represents the fraction of particles in contact within some averaging fluid element, was shown to be congruent with the microscopic description of contact force networks.

A defining feature of DST, as opposed to continuous shear thickening (CST), is the characteristic “S” shape of the flow curve. This unique non-monotonic shape facilitates unstable behaviors in rheometric14,15 and complex flows. For example, Texier et al.16,17 demonstrated that the negative gradients of the flow curve ( d γ ̇ d σ < 0) present in DST materials cause instability in free surface flows of cornstarch suspension, leading to roll-waves. In the context of simple shear, Cho et al.18 investigated relaxation dynamics in a parallel-plate setup and proposed a non-trivial two-step mechanism based on a viscoelastic response and force network disintegration, although their system never fully transitioned back to the entirely frictionless state. Rathee et al.19 employed boundary stress microscopy (BSM) to study stress distribution in rheometric flows. Their results revealed a high degree of heterogeneity in stress fields,20 including two separate high-stress regions in contact with upper and lower plates, respectively, and propagation of high-stress fronts.21 

WC-type models are inherently microstructural—the key parameter governing the state of the suspension at any point is the arrangement of the particle contact network. Typically, works in the context of WC theory propose a relationship between development of the microstructure and macro-hydrodynamic flow properties through a transient evolution equation. Nakanishi et al.22 employed a scalar model with an assumed S-shape rheology curve in a range of geometries, including simple shear, Poiseuille flow, and free surface inclined plane. The results in simple shear exhibited periodic solutions, accompanied by formation of inhomogeneous fields of the microstructural parameter across the gap. Baumgarten and Kamrin23 proposed a scalar model based on the balance of strain-dependent hardening and softening attributed to load buckling, shearing, and electrostatic repulsion. Their model predicted rheological behavior in good agreement with the experimental evidence.

Another development comes from the work of Gillissen et al.,24 who proposed a full tensorial evolution equation. In addition to predicting DST rheology, they were able to predict negligible first normal stress difference, significant negative second normal stress difference, and flow behavior in more complicated flow arrangements, including flow reversal and superposed transverse oscillations,25 although the values of normal stress differences predicted by their model were quantitatively too high (N1) and too low (N2) relative to the discrete element method (DEM) values, along with the increases during thickening being underestimated. Their model also did not consider the role of contact friction in the development of microstructure and excluded tangential forces and torque balance.

So far, nearly all continuum microstructural models have been implemented in Eulerian frameworks. In recent times, there has been much interest in simulations of complex fluids with grid-less methods, primarily smoothed particle hydrodynamics (SPH).26,27 SPH was originally developed by Gingold and Monaghan28 and Lucy29 in 1977 in the study of astrophysical problems. It is a fully Lagrangian scheme, where the fluid is represented by a set of discrete “particles” carried by flow. Hydrodynamic properties are computed at the position of each particle by averaging via a smoothing kernel at each time step, and the trajectories of particles evolve over time. In comparison with Eulerian and mixed Lagrangian–Eulerian methods, SPH has a number of advantages and disadvantages. Generally, grid-based schemes benefit from superior convergence rates and easier implementation of boundary conditions.26 On the other hand, SPH is inherently conservative and handles free surface boundaries and interfacial multiphase flows naturally30—which is of importance in a number of studies, including the behavior of granular bed loads,31 mixing,32 and cavitating flows.33 Moreover, SPH allows for direct access to the flow history—a crucial feature for microstructural models of complex fluids. This can be leveraged in heterogeneous multi-scale modeling (HMM),34 with an example application in biologically relevant flows,35 or by directly solving integral constitutive equations.36 

Beyond Newtonian fluids, SPH has been used in a host of fields, including solid and complex fluid mechanics. In 2002, Ellero et al.37 simulated viscoelastic flow by incorporating the corotational Maxwell model in SPH, followed by Oldroyd-B.38 Since then, Fang et al.39,40 studied the behavior of free surface flows of an Oldroyd-B fluid, with a focus on addressing tensile instability. Rafiee et al.41 implemented both Maxwell and Oldroyd-B models in SPH free-surface flows, including impacting drop and jet buckling. Xu et al.42 simulated three-dimensional (3D) injection molding of a cross-model fluid. Implementation of Oldroyd-B was studied in a periodic array of cylinders by Vázquez-Quesada and Ellero43 and Grilli et al.,44 and in extrudate swell by Xu and Deng.45 

Furthermore, a range of viscoplastic and elasto-viscoplastic models have been implemented in SPH. Rodriguez-Paz and Bonet46 simulated debris flow down an incline using the Bingham model and the generalized viscoplastic model. Minatti and Paris47 considered free surface granular flows, with a constitutive relation based on the works of Pouliquen et al.48 and Jop et al.,49 and validated their model against granular column collapse experiments. Furthermore, Bingham-like viscoplastic models have been applied to sedimenting flows,50,51 rheometric flows,52 and mixed fluid–structure cases.53,54 Recently, Rossi et al.55 implemented a modified microstructural Papanastasiou model to study thixo-viscoplastic flows around a cylinder. The first application of SPH to DST flow was carried out in the work of Vázquez-Quesada et al.,56 wherein they considered a non-microstructural inverse biviscous model in a simple planar flow.

In this work, we present in detail the first implementation of a scalar non-local microstructural DST rheological model in an SPH scheme. Unlike previous models,56 we obtain typical S-shaped DST rheology in the simple shear and characteristically inflexed velocity profiles accompanied by flow reduction behavior in channel flow. In both cases, stress-splitting is observed as a result of the local term in the evolution equation acting on a perturbed microstructure field. Section II lays out the details of the model used in this study. Section III presents the details of the numerical scheme and boundary condition treatment.

In Sec. IV A, the origin of stress-splitting and the effect of non-locality on establishing spatial correlations in the structure are discussed in context of simple shear. In Sec. IV B, we derive a solution for planar channel flow for a fluid following the Wyart–Cates model and use it to validate our SPH simulations. Finally, we compare our results to a previous non-microstructural model.

In this section, we outline the main equations used to model our DST fluid system. The continuum system can be written in the Lagrangian frame as follows:
d ρ d t = ρ · v ,
(1)
d v d t = 1 ρ · σ + F ,
(2)
σ = p I + η ( v + v T ) ,
(3)
where d / d t / t + v · is the material derivative, ρ ( x , t ) is the density field, v ( x , t ) is the velocity field, σ ( x , t ) is the stress field, p is the pressure, and F is a body force acting on fluid. In this DST model, the local suspension viscosity is determined using the Maron–Pierce expression,57 
η ( ϕ , ϕ J , t ) = η s ( 1 ϕ ϕ J ) 2 ,
(4)
where η ( ϕ , ϕ J , t ) and η s are, respectively, the suspension and constant solvent viscosities. ϕ is the particle volume fraction, and ϕ J is the jamming volume fraction. Assuming particle migration is negligible for the problems studied in the present work, the value of ϕ is kept constant. The value of ϕ J is inferred dynamically as
ϕ J ( f ( x , t ) ) = ϕ m f ( x , t ) + ϕ 0 { 1 f ( x , t ) } ,
(5)
where ϕ m and ϕ 0 are frictional and frictionless divergence volume fractions, respectively.

The discontinuous shear-thickening phenomenon wherein the viscosity of the medium undergoes a sharp transition from a lower to a higher value is modeled with the help of a single microstructure scalar field called the fraction of frictional contacts f ( x , t ). The use of this scalar field model follows from the Wyart–Cates model,8 in which the value of f is determined as exp ( σ / σ ). The repulsive force between particles, F0, sets the characteristic stress scale at the onset of shear thickening σ = F 0 / 6 π a 2.11  Figure 1 shows (a) typical flow curves produced by the model and (b) the corresponding viscosity dependence on shear stress. For sufficiently low volume fractions ( ϕ < 0.50), continuous shear thickening (CST) is achieved. Increasing volume fraction above 0.50 leads to a transition from CST to a DST regime, marked by the appearance of the characteristic S-shape, where stress is a multi-valued function of shear rate. Further increases in the volume fraction sharpen the S-curve up to a point of divergence ( ϕ = ϕ m).

FIG. 1.

WC model of the DST fluids: (a) flow curves in the stress-shear plane and (b) viscosity dependence on shear stress. Plots for σ = 0.005, ϕ m = 0.562, and ϕ 0 = 0.693.

FIG. 1.

WC model of the DST fluids: (a) flow curves in the stress-shear plane and (b) viscosity dependence on shear stress. Plots for σ = 0.005, ϕ m = 0.562, and ϕ 0 = 0.693.

Close modal
The value of f ( x , t ) at a point in the flow field depends on the local value of stress, and the dependency follows a sigmoid function. Time dependence of microstructure is captured in the microstructure evolution equation (6), which in the present work consists of a local and a non-local term. The local term has been widely used in the DST literature17,22,23,58,59 and represents approach of the system toward a shear-stress dependent equilibrium state prescribed by the Wyart–Cates model (7) with a characteristic timescale set by the shear-rate and a microstructural rate constant Kf. In addition to the usual local term, we introduce a non-local contribution in the dynamics of structure, where the microstructural parameter f diffuses with a diffusion coefficient α. In this work, we follow the approach of Kamrin and Henann60 and include microstructure diffusivity directly in the evolution equation. The non-local feature of our model has the benefit of stabilizing the local-driven instability discussed in Sec. IV A,
d f ( x , t ) d t = K f γ ̇ ( f ̂ ( x , t ) f ( x , t ) ) + α 2 f ( x , t ) ,
(6)
f ̂ ( x , t ) = exp ( [ σ σ x y ( x , t ) ] β ) ,
(7)
where γ ̇ ( x , t ) = γ ̇ : γ ̇ / 2 is the local shear calculated as the second invariant of the strain rate tensor γ ̇ ( x , t ) = ( v + v T ) / 2. The first term on the right-hand side of Eq. (6) accounts for the local evolution of microstructure, whereas the second term introduces non-local effects via microstructure diffusivity, where α is the diffusion coefficient.
The Navier–Stokes equations are solved using the SPH methodology. A possible discretized form of the governing equations is given as follows:61 
d ρ i d t = j m j v i j · i W i j ,
(8)
d v i d t = j m j p i + p j ρ i ρ j i W i j + j m j ( η i + η j ) x i j · i W i j ρ i ρ j r i j 2 v i j + g i ,
(9)
where mj is the mass of a particle, W i j = W ( | x i j | ) is the smoothing kernel, and v i j = v i v j , x i j = x i x j, and r i j = x i j · x i j.
Following the weakly compressible SPH (WCSPH) formulation, a 1% variation in density is permitted. The fluid pressure is then calculated based on the Tait's equation of state,
p = ρ 0 c 2 γ [ ( ρ ρ 0 ) γ 1 ] + p b ,
(10)
where ρ0 is the density of the fluid, p b is the background pressure, and c is the numerical speed of sound typically set as 10 times the maximum flow velocity v max.
So far, we have described details of SPH for general flows with Navier–Stokes equations. In this section, we discuss the closure SPH model for the definition of viscosity η i ( f i , t ) in the SPH–DST model. The constitutive equations are written as follows, where the subscript i indicates the particle index (that is f ( x i , t ) = f i),
f ̇ i = K f γ ̇ i ( f ̂ i f i ) + α 2 f i ,
(11)
ϕ i J = ϕ m f i + ϕ 0 ( 1 f i ) ,
(12)
η i = η s ( 1 ϕ i ϕ i J ) 2 .
(13)
In Eq. (11), fi is calculated to evaluate ϕ J in Eq. (12), which in turn is used to obtain viscosity via Eq. (13). The value of γ ̇ in Eq. (11) is computed as the contraction of the strain rate tensor γ ̇ as follows:
γ ̇ = 1 2 γ ̇ : γ ̇ , γ ̇ l k = 1 2 ( u l x k + u k x l ) .
(14)
The SPH approximation of the gradient of the α-component of the velocity vector in β-direction is obtained as
( u l x k ) i = j m j ρ j ( u j l u i l ) i k W i j .
(15)
The non-local term in Eq. (11) is discretized in a manner consistent with the Morris formulation of the moment equation,26 
( 2 f ) i f = j f 2 m j ( f i f j ) x i j · i W i j ρ j r i j 2 .
(16)
To simulate the S-curve (see Fig. 1) with the above-mentioned definition of order parameter f, the choice of ϕ has to be sufficiently close to ϕ m which in the above-mentioned case is ϕ > 0.50. However, the value of η / η s for ϕ = 0.54 at a high shear rate is approximately 3000, which drastically lowers the time step size required for stable integration of the governing equations. A simple idea used to address this was to modify the definition of f so that we can anticipate the S at a lower volume fraction. Therefore, the modified definition of f is as follows:
f ̂ i = exp { ( σ / σ i ) β } ,
(17)
where β is the parameter setting the “sharpness” of the transition between the frictional and frictionless plateaus, without changing the plateau values. We set β = 2 in this work. Increasing it allows us to introduce the previously discussed S-shaped rheology to an otherwise CST volume fraction, resulting in DST behavior without increasing viscosity to computationally prohibitive values. Furthermore, we choose the rest of the parameters, primarily ϕ m and ϕ 0, such that all simulations lie within the prescribed region of 10 3 < R e < 1. The lower bound was chosen to facilitate feasible run time by preventing excessively small time-steps (36) due to shear thickening. The upper bound is used to exclude fluid inertial effects from all cases in this work. Once the values of fi and ϕ J ( f i ) are obtained, the dynamic viscosity of a SPH particle is updated using Eq. (13).

The dummy particle method of Adami et al.62 is employed for imposing no-slip velocity boundary conditions and impermeability conditions. For a full explanation behind this method, we direct the reader to their work. Below we present their method in notation consistent with the present model.

Pressure, density, and velocity are assigned to the dummy particles as
p i w = j f p j W i j + ( g a i w ) · j f ρ j x i j W i j j f W i j ,
(18)
a i w = p i f ρ i f + g ,
(19)
ρ i w = ρ 0 ( p i w p b p 0 + 1 ) 1 γ ,
(20)
p 0 = ρ 0 c 2 γ ,
(21)
v i w = 2 v wall v ̃ i w ,
(22)
where subscript w and f denote the dummy wall and fluid particles, respectively, v i w is the prescribed wall velocity, and
v ̃ i w = j f v j W i j j f W i j .
(23)

A semi-implicit predictor–corrector type integration scheme is used for time marching as follows.

  1. Predictor step:
    v i n + 1 2 = v i n + Δ t 2 ( d v i d t ) n ,
    (24)
    ρ i n + 1 2 = ρ i n + Δ t 2 ( d ρ i d t ) n ,
    (25)
    x i n + 1 2 = x i n + Δ t 2 ( v i ) n ,
    (26)
    f i n + 1 2 = f i n + Δ t 2 ( d f i d t ) n .
    (27)
  2. Corrector step:
    v i n + 1 2 = v i n + Δ t 2 ( d v i d t ) n + 1 2 ,
    (28)
    ρ i n + 1 2 = ρ i n + Δ t 2 ( d ρ i d t ) n + 1 2 ,
    (29)
    x i n + 1 2 = x i n + Δ t 2 ( v i ) n + 1 2 ,
    (30)
    f i n + 1 2 = f i n + Δ t 2 ( d f i d t ) n + 1 2 .
    (31)

The flow variables at the subsequent time step are obtained as follows:
v i n + 1 = 2 v i n + 1 2 v i n ,
(32)
ρ i n + 1 = 2 ρ i n + 1 2 ρ i n ,
(33)
x i n + 1 = 2 x i n + 1 2 x i n ,
(34)
f i n + 1 = 2 f i n + 1 2 f i n ,
(35)
where n represents the current time instant, n + 1 2 represents the predicted time instant (variables at an intermediate time level), and n + 1 represents the corrected time instant. This scheme is second-order accurate in time. Finally, the time step size Δ t is determined based on the following condition:
Δ t = min { 0.25 h sl c , 0.125 h sl 2 ν , 0.25 ( h sl | g | ) 1 / 2 } ,
(36)
where h sl is the smoothing length and ν is the kinematic viscosity. The adaptive time step approach ensures that the CFL condition63 is satisfied, along with additional constraints due to the viscous diffusion and particle acceleration.61 Due to the extreme shear-thickening nature of the material studied in the present work, viscous diffusion is likely to be the most restrictive condition (Fig. 2).
FIG. 2.

Comparison of WC model flow curves, i.e., (a) rate dependence of viscosity and (b) stress dependence of viscosity for β = 1 and β = 2.

FIG. 2.

Comparison of WC model flow curves, i.e., (a) rate dependence of viscosity and (b) stress dependence of viscosity for β = 1 and β = 2.

Close modal

The behavior of the previously introduced DST model is explored by simulating the rheology in a simple shear geometry. The fluid is bound between two solid walls vertically and two periodic boundaries horizontally at a distance l = 0.01 m. The domain is planar with no depth in the vorticity direction. Unlike previous shear-imposed simulations,64 input stress σ in is set here on the upper wall by assigning wall particles an appropriate velocity.

To impose specified stress on the top wall of the Couette geometry, the following approach is used. The wall particles are assigned with velocity ( v i w) computed from the wall force ( F diff), which is in turn determined from a predictor–corrector method,
v i w n + 1 = v i w n + K p F diff Δ t .
(37)
The net wall force ( F diff) is computed as the difference between the applied wall force ( F in) and the resistive shear force ( F r) exerted by the fluid on the top wall,
F diff = F in F r ,
(38)
F r = j w F j ,
(39)
F i w = j f F i j ν ,
(40)
where subscripts w and f denote wall and fluid particles, respectively. F i j ν is the interparticle viscous force in Eq. (9). All simulations use a quintic kernel, and smoothing length h sl = 1.1 Δ x, where Δ x = 0.02 l is the initial spatial resolution. Simulations are carried out for the same material parameters (see Table I) with different imposed stress values. Three microstructure diffusion coefficients are tested to probe the effects of the non-local term—strong non-locality ( α = 10 5 m 2 s 1), moderate non-locality ( α = 10 8 m 2 s 1), and zero non-locality ( α = 0 m 2 s 1). In addition, three values of volume fraction were simulated for the case of strong non-locality: CST ( ϕ = 0.48), moderate DST ( ϕ = 0.50), and strong DST ( ϕ = 0.54). Steady-state values are used to construct the flow curve (Fig. 4). All rheometric values are obtained by averaging over the entire domain and across time, with an equilibration time of 100 s. Convergence study was carried out in the context of purely local channel flow (Sec. IV B), where the obtained SPH velocity profiles were compared to the theoretical solution. Significant improvements in the solution accuracy were seen up to the resolution of 50 × 50 particles, beyond which only marginal improvements were observed. Hence, in both sections, resolution of 50 × 50 is employed. Both walls were constructed with 200 dummy particles each. Simulations were run for up to 72 h of computing time and 10 min of simulation time.
TABLE I.

Simulation parameters.

h (m) l (m) ρ (kg m−3) ηs (Pa s) ϕ ϕ m ϕ 0 σ ( Pa )
0.01  0.01  1000  0.001  0.48/0.50/0.54  0.562  0.693  0.005 
h (m) l (m) ρ (kg m−3) ηs (Pa s) ϕ ϕ m ϕ 0 σ ( Pa )
0.01  0.01  1000  0.001  0.48/0.50/0.54  0.562  0.693  0.005 

Far below onset stress, behavior is essentially Newtonian—all three measured properties (shear stress, shear rate, and microstructural parameter f) approach their respective steady-state values asymptotically and within relatively short timescales. Once the stress is increased, the time required to achieve a steady state is significantly increased, and over/under-shoots are present prior to achieving a steady state. This is likely due to the formation of microstructure being relatively slow compared to purely inertial transient effects—in the initial stages, the lack of structure allows shear rate to peak above the equilibrium value. This overshoot then accelerates formation of microstructure, increasing viscosity, which in turn adjusts the shear rate (or undershoots). To mitigate this behavior, we choose to initialize all simulations from a pre-defined homogeneous field f corresponding to the correct prediction of the WC at a given input stress, along with the appropriate shear rate via the upper wall velocity. For all values of α, simulation results [Figs. 4(a)–4(c)] on stable branches are in excellent agreement with the theoretical model predictions.

In the case of strong non-locality, the simulation results for all three volume fractions [Fig. 3(a)] match the WC model very well across the entire flow curve, including the region of negative gradient for the DST volume fractions. Both the stress and microstructure fields are entirely homogeneous [Figs. 4(d) and 4(g)] due to the strong non-local term, with minor deviations in the computed shear rate field, likely due to the imprecision associated with numerical approximations and low resolution.

FIG. 3.

Results for three volume fractions: 0.54 (circles), 0.50 (triangles), and 0.48 (diamonds). Strong non-locality ( α = 10 5 m2 s−1).

FIG. 3.

Results for three volume fractions: 0.54 (circles), 0.50 (triangles), and 0.48 (diamonds). Strong non-locality ( α = 10 5 m2 s−1).

Close modal
FIG. 4.

The results of the SPH simulation compared against the WC model (black) for strong [(a), (d), and (g)], moderate [(b), (e), and (h)], and zero [(c), (f), and (i)] non-locality. The top row (a)–(c) consists of domain-time averaged shear stress and shear rate signal (SPH in blue and red). A single simulation (red) in each diagram ( Σ E = 1.25) was chosen, and values of all SPH particles were plotted over the WC flow curve (d)–(f). Fields of the microstructure parameter f associated with the previously chosen simulations in (d)–(f) visualized in (g)–(i). For high non-locality (g), entire smooth fields of f and stress are achieved. In the case of moderate non-locality (h), two high stress regions are attached to the upper and lower walls with bands propagating in the flow direction. Zero non-locality (i) yields discontinuous microstructure fields.

FIG. 4.

The results of the SPH simulation compared against the WC model (black) for strong [(a), (d), and (g)], moderate [(b), (e), and (h)], and zero [(c), (f), and (i)] non-locality. The top row (a)–(c) consists of domain-time averaged shear stress and shear rate signal (SPH in blue and red). A single simulation (red) in each diagram ( Σ E = 1.25) was chosen, and values of all SPH particles were plotted over the WC flow curve (d)–(f). Fields of the microstructure parameter f associated with the previously chosen simulations in (d)–(f) visualized in (g)–(i). For high non-locality (g), entire smooth fields of f and stress are achieved. In the case of moderate non-locality (h), two high stress regions are attached to the upper and lower walls with bands propagating in the flow direction. Zero non-locality (i) yields discontinuous microstructure fields.

Close modal

In the purely local (α = 0) case, for all imposed stresses falling within the unstable region of negative flow curve gradient, seemingly spurious “jumps” in microstructure, local stress, and viscosity were observed, where individual particles attain either much higher or lower stress and microstructural parameter with respect to the imposed stress. This behavior had no impact on the numerical stability of the simulations, and a solution could be computed even in this extreme case. These jumps occurred with no clear spatial correlation, yet both the domain averaged properties and the properties measured at the wall remained in close agreement with the correct theoretically predicted values. By plotting all individual SPH particles for a given applied stress over the expected flow diagram [Fig. 4(f)], it is apparent that these jumps are directed toward the valid stable solutions of the WC model.

In principle, a DST fluid could undergo vorticity banding under stress-controlled regime, which, in the simplest case, would manifest as any stress applied within the unstable branch (the branch joining frictional and frictionless branches) being split vertically into two bands organized along the vorticity direction each with either higher or lower stress, but the same shear rate. The width of the bands is then decided by a lever rule such that the total average stress is equal to the applied stress. Such phenomena in DST materials have been experimentally observed in the work of Herle et al.,14 but Hermes et al.65 excluded the possibility of steady-state vorticity bands in non-Brownian suspensions on the account of particle migration and normal stresses. For a purely local model and a constant volume fraction field, no effects counteract such splitting, giving rise to a new steady state configuration, where all individual SPH particles occupy stable branches of the flow curve ( d f i d t = 0), while being split such that the stress control equation is also satisfied ( d v wall d t = 0).

The lack of spatial correlation in structure [Fig. 4(i)] is caused by the action of the microstructure evolution equation on the numerical noise introduced in the SPH solution. Suppose all particles are initialized exactly at the steady state solution, with small perturbations in a random direction in the rate-stress plane. The particles that find themselves to the right of the flow curve will experience d f i d t > 0, leading to an increase in the local viscosity. This will, in turn, increase the local shear stress, increasing the value of f i ̂, ensuring that the particle will continue to move upward until it reaches a stable upper branch where d f i d t = 0. An identical argument applies to the particles to the left of the curve, with the resultant downward movement instead. In the case of our simulations, the small perturbations are introduced by the inherent numerical noise associated with meshless calculation of the local shear rate γ ̇ i via Eq. (14). These small imprecisions do not have a spatial correlation; hence after stress-splitting amplification, we arrive at effectively random spatial distribution of microstructural parameter. The left-right split of particles should approach an even split ( ψ 1 = ψ 2 = 0.5); however, the flow curve is not symmetric about the horizontal imposed stress line, and such split would result in an increase in measured shear stress. This means that as the stress-splitting takes place, particles are simultaneously rearranged between the lower and upper branches via the stress-control mechanism of Eqs. (37)–(40) to achieve correct average stress.

To further explore both the microstructural transients and steady states, new dimensionless parameters are defined as follows:
ψ 1 = number on the lower branch total number ,
(41)
ψ 2 = number on the upper branch total number ,
(42)
ψ u = number on the unstable branch total number .
(43)
The number of particles on either the upper or lower branch is measured with reference to the two threshold stress values σ α and σ β, which occur at the intersection of the unstable branch with the upper and lower branches, respectively. Typical response consists of a 50 s period for banding process to begin and almost 150–200 s to reach a steady state. Different initial configurations were tested, including (a) relaxed microstructure, (b) correct steady-state microstructure, and (c) excessive level of structure. Due to the local nature of this process, the splitting is independent of the resolution, as there is no finite characteristic length scale associated with this process. All three configurations yielded the same steady-state banding ratios, i.e., the above-mentioned steady state is broadly independent of the path taken to achieve it. Furthermore, shear-controlled cases were tested both within and outside the multivariate region and yielded results with no apparent banding, a result consistent with the theory where only stress-controlled flow of DST fluids should experience vorticity banding.
To better understand the exact split between ψ 1 and ψ 2 at steady state, a simple model of DST vorticity banding is considered. Suppose that stress is applied within the unstable region of the flow curve at a coordinate ( γ ̇ u , σ u ). This unstable point will split directly in the vertical direction until it meets the lower and upper branches, resulting in a new steady-state consisting of two stable points ( γ ̇ u , σ 1 ) and ( γ ̇ u , σ 2 ), respectively. We also assume that the total average stress is preserved such that
σ u = ψ 1 σ 1 + ψ 2 σ 2 .
(44)
By considering that the system is now in a steady state ( ψ 1 + ψ 2 = 1 and ψ u = 0), both fractions can be uniquely determined for any input stress. It is immediately apparent that the exact split will be governed by the exact shape of the flow curve. Finally, we define a dimensionless stress coordinate q using critical stress values σ β and σ α [see Fig. 4(f)],
q = σ u σ β σ α σ β ,
(45)
such that, for q < 0, input stress is on the lower stable branch, 0 q 1 input stress is on the unstable branch, and for q > 1, input stress is on the upper stable branch. The results along with simulations are presented in Fig. 5. The simulation results exhibit a significantly lower number of particles on the lower branch and a higher number of particles on the upper branch relative to the theoretical predictions. This is due to how the split rearrangement is achieved. Since an even split would result in an increase in the measured shear stress, initial splitting causes the stress-control scheme to reduce shear rate. The shear rate reduction is set by the decrease in wall velocity, which affects all particles, causing a shift to the left in the ( γ ̇σ) phase diagram. Particles on the upper branch can only shift to the lower branch by crossing the critical stress σ β, where particles can no longer follow the linear relationship between shear rate and shear stress and “fall” down to the stable branch, due to the combined effects of the microstructure evolution equation and the stress-control scheme. This results in the particles on the upper stable branch, occupying shear stress and shear rate values lower than anticipated in the naive vertical splitting scenario. A lower stress value of σ 1 in Eq. (44) necessitates a higher value of ψ 2 (and a lower value of ψ 1) at steady state. This mechanism also accounts for the poor agreement between the simulations and the WC model in Fig. 4(c)—average shear stress is close to the correct input value, but the splitting causes much lower shear rate values.
FIG. 5.

Theoretically predicted branch split ratios (lines) vs simulations. Simulation values represented by black circles ( ψ 1) and triangles ( ψ 2).

FIG. 5.

Theoretically predicted branch split ratios (lines) vs simulations. Simulation values represented by black circles ( ψ 1) and triangles ( ψ 2).

Close modal

Finally, the case of moderate non-locality is considered. The averaged results show good agreement with the WC model Fig. 4(b) along the entire flow curve, just as in the homogeneous case. Unlike the homogeneous case, the plot of individual SPH particles in Fig. 4(e) shows a clear and significant spread in both shear stress and shear rate. However, this spread is different from the local case—particles occupy mainly the unstable region between the two stable branches. The bands in the microstructure parameter are also accompanied by bands in shear stress and shear rate. This steady state occurs as a result of a balance between the split-inducing local term and the smoothing non-local term in the microstructure evolution equation. Furthermore, the competition between the local and non-local effects yields an apparent spatial correlation in the microstructure field in Fig. 4(h) reminiscent of the bands observed in the work of Nakanishi et al.22 

1. Theoretical analysis

In this section, we consider flow in a planar channel (Fig. 6), with gap d = 0.01 m and periodic boundary conditions in the flow direction. Body force F parallel to the flow direction is applied to all particles. In the following, a derivation of the theoretical solution is provided. We begin by assuming a steady state, well-developed unidirectional [ u = ( u , 0 , 0 )] axial flow driven by a body force F. Equation (2) simplifies to
σ x y y = F .
(46)
The constitutive equation for the only non-zero component of the stress tensor is σ x y = η d u d y subject to no-slip and vanishing stress boundary conditions u = 0 at y = 0 and σ x y = 0 at y = h.
FIG. 6.

Sketch for the channel geometry. The flow is driven solely by a constant body force F.

FIG. 6.

Sketch for the channel geometry. The flow is driven solely by a constant body force F.

Close modal
Integrating Eq. (46) with appropriate boundary conditions yields
σ x y = F ( h y ) .
(47)
Using Eqs. (4) and (5), the constitutive relation can be rewritten to yield an ODE,
d u d y = F ( h y ) η s ( 1 ϕ ϕ m exp ( [ σ * F ( h y ) ] β ) + ϕ 0 ( 1 exp ( [ σ * F ( h y ) ] β ) ) 2 .
(48)

Equation (48) can be solved numerically to obtain shear rate, velocity, and viscosity profiles. In this work, the solution was obtained via ODE45 MATLAB algorithm with a relative tolerance of 10 8.

2. Numerical results

With the chosen parameters, we focus primarily on transitional DST regime where both frictional and frictionless states are meaningfully present within the flow and consider only the purely local case. At the steady state, the velocity shape of DST flow displays a sharper profile with an upward inflection near the surface [Fig. 7(c)] and is in good agreement with the theoretical solutions of our model. The inflection corresponds to the threshold transition region between primarily frictional and frictionless regimes.

FIG. 7.

Velocity profiles (left) and corresponding microstructural profiles (right) for dimensionless driving forces F h / σ * of (a) and (b), 10 (c) and (d), and 100 (e) and (f) for ϕ = 0.54.

FIG. 7.

Velocity profiles (left) and corresponding microstructural profiles (right) for dimensionless driving forces F h / σ * of (a) and (b), 10 (c) and (d), and 100 (e) and (f) for ϕ = 0.54.

Close modal

It is worth noting that for these values, the ratio between frictionless and frictional viscosity, η ¯, is approximately 30—far lower than typical magnitudes observed in experimental model systems such as cornstarch, where the viscosity can span multiple orders of magnitude during shear thickening.66,67 This low value of η ¯ was chosen to facilitate feasible simulation run time by avoiding excessively low or high viscosity at any point during shear thickening.

The resulting velocity profiles are presented in Fig. 7. The driving force (and thus the stress profile) was increased to probe the influence of shear thickening on the shape of the velocity profiles.

For sufficiently low driving forces, the velocity profile remains parabolic (a) due to the lack of developed microstructure (b)—viscosity remains close to the lower branch viscosity η0. The microstructural profile is in good agreement with the theoretical profile, with minor deviation near the centerline. This is due to the initialization of the structure from a finite value—particles near centerline experience near-zero shear rate, and thus, the structure cannot evolve to the appropriate steady state. This deviation is not significantly reflected in the velocity profile due to the weak scaling of viscosity with the microstructure at such low values, resulting in an excellent agreement in the velocity profile.

At intermediate driving forces, flow in the center ( y / D = 0.5) remains unconstrained; however, the microstructure builds up from the walls (d) toward the center leading to a sharpened velocity profile (c) with an upward inflection—feature characteristic of the S-shaped rheology. It is apparent that the stress-splitting is present here, with large deviations in the simulated microstructure from the theoretical values. Despite this discrepancy, the simulated velocity profile has very good agreement with the theory. The banding behavior seems to preserve average velocity, stress, and microstructural profiles within the channel.

Further increase in the driving force leads to almost completely shear frictional flow (f)—the velocity becomes broadly parabolic (e) again with the velocity corresponding to the upper branch viscosity η1; however, small inflection persists at the tip, where the transition takes place over a very thin region. Any banding behavior is constrained the small transitional slice, leading to a deteriorated microstructural solution near the centerline.

The model was compared to the simplest DST model, i.e., inverse biviscous.56 The biviscous model does not predict S-curves; d γ / d σ > 0 for all cases. The two models were compared by plotting the average velocity for a range of driving forces (Fig. 8) for an equivalent viscosity ratio. The biviscous model predicts two branches corresponding to lower and upper viscosities, joined together by a plateau where increasing the driving force does not affect the average flow. The model presented in this work exhibits similar low and high viscosity branches; however, the conjoining region exhibits a flow reduction pattern unique to the S-shape rheology. It is noteworthy that the presence of a negative flow curve gradient is not always sufficient on its own—a certain minimum viscosity ratio is required to achieve flow reduction behavior.

FIG. 8.

Comparison between the present model (blue) and biviscous model56 (red). Simulations (black) carried out for ϕ = 0.54 , ϕ 0 = 0.64, and ϕ m = 0.56 , u ¯ is the average velocity in the channel, and u ¯ * is an arbitrary average velocity chosen at F / F max = 1. The upper and lower viscosity in the biviscous model re-scaled to match frictionless and frictional viscosity.

FIG. 8.

Comparison between the present model (blue) and biviscous model56 (red). Simulations (black) carried out for ϕ = 0.54 , ϕ 0 = 0.64, and ϕ m = 0.56 , u ¯ is the average velocity in the channel, and u ¯ * is an arbitrary average velocity chosen at F / F max = 1. The upper and lower viscosity in the biviscous model re-scaled to match frictionless and frictional viscosity.

Close modal

In this work, we have presented the first implementation of a microstructural model of discontinuous shear thickening in SPH simulations. A simple non-local scalar microstructural evolution equation based on the work of Baumgarten and Kamrin23 was proposed and implemented in simple shear and channel geometries. Measured transients in shear rate were found to be in qualitative agreement with experimentally measured shear rates under constant imposed stress.

The simple shear results were used to construct flow curves and compare the macroscopic results with the microscopic WC model which we employ. Simulation results were in good agreement along lower and upper viscosity branches and predicted the characteristic S-shape transition between them. In the absence of non-locality, the stress-splitting instability has been observed as a result of the action of the local term in the microstructure evolution equation upon small numerical perturbations in the shear rate field. The resulting steady state was found to be characterized by correct domain-averaged shear stress, despite all SPH particles occupying stable upper and lower branches, lack of spatial correlation, and a low shear rate relative to the WC model. Both the low shear rate and the discrepancy between theoretical and measured split ratios were contextualized in terms of the particle rearrangement mechanism resulting from the stress-control scheme. Both moderate and strong non-locality yielded good agreement in average shear rate and shear stress, with banding spatial correlation being observed in the case of moderate non-locality.

Channel flow simulations showed excellent agreement in velocity profiles with the theoretical solution of the WC model across a range of driving forces. Once channel flow is forced into the intermediate DST regime, the velocity profile forms an upward inflection as a result of the S-shaped rheology. Stress-splitting was observed in the channel flow, with no discernable impact on the velocity profiles. Comparison in channel flow was made with the inverse biviscous mode.56 The present S-shape microstructure model produced a flow reduction behavior, in contrast to the flow plateau displayed by the biviscous model.

Further study is required to assess stability of the present model and implementation. Exhaustive assessment of the stress-splitting phenomena, band spatial correlation, and the associated characteristic length scale is required. Reproduction of known instability phenomena15,58 along with introduction of depth in the vorticity direction could provide further insight into the chaotic rheology of the DST materials.

We thank Zhongqiang Xiong for discussions. P.A. and B.S. acknowledge funding from the Engineering and Physical Sciences Research Council (EP/S034587/1). This research is partially supported by the Basque Government through the BERC 2022-2025 program and by the Spanish State Research Agency through BCAM Severo Ochoa excellence accreditation CEX2021-0011 42-S/MICIN/AEI/10.13039/501100011033 and through the project PID2020-117080RB-C55 (“Microscopic foundations of softmatter experiments: computational nano-hydrodynamics” and acronym “Compu-Nano-Hydro”). R.S. acknowledges funding from the National Natural Science Foundation of China (12174390 and 12150610463) and Wenzhou Institute, University of Chinese Academy of Sciences (WIUCASQD2020002).

The authors have no conflicts to disclose.

Peter Angerman: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead); Sagaya S. Prasanna Kumar: Conceptualization (equal); Resources (equal); Software (equal); Ryohei Seto: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Visualization (supporting); Writing – review & editing (equal). Bjornar Sandnes: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Marco Ellero: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (lead); Resources (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

Following the work of Lind et al.,68 we introduce shifting of particle position to disallow SPH particles from following the streamlines of the flow, which might lead to anisotropic particle distribution under certain conditions. It is pointed out that the accuracy of SPH approximations diminishes with irregular particle distribution. The shifting of particle position is completely numerical in origin and not physical. Therefore, we keep the shifting to a minimal level by setting the condition that maximum particle shifting x i , shift max < 0.1 h, where h is the smoothing length (approximately equal to the average interparticle spacing).

The particles are shifted based on a shifting velocity v i , s, which depends on the local particle concentration gradient vector C i,
v i , s = D i C i .
(A1)
Here, Di is the shifting velocity coefficient. The normalized particle concentration Ci at a given particle location is computed as
C i = j V j W i j ,
(A2)
where V is the volume of the particles. The gradient of the particle concentration is then calculated as
C i = j { 1 + 1 4 ( W i j W i i ) 4 } V j W i j .
(A3)
For more details on the treatment of free surface boundary conditions in SPH, the reader is referred to the work of Lind et al.68 
1.
J.
Mewis
and
N. J.
Wagner
,
Colloidal Suspension Rheology
, Cambridge Series in Chemical Engineering (
Cambridge University Press
,
2011
).
2.
R. I.
Tanner
, “
Review Article: Aspects of non-colloidal suspension rheology
,”
Phys. Fluids
30
(
10
),
101301
(
2018
).
3.
M. M.
Denn
and
J. F.
Morris
, “
Rheology of non-Brownian suspensions
,”
Annu. Rev. Chem. Biomol. Eng.
5
(
1
),
203
228
(
2014
).
4.
R. I.
Tanner
,
Engineering Rheology
, Oxford Engineering Science Series (
Oxford University Press
,
Oxford
,
2000
).
5.
J. F.
Morris
, “
Shear thickening of concentrated suspensions: Recent developments and relation to other phenomena
,”
Annu. Rev. Fluid Mech.
52
(
1
),
121
144
(
2020
).
6.
E.
Brown
and
H. M.
Jaeger
, “
Shear thickening in concentrated suspensions: Phenomenology, mechanisms and relations to jamming
,”
Rep. Prog. Phys.
77
(
4
),
046602
(
2014
).
7.
R.
Seto
,
R.
Mari
,
J. F.
Morris
, and
M. M.
Denn
, “
Discontinuous shear thickening of frictional hard-sphere suspensions
,”
Phys. Rev. Lett.
111
,
218301
(
2013
).
8.
M.
Wyart
and
M. E.
Cates
, “
Discontinuous shear thickening without inertia in dense non-Brownian suspensions
,”
Phys. Rev. Lett.
112
,
098302
(
2014
).
9.
J. E.
Thomas
,
A.
Goyal
,
D. S.
Bedi
,
A.
Singh
,
E.
Del Gado
, and
B.
Chakraborty
, “
Investigating the nature of discontinuous shear thickening: Beyond a mean-field description
,”
J. Rheol.
64
(
2
),
329
341
(
2020
).
10.
R.
Seto
,
A.
Singh
,
B.
Chakraborty
,
M.
Denn
, and
J.
Morris
, “
Shear jamming and fragility in dense suspensions
,”
Granul. Matter
21
,
82
(
2019
).
11.
M.
Gameiro
,
A.
Singh
,
L.
Kondic
,
K.
Mischaikow
, and
J. F.
Morris
, “
Interaction network analysis in shear thickening suspensions
,”
Phys. Rev. Fluids
5
,
034307
(
2020
).
12.
M.
Nabizadeh
,
A.
Singh
, and
S.
Jamali
, “
Structure and dynamics of force clusters and networks in shear thickening suspensions
,”
Phys. Rev. Lett.
129
,
068001
(
2022
).
13.
O.
Sedes
,
H. A.
Makse
,
B.
Chakraborty
, and
J. F.
Morris
, “
K-core analysis of shear-thickening suspensions
,”
Phys. Rev. Fluids
7
,
024304
(
2022
).
14.
V.
Herle
,
P.
Fischer
, and
E. J.
Windhab
, “
Stress driven shear bands and the effect of confinement on their structures—A rheological, flow visualization, and Rheo-Sals study
,”
Langmuir
21
(
20
),
9051
9057
(
2005
).
15.
J.
Richards
,
J.
Royer
,
B.
Liebchen
,
B.
Guy
, and
W.
Poon
, “
Competing timescales lead to oscillations in shear-thickening suspensions
,”
Phys. Rev. Lett.
123
,
038004
(
2019
).
16.
B.
Texier
,
H.
Lhuissier
,
Y.
Forterre
, and
B.
Metzger
, “
Surface-wave instability without inertia in shear-thickening suspensions
,”
Commun. Phys.
3
,
232
(
2020
).
17.
B.
Darbois Texier
,
H.
Lhuissier
,
B.
Metzger
, and
Y.
Forterre
, “
Shear-thickening suspensions down inclines: From Kapitza to Oobleck waves
,”
J. Fluid Mech.
959
,
A27
(
2023
).
18.
J. H.
Cho
,
A. H.
Griese
,
I. R.
Peters
, and
I.
Bischofberger
, “
Lasting effects of discontinuous shear thickening in cornstarch suspensions upon flow cessation
,”
Phys. Rev. Fluids
7
,
063302
(
2022
).
19.
V.
Rathee
,
D. L.
Blair
, and
J. S.
Urbach
, “
Dynamics and memory of boundary stresses in discontinuous shear thickening suspensions during oscillatory shear
,”
Soft Matter
17
,
1337
1345
(
2021
).
20.
V.
Rathee
,
D. L.
Blair
, and
J. S.
Urbach
, “
Localized transient jamming in discontinuous shear thickening
,”
J. Rheol.
64
(
2
),
299
308
(
2020
).
21.
V.
Rathee
,
J.
Miller
,
D. L.
Blair
, and
J. S.
Urbach
, “
Structure of propagating high-stress fronts in a shear-thickening suspension
,”
Proc. Natl. Acad. Sci. U. S. A.
119
(
32
),
e2203795119
(
2022
).
22.
H.
Nakanishi
,
S.
Nagahiro
, and
N.
Mitarai
, “
Fluid dynamics of dilatant fluids
,”
Phys. Rev. E
85
,
011401
(
2012
).
23.
A. S.
Baumgarten
and
K.
Kamrin
, “
A general constitutive model for dense, fine-particle suspensions validated in many geometries
,”
Proc. Natl. Acad. Sci. U. S. A.
116
(
42
),
20828
20836
(
2019
).
24.
J. J. J.
Gillissen
,
C.
Ness
,
J. D.
Peterson
,
H. J.
Wilson
, and
M. E.
Cates
, “
Constitutive model for time-dependent flows of shear-thickening suspensions
,”
Phys. Rev. Lett.
123
,
214504
(
2019
).
25.
J. J. J.
Gillissen
,
C.
Ness
,
J. D.
Peterson
,
H. J.
Wilson
, and
M. E.
Cates
, “
Constitutive model for shear-thickening suspensions: Predictions for steady shear with superposed transverse oscillations
,”
J. Rheol.
64
(
2
),
353
365
(
2020
).
26.
S. J.
Lind
,
B. D.
Rogers
, and
P. K.
Stansby
, “
Review of smoothed particle hydrodynamics: Towards converged Lagrangian flow modelling
,”
Proc. R. Soc. London A
476
,
20190801
20192020
(
2241
).
27.
M.
Ellero
,
M.
Serrano
, and
P.
Español
, “
Incompressible smoothed particle hydrodynamics
,”
J. Comput. Phys.
226
(
2
),
1731
1752
(
2007
).
28.
R. A.
Gingold
and
J. J.
Monaghan
, “
Smoothed particle hydrodynamics: Theory and application to non-spherical stars
,”
Mon. Not. R. Astron. Soc.
181
(
3
),
375
389
(
1977
).
29.
L. B.
Lucy
, “
A numerical approach to the testing of the fission hypothesis
,”
Astron. J.
82
,
1013
1024
(
1977
).
30.
J. J.
Monaghan
, “
Smoothed particle hydrodynamics
,”
Rep. Prog. Phys.
68
(
8
),
1703
(
2005
).
31.
A. M.
Abdelrazek
,
I.
Kimura
, and
Y.
Shimizu
, “
Simulation of three-dimensional rapid free-surface granular flow past different types of obstructions using the SPH method
,”
J. Glaciol.
62
(
232
),
335
347
(
2016
).
32.
G.
Reece
,
B. D.
Rogers
,
S.
Lind
, and
G.
Fourtakas
, “
New instability and mixing simulations using SPH and a novel mixing measure
,”
J. Hydrodyn.
32
(
4
),
684
698
(
2020
).
33.
F.
Kalateh
and
A.
Koosheh
, “
Simulation of cavitating fluid-structure interaction using SPH-FE method
,”
Math. Comput. Simul.
173
,
51
70
(
2020
).
34.
N.
Moreno
and
M.
Ellero
, “
Generalized Lagrangian heterogeneous multiscale modelling of complex fluids
,”
J. Fluid Mech.
969
,
A2
(
2023
).
35.
N.
Moreno
,
S.
Nunes
, and
V.
Calo
, “
Morphological transitions of block copolymer micelles: Implications for isoporous membrane ordering
,” arXiv:2304.03900 (
2023
).
36.
L.
Santelli
,
A.
Vázquez-Quesada
, and
M.
Ellero
, “
SPH simulations of integral multimode and fractional viscoelastic models
” (unpublished).
37.
M.
Ellero
,
M.
Kröger
, and
S.
Hess
, “
Viscoelastic flows studied by smoothed particle dynamics
,”
J. Non-Newtonian Fluid Mech.
105
(
1
),
35
51
(
2002
).
38.
M.
Ellero
and
R. I.
Tanner
, “
SPH simulations of transient viscoelastic flows at low Reynolds number
,”
J. Non-Newtonian Fluid Mech.
132
(
1
),
61
72
(
2005
).
39.
J.
Fang
,
R. G.
Owens
,
L.
Tacher
, and
A.
Parriaux
, “
A numerical study of the SPH method for simulating transient viscoelastic free surface flows
,”
J. Non-Newtonian Fluid Mech.
139
,
68
84
(
2006
).
40.
J.
Fang
,
A.
Parriaux
,
M.
Rentschler
, and
C.
Ancey
, “
Improved SPH methods for simulating free surface flows of viscous fluids
,”
Appl. Numer. Math.
59
(
2
),
251
271
(
2009
).
41.
A.
Rafiee
,
M. T.
Manzari
, and
M.
Hosseini
, “
An incompressible SPH method for simulation of unsteady viscoelastic free-surface flows
,”
Int. J. Non Linear Mech.
42
(
10
),
1210
1223
(
2007
).
42.
X.
Xu
,
J.
Ouyang
,
B.-X.
Yang
, and
Z.
Liu
, “
SPH simulations of three-dimensional non-Newtonian free surface flows
,”
Comput. Methods Appl. Mech. Eng.
256
,
101
116
(
2013
).
43.
A.
Vázquez-Quesada
and
M.
Ellero
, “
SPH simulations of a viscoelastic flow around a periodic array of cylinders confined in a channel
,”
J. Non-Newtonian Fluid Mech.
167–168
,
1
8
(
2012
).
44.
M.
Grilli
,
A.
Vázquez-Quesada
, and
M.
Ellero
, “
Transition to turbulence and mixing in a viscoelastic fluid flowing inside a channel with a periodic array of cylindrical obstacles
,”
Phys. Rev. Lett.
110
,
174501
(
2013
).
45.
X.
Xu
and
X.-L.
Deng
, “
An improved weakly compressible SPH method for simulating free surface flows of viscous and viscoelastic fluids
,”
Comput. Phys. Commun.
201
,
43
62
(
2016
).
46.
M. X.
Rodriguez-Paz
and
J.
Bonet
, “
A corrected smooth particle hydrodynamics method for the simulation of debris flows
,”
Numer. Methods Partial Differ. Equations
20
(
1
),
140
163
(
2004
).
47.
L.
Minatti
and
E.
Paris
, “
A SPH model for the simulation of free surface granular flows in a dense regime
,”
Appl. Math. Model.
39
(
1
),
363
382
(
2015
).
48.
O.
Pouliquen
,
C.
Cassar
,
P.
Jop
,
Y.
Forterre
, and
M.
Nicolas
, “
Flow of dense granular material: Towards simple constitutive laws
,”
J. Stat. Mech.
2006
(
7
),
P07020
.
49.
P.
Jop
,
Y.
Forterre
, and
O.
Pouliquen
, “
A constitutive law for dense granular flows
,”
Nature
441
(
7094
),
727
730
(
2006
).
50.
G.
Fourtakas
and
B. D.
Rogers
, “
Modelling multi-phase liquid-sediment scour and resuspension induced by rapid flows using smoothed particle hydrodynamics (SPH) accelerated with a graphics processing unit (GPU)
,”
Adv. Water Resour.
92
,
186
199
(
2016
).
51.
M.
Khanpour
,
A. R.
Zarrati
,
M.
Kolahdoozan
,
A.
Shakibaeinia
, and
S. M.
Amirshahi
, “
Mesh-free SPH modeling of sediment scouring and flushing
,”
Comput. Fluids
129
,
67
78
(
2016
).
52.
H.
Zhu
,
N. S.
Martys
,
C.
Ferraris
, and
D.
De Kee
, “
A numerical study of the flow of Bingham-like fluids in two-dimensional vane and cylinder rheometers using a smoothed particle hydrodynamics (SPH) based method
,”
J. Non-Newtonian Fluid Mech.
165
(
7
),
362
375
(
2010
).
53.
A.
Paiva
,
F.
Petronetto
,
T.
Lewiner
, and
G.
Tavares
, “
Particle-based viscoplastic fluid/solid simulation
,”
Comput.-Aided Des.
41
(
4
),
306
314
(
2009
).
54.
C.
Zhu
,
Y.
Huang
, and
L.
Zhan
, “
SPH-based simulation of flow process of a landslide at Hongao landfill in China
,”
Nat. Hazards
93
(
3
),
1113
1126
(
2018
).
55.
E.
Rossi
,
I.
Garcia de Beristain
,
A.
Vazquez-Quesada
,
J. E.
López-Aguilar
, and
M.
Ellero
, “
SPH simulations of thixo-viscoplastic fluid flow past a cylinder
,”
J. Non-Newtonian Fluid Mech.
308
,
104891
(
2022
).
56.
A.
Vázquez-Quesada
,
N. J.
Wagner
, and
M.
Ellero
, “
Planar channel flow of a discontinuous shear-thickening model fluid: Theory and simulation
,”
Phys. Fluids
29
(
10
),
103104
(
2017
).
57.
S. H.
Maron
and
P. E.
Pierce
, “
Application of ree-eyring generalized flow theory to suspensions of spherical particles
,”
J. Colloid Sci.
11
(
1
),
80
95
(
1956
).
58.
R.
Chacko
,
R.
Mari
,
M.
Cates
, and
S.
Fielding
, “
Dynamic vorticity banding in discontinuously shear thickening suspensions
,”
Phys. Rev. Lett.
121
,
108003
(
2018
).
59.
G.
Bossis
,
P.
Boustingorry
,
Y.
Grasselli
,
A.
Meunier
,
R.
Morini
,
A.
Zubarev
, and
O.
Volkova
, “
Discontinuous shear thickening in the presence of polymers adsorbed on the surface of calcium carbonate particles
,”
Rheol. Acta
56
(
5
),
415
430
(
2017
).
60.
K.
Kamrin
and
D. L.
Henann
, “
Nonlocal modeling of granular flows down inclines
,”
Soft Matter
11
,
179
185
(
2015
).
61.
J. P.
Morris
,
P. J.
Fox
, and
Y.
Zhu
, “
Modeling low Reynolds number incompressible flows using SPH
,”
J. Comput. Phys.
136
(
1
),
214
226
(
1997
).
62.
S.
Adami
,
X. Y.
Hu
, and
N. A.
Adams
, “
A generalized wall boundary condition for smoothed particle hydrodynamics
,”
J. Comput. Phys.
231
(
21
),
7057
7075
(
2012
).
63.
R.
Courant
,
K.
Friedrichs
, and
H.
Lewy
, “
Über die partiellen differenzengleichungen der mathematischen physik
,”
Math. Ann.
100
(
1
),
32
74
(
1928
).
64.
A.
Vázquez-Quesada
,
P.
Español
,
R. I.
Tanner
, and
M.
Ellero
, “
Shear thickening of a non-colloidal suspension with a viscoelastic matrix
,”
J. Fluid Mech.
880
,
1070
1094
(
2019
).
65.
M.
Hermes
,
B. M.
Guy
,
W. C. K.
Poon
,
G.
Poy
,
M. E.
Cates
, and
M.
Wyart
, “
Unsteady flow and particle migration in dense, non-Brownian suspensions
,”
J. Rheol.
60
(
5
),
905
916
(
2016
).
66.
I. R.
Peters
,
S.
Majumdar
, and
H. M.
Jaeger
, “
Direct observation of dynamic shear jamming in dense suspensions
,”
Nature
532
(
7598
),
214
217
(
2016
).
67.
D.
Ozturk
,
M. L.
Morgan
, and
B.
Sandnes
, “
Flow-to-fracture transition and pattern formation in a discontinuous shear thickening fluid
,”
Commun. Phys.
3
,
119
(
2020
).
68.
S. J.
Lind
,
R.
Xu
,
P. K.
Stansby
, and
B. D.
Rogers
, “
Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves
,”
J. Comput. Phys.
231
(
4
),
1499
1523
(
2012
).
Published open access through an agreement withJISC Collections