Microstructural Smoothed Particle Hydrodynamics Model and Simulations of Discontinuous Shear-Thickening Fluids

Despite the recent interest in the discontinuous shear-thickening (DST) behaviour, 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: weak, 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. Weak non-locality produced stress-splitting instability, resulting in discontinuous stress fields and poor agreement with the WC model. The mechanism of the stress-splitting has been explored and contextualised by the interaction of local microstructure evolution and the stress-control scheme. Velocity profiles obtained in body force-driven channel flow were found to be in excellent agreement with the analytical solution, yielding an upward inflection corresponding to the typical S-curve. 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 noncolloidal [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) behaviour [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 behaviour 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 a) Electronic mail: 2037484@swansea.ac.uk 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.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 behaviours in rheometric [9,10] and complex flows.For example, Texier et al. [11,12] demonstrated that the negative gradients of flow curve ( d γ dσ < 0) present in DST materials cause instability in free surface flows of cornstarch suspension, leading to roll-waves.
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 macrohydrodynamic flow properties through a transient evolution equation.Nakanishi et al. [13] 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 [14] proposed a scalar model based on the balance of strain-dependant hardening and softening attributed to load buckling, shearing, and electrostatic repulsion.Their model predicted rheological behaviour in good agreement with the experimental evidence.
Another development comes from the work of Gillissen et al. [15], 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 behaviour in more complicated flow arrangements, including flow reversal and superposed transverse oscillations [16], 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 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) [17,18].SPH was originally developed by Gingold and Monaghan [19] and Lucy [20] 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 [17].On the other hand, SPH is inherently conservative and handles free surface boundaries and interfacial multiphase flows naturally [21] -which is of importance in a number of studies, including the behaviour of granular bed loads [22], mixing [23], and cavitating flows [24].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) [25], with an example application in biologically relevant flows [26], or by directly solving integral constitutive equations [27].
Beyond Newtonian fluids, SPH has been used in a host of fields, including solid and complex fluid mechanics.In 2002, Ellero et al. [28] simulated viscoelastic flow by incorporating the corotational Maxwell model in SPH, followed by Oldroyd-B [29].Since then, Fang et al. [30,31] studied the behaviour of free surface flows of an Oldroyd-B fluid, with a focus on addressing tensile instability.Rafiee et al. [32] implemented both Maxwell and Oldroyd-B models in SPH free surface flows, including impacting drop and jet buckling.Xu et al. [33] simulated 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 [34] and Grilli et al. [35], and in extrudate swell by Xu and Deng [36].
Furthermore, a range of viscoplastic and elastoviscoplastic models have been implemented in SPH.Rodriguez-Paz and Bonet [37] simulated debris flow down an incline using the Bingham model and the generalized viscoplastic model.Minatti and Paris [38] considered free surface granular flows, with a constitutive relation based on the works of Pouliquen et al. [39] and Jop et al. [40], and validated their model against granular column collapse experiments.Furthermore, Bingham-like viscoplastic models have been applied to sedimenting flows [41,42], rheometric flows [43] and mixed fluid-structure cases [44,45].Recently, Rossi et al. [46] implemented a modified microstructural Papanastasiou model to study thixoviscoplastic flows around a cylinder.The first application of SPH to DST flow was carried out in the work of Vázquez-Quesada et al. [47], 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 microstructural DST rheological model in an SPH scheme.Unlike previous models [47], we obtain typical S-shaped DST rheology in the simple shear and characteristically inflexed velocity profiles accompanied by flow reduction behaviour in channel flow.In both cases, stress-splitting is observed.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.Finally, results for simple shear and planar channel flow are presented in sections IV A and IV B, respectively.

II. GOVERNING EQUATIONS
In this section, we outline the main equations used to model our DST fluid system.The continuum sys-tem can be written in the Lagrangian frame: where d/dt ≡ ∂/∂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 [48].
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 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 Wyarts-Cates model [8], in which the value of f is determined as exp(−σ * /σ).Here, σ * is the onset stress at which the suspension begins to shearthicken.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 shearthickening (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 volume fraction sharpen the S-curve up to a point of divergence (ϕ = ϕ m ).
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.Baumgarten and Kamrin [14] proposed an evolution equation for f considering the influence of different physical mechanisms towards the formation and destruction of the frictional contacts between suspended particles.In the present study, we propose a modified evolution equation combining some features of WC and Baumgarten's model.The evolution equation for f is set as where K f is the microstructural rate constant governing formation and destruction of structure, and γ(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 ( 6) accounts for the local evolution of microstructure, whereas the second term introduces non-local effects via microstructure diffusivity, where α is the diffusion coefficient.

III. SPH-DST MODEL A. SPH equations of motions
The Navier-Stokes equations are solved using the SPH methodology.A possible discretized form of the governing equations is given below [49]: where m j is the mass of a particle, 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: 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 .

B. SPH-DST microstructure model
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 In (11), f i is calculated to evaluate ϕ J in (12), which in turn is used to obtain viscosity via (13).The value of γ in (11) is computed as the contraction of the strain rate tensor γ as follows: The SPH approximation of the gradient of the αcomponent of the velocity vector in β-direction is obtained as The non-local term in equation ( 11) is discretised in a manner consistent with the Morris formulation of the moment equation [17], To simulate the 'S' -curve (see Fig. 1) with the above definition of order parameter f , the choice of ϕ has to be sufficiently close to ϕ m which in the above 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.
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 behaviour without increasing viscosity to computationally prohibitive values.The rest of the parameters are obtained by fitting the model to experimental data and other models.Once the value of f i and ϕ J (f i ) are obtained, the dynamic viscosity of a SPH particle is updated using (13).

C. Solid wall boundary modeling
The dummy particle method of Adami et al. [50] is employed for imposing no-slip velocity boundary conditions and impermeability conditions.Pressure, density, and velocity are assigned to the dummy par- ticles as where subscript w and f denote the dummy wall and fluid particles, respectively, v i∈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.
(1) Predictor step: ρ x (2) Corrector step: x The flow variables at the subsequent time step are obtained as follows: where n represents the current time instant, n + , (36) where h sl is the smoothing length, and ν is the kinematic viscosity.

IV. NUMERICAL SIMULATIONS A. Stress-Imposed Shear Flow
The behaviour 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 [51], 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.
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 where subscripts w and f denote wall and fluid particles, respectively.F ν ij is the inter-particle viscous force in eq. ( 9).All simulations use a quintic kernel, and smoothing length h sl = 1.1∆x,where ∆x = 0.02l 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 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 weak nonlocality (α = 0 m 2 s −1 ).In addition, three volume Far below onset stress, behaviour 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 behaviour, we choose to initialise 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 (Fig. 4(a)-(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 Fig. 4(d) and (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 behaviour 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 towards 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 organised 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. [9], but Hermes et al. [52] excludes 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 ( dfi dt = 0), while being split such that the stress control equation is also satisfied ( dv wall dt = 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 initialised exactly at the steady state solution, with small perturbation in a random direction in the rate-stress plane.The particles that find themselves to the right of the flow curve will experience dfi dt > 0, leading to an increase in local viscosity.This will, in turn, increase the local shear stress, increasing the value of fi , ensuring that the particle will continue to move upwards until it reaches a stable upper branch where dfi dt = 0. 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 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 leftright 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 of measured shear stress.This means that as the stress-splitting takes place, particles are simultaneously rearranged between the lower and upper branches via the stresscontrol 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: number on the lower branch total number , ψ 2 = number on the upper branch total number , 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 s to 200 s to reach a steady state.Different initial configurations were tested including a) relaxed microstructure, b) correct steadystate 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, shearcontrolled 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: By considering that the system is now in a steady FIG.5: Theoretically predicted branch split ratios (lines) vs. simulations.Simulation values represented by black circles (ψ 1 ) and triangles (ψ 2 ) .
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)): 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.Lower stress value of σ 1 in eq. ( 44) necessitates higher value of ψ 2 (and lower value of ψ 1 ) at steady state.This mechanism also accounts for the poor agreement between the simulations and the WC model 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 Fig. 4(h) reminiscent of the bands observed in the work of Nakanishi et al. [13].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 .Eq. ( 2) simplifies to The constitutive equation for the only non-zero component of the stress tensor is σ xy = η du dy subject to no-slip and vanishing stress boundary conditions u = 0 at y = 0 and σ xy = 0 at y = h.Integrating eq. ( 46) with appropriate boundary conditions yields Using eqs.( 4) and ( 5), the constitutive relation can be rewritten to yield an ODE: Eq. ( 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 .

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 are 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 [53,54].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 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 centre line.This is due to the initialisation of the structure from a finite value -particles near centre line experience nearzero shear rate, and thus the structure can not 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 centre (y/D = 0.5) remains unconstrained; however, the microstructure builds up from the walls (d) towards the centre leading to a sharpened velocity profile (c) with an upward inflection.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 behaviour 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 ve-locity 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 behaviour is constrained the small transitional slice, leading to a deteriorated microstructural solution near the centre line.The model was compared to the simplest DST model, i.e., inverse biviscous [47].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 behaviour.

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 [14] 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 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 characterised 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 low shear rate and the discrepancy between theoretical and measured split ratios were contextualised in terms of the particle rearrangment 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.Stresssplitting 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 [47].The present S-shape microstructure model produced a flow reduction behaviour, 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 [10,55] along with introduction of depth in the vorticity direction could provide further insight into the chaotic rheology of the DST materials.
The data that support the findings of this study are available from the corresponding author upon reasonable request.

FIG. 4 :
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 weak ((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 visualised in (g)-(i).

TABLE I :
Simulation Parameters