Microstructural smoothed particle hydrodynamics model and simulations of discontinuous shear-thickening fluids Open Access

Suspensions of both colloidal¹ and non-colloidal^2,3 particles are ubiquitous in industrially relevant flows,⁴ 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,⁵ 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,⁶ with the most widely adopted being the mechanism of Seto et al.⁷ or Wyart–Cates (WC) model,⁸ 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.,¹³ 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 rheometric^14,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.¹⁸ 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.¹⁹ employed boundary stress microscopy (BSM) to study stress distribution in rheometric flows. Their results revealed a high degree of heterogeneity in stress fields,²⁰ including two separate high-stress regions in contact with upper and lower plates, respectively, and propagation of high-stress fronts.²¹ 

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.²² 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 Kamrin²³ 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.,²⁴ 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,²⁵ although the values of normal stress differences predicted by their model were quantitatively too high (N₁) and too low (N₂) 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 Monaghan²⁸ and Lucy²⁹ 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.²⁶ On the other hand, SPH is inherently conservative and handles free surface boundaries and interfacial multiphase flows naturally³⁰—which is of importance in a number of studies, including the behavior of granular bed loads,³¹ mixing,³² and cavitating flows.³³ 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),³⁴ with an example application in biologically relevant flows,³⁵ or by directly solving integral constitutive equations.³⁶ 

Beyond Newtonian fluids, SPH has been used in a host of fields, including solid and complex fluid mechanics. In 2002, Ellero et al.³⁷ simulated viscoelastic flow by incorporating the corotational Maxwell model in SPH, followed by Oldroyd-B.³⁸ 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.⁴¹ implemented both Maxwell and Oldroyd-B models in SPH free-surface flows, including impacting drop and jet buckling. Xu et al.⁴² 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 Ellero⁴³ and Grilli et al.,⁴⁴ and in extrudate swell by Xu and Deng.⁴⁵ 

Furthermore, a range of viscoplastic and elasto-viscoplastic models have been implemented in SPH. Rodriguez-Paz and Bonet⁴⁶ simulated debris flow down an incline using the Bingham model and the generalized viscoplastic model. Minatti and Paris⁴⁷ considered free surface granular flows, with a constitutive relation based on the works of Pouliquen et al.⁴⁸ and Jop et al.,⁴⁹ and validated their model against granular column collapse experiments. Furthermore, Bingham-like viscoplastic models have been applied to sedimenting flows,^50,51 rheometric flows,⁵² and mixed fluid–structure cases.^53,54 Recently, Rossi et al.⁵⁵ 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.,⁵⁶ 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,⁵⁶ 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.

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,⁸ in which the value of f is determined as $exp ( − σ ∗ / σ )$⁠. The repulsive force between particles, F₀, sets the characteristic stress scale at the onset of shear thickening $σ ∗ = F 0 / 6 π a 2$⁠.¹¹  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$⁠).

The dummy particle method of Adami et al.⁶² 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.

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 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,⁶⁴ input stress $σ in$ is set here on the upper wall by assigning wall particles an appropriate velocity.

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.

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.,¹⁴ but Hermes et al.⁶⁵ 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.

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.²² 

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$⁠.

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.

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 η₀. 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 η₁; 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.⁵⁶ 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.

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 Kamrin²³ 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.⁵⁶ 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 phenomena^15,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.,⁶⁸ 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).