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.

## I. INTRODUCTION

Suspensions of both colloidal^{1} and non-colloidal^{2,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 rheometric^{14,15} and complex flows. For example, Texier *et al.*^{16,17} demonstrated that the negative gradients of the flow curve ( $ d \gamma \u0307 d \sigma < 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 Kamrin^{23} 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 (*N*_{1}) and too low (*N*_{2}) 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^{28} and Lucy^{29} 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 naturally^{30}—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 Ellero^{43} 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 Bonet^{46} simulated debris flow down an incline using the Bingham model and the generalized viscoplastic model. Minatti and Paris^{47} 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.

## II. GOVERNING EQUATIONS

*p*is the pressure, and

**is a body force acting on fluid. In this DST model, the local suspension viscosity is determined using the Maron–Pierce expression,**

*F*^{57}

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 \u2009 ( \u2212 \sigma \u2217 / \sigma )$. The repulsive force between particles, *F*_{0}, sets the characteristic stress scale at the onset of shear thickening $ \sigma \u2217 = F 0 / 6 \pi 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 ( $ \varphi < 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 ( $ \varphi = \varphi m$).

^{17,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

*K*. In addition to the usual local term, we introduce a non-local contribution in the dynamics of structure, where the microstructural parameter

_{f}*f*diffuses with a diffusion coefficient

*α*. In this work, we follow the approach of Kamrin and Henann

^{60}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,

*α*is the diffusion coefficient.

## III. SPH–DST MODEL

### A. SPH equations of motions

^{61}

*m*is the mass of a particle, $ W i j = W ( | x i j | )$ is the smoothing kernel, and $ v i j = v i \u2212 v j , \u2009 x i j = x i \u2212 x j$, and $ r i j = x i j \xb7 x i j$.

_{j}*ρ*

_{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$.

### B. SPH–DST microstructure model

*i*indicates the particle index (that is $ f ( x i , t ) = f i$),

*f*is calculated to evaluate $ \varphi J$ in Eq. (12), which in turn is used to obtain viscosity via Eq. (13). The value of $ \gamma \u0307$ in Eq. (11) is computed as the contraction of the strain rate tensor $ \gamma \u0307$ as follows:

_{i}*α*-component of the velocity vector in

*β*-direction is obtained as

^{26}

*f*, the choice of $\varphi $ has to be sufficiently close to $ \varphi m$ which in the above-mentioned case is $ \varphi > 0.50$. However, the value of $ \eta / \eta s$ for $ \varphi = 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:

*β*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 $ \varphi m$ and $ \varphi 0$, such that all simulations lie within the prescribed region of $ 10 \u2212 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

*f*and $ \varphi J ( f i )$ are obtained, the dynamic viscosity of a SPH particle is updated using Eq. (13).

_{i}### C. Solid wall boundary modeling

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.

*w*and

*f*denote the dummy wall and fluid particles, respectively, $ v i \u2208 w$ is the prescribed wall velocity, and

### D. Time integrators

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

- Predictor step:$ v i n + 1 2 = v i n + \Delta t 2 ( d v i d t ) n ,$$ \rho i n + 1 2 = \rho i n + \Delta t 2 ( d \rho i d t ) n ,$$ x i n + 1 2 = x i n + \Delta t 2 ( v i ) n ,$$ f i n + 1 2 = f i n + \Delta t 2 ( d f i d t ) n .$
- Corrector step:$ v i n + 1 2 = v i n + \Delta t 2 ( d v i d t ) n + 1 2 ,$$ \rho i n + 1 2 = \rho i n + \Delta t 2 ( d \rho i d t ) n + 1 2 ,$$ x i n + 1 2 = x i n + \Delta t 2 ( v i ) n + 1 2 ,$$ f i n + 1 2 = f i n + \Delta t 2 ( d f i d t ) n + 1 2 .$

*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 $ \Delta t$ is determined based on the following condition:

*ν*is the kinematic viscosity. The adaptive time step approach ensures that the CFL condition

^{63}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).

## IV. NUMERICAL SIMULATIONS

### A. Stress-imposed shear flow

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 $ \sigma in$ is set here on the upper wall by assigning wall particles an appropriate velocity.

*w*and

*f*denote wall and fluid particles, respectively. $ F i j \nu $ is the interparticle viscous force in Eq. (9). All simulations use a quintic kernel, and smoothing length $ h sl = 1.1 \Delta x$, where $ \Delta 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 ( $ \alpha = 10 \u2212 5 m 2 s \u2212 1$), moderate non-locality ( $ \alpha = 10 \u2212 8 m 2 s \u2212 1$), and zero non-locality ( $ \alpha = 0 m 2 s \u2212 1$). In addition, three values of volume fraction were simulated for the case of strong non-locality: CST ( $ \varphi = 0.48$), moderate DST ( $ \varphi = 0.50$), and strong DST ( $ \varphi = 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.

h (m)
. | l (m)
. | ρ (kg m^{−3})
. | η (Pa s)
. _{s} | $\varphi $ . | $ \varphi m$ . | $ \varphi 0$ . | $ \sigma \u2217 \u2009 ( 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})
. | η (Pa s)
. _{s} | $\varphi $ . | $ \varphi m$ . | $ \varphi 0$ . | $ \sigma \u2217 \u2009 ( 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.

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 \u0302$, 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 $ \gamma \u0307 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 ( $ \psi 1 = \psi 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.

*q*using critical stress values $ \sigma \beta $ and $ \sigma \alpha $ [see Fig. 4(f)],

*q*< 0, input stress is on the lower stable branch, $ 0 \u2264 q \u2264 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 ( $ \gamma \u0307$–

*σ*) phase diagram. Particles on the upper branch can only shift to the lower branch by crossing the critical stress $ \sigma \beta $, 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 $ \sigma 1$ in Eq. (44) necessitates a higher value of $ \psi 2$ (and a lower value of $ \psi 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.

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}

### B. Channel flow

#### 1. Theoretical analysis

*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

*u*= 0 at

*y*= 0 and $ \sigma x y = 0$ at

*y*=

*h*.

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 \u2212 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.

It is worth noting that for these values, the ratio between frictionless and frictional viscosity, $ \eta \xaf$, 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 $ \eta \xaf$ 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 \gamma / d \sigma > 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.

## V. CONCLUSIONS

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^{23} 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 phenomena^{15,58} along with introduction of depth in the vorticity direction could provide further insight into the chaotic rheology of the DST materials.

## ACKNOWLEDGMENTS

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX: PARTICLE SHIFTING TECHNIQUE

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

*D*is the shifting velocity coefficient. The normalized particle concentration

_{i}*C*at a given particle location is computed as

_{i}*V*is the volume of the particles. The gradient of the particle concentration is then calculated as

*et al.*

^{68}

## REFERENCES

*Colloidal Suspension Rheology*

*Engineering Rheology*

*K*-core analysis of shear-thickening suspensions