Small- and large-amplitude oscillatory rheometry with bead–spring dumbbells in Stokesian Dynamics to mimic viscoelasticity

In many areas of suspension mechanics, such as filled polymer fluids or household products such as toothpaste, the suspending fluid itself is inherently non-Newtonian and may exhibit viscoelastic properties. In this paper, we extend the Stokesian Dynamics formalism to incorporate a simple model of viscoelasticity by using small spheres as ‘beads’ in a bead–spring dumbbell (such as is found in the derivation of Oldroyd and FENE constitutive models for dilute polymer solutions). Various different spring laws are then tested in both small-amplitude and largeamplitude oscillatory shear, and their rheological behaviour is compared to continuum constitutive models.


Introduction
Suspensions of particles in fluids can be found both in nature and as the basis of many products in industry. Blood, ceramics, paper pulp, paint and adhesives, to name just a few, can all be characterised as a background fluid in which small particles are distributed. A popular simulation technique for these suspensions is Stokesian Dynamics [1] : a microhydrodynamic, low Reynolds number approach to modelling suspensions which considers the interaction of particles with each other against a Newtonian background solvent. Typically it is chosen for its suitability for three-dimensional simulation with low calculation and time penalty.
However, in many areas of suspension mechanics, the suspending background fluid itself is inherently non-Newtonian and may exhibit viscoelastic properties. A sensible step, then, is to extend the Stokesian Dynamics formalism to incorporate a simple model of viscoelasticity. As first seen in Binous and Phillips [2] , we do this by using small spheres as 'beads' in bead-spring dumbbells.
Having done this, we observe the performance of these dumbbell suspensions in simulation by submitting them to oscillatory shear. Oscillatory rheometry, both with small-amplitude shear and large-amplitude shear, has become a standard tool in the classification of viscoelastic fluids [3] . We can compare the measurements from the simulations under oscillatory shear with established constitutive models and experimental results from the literature to establish the validity of this extension of the Stokesian Dynamics method.
Alternative simulation techniques for viscoelastic suspensions have been developed, using finite element [4] , finite volume [5] , and fic-titious domain methods [6] , as well as methods which treat all suspended particles as passive, allowing the fluid flow to be computed separately [7] . These methods require meshing a (normally periodic) domain, which can place computational limits on the size of the domain and the smallest particle size. Stokesian Dynamics does not require meshing, and so can offer larger domains, a wider choice of particle sizes, flexibility with geometry -walls can be created from fixed particles -and different imposed flows.
In Section 2 of this paper, we summarise the theory of oscillatory rheometry and describe how we extend the Stokesian Dynamics formalism to incorporate the bead-spring dumbbell model of viscoelasticity. We then describe how we extract the rheological measurements from our simulations. In Section 3 , we compare the behaviour of different spring laws under small-amplitude oscillatory shear, while in Section 4 , we compare different models under large-amplitude oscillatory shear. In both of these sections, we investigate the effect of altering the parameters in the model, and compare their rheological behaviour to continuum constitutive models.

Theory
Ideal rheometrical measurements are taken in simple shear, or viscometric, flow,  For SAOS, we expect the former, but for LAOS, we expect the latter.
In a linear viscoelastic fluid, the stress response to this applied shear is dictated not just by the current rate of strain, but also by historical rate of strain, The function G ( t ) is the relaxation modulus of the fluid, and represents the importance of the rate of strain from a time t ago on the current stress in the system. Determining the form of the relaxation modulus is the goal of linear rheology, as it allows for the classification of viscoelastic fluids. For example, a purely viscous fluid of viscosity has a relaxation modulus of ( ) = ( ) , where ( t ) is the Dirac delta function, and a linearly elastic solid has a constant relaxation modulus: ( ) = 0 . The relaxation modulus of a fluid can be determined by applying oscillatory shear, where the shear, , and shear rate, ̇, are given by ( ) = 0 sin ( ) , ̇( ) = 0 cos ( ) , for an amplitude 0 and frequency . This can be realised in experiments by placing the fluid in a Couette cell and rotating the inner cylinder so as to impose a shear on the fluid. So long as the gap is narrow compared to the cylinder radii, and there are no instabilities or shear inhomogeneities, this is equivalent to simple shear flow ( Fig. 1 ). In practice the amplitude of the oscillation must be small enough so that the stress response of the fluid is also sinusoidal, i.e. the fluid must remain in its linear regime. At these amplitudes the stress is linear in the amplitude [3] . These tests are called small-amplitude oscillatory shear (SAOS).
If the amplitude is increased, the stress response of a fluid may no longer be sinusoidal. For these large-amplitude oscillatory shear (LAOS) tests, a typical nonlinear response is demonstrated in Fig. 2 . Although the following definitions are only defined for small-amplitude oscillatory shear, their large-amplitude analogues provide useful rheological data [3] , as discussed in Section 4 .
Imposing an oscillatory shear, Eq. (3) , if we stay in the linear regime the stress can be written as [8] where G ′ is the storage modulus and G ′′ the loss modulus. This form is powerful because it splits the viscous and elastic contributions: the storage modulus G ′ is associated with the total shear , and thus repre- sents elasticity. The loss modulus G ′′ is associated with the instantaneous shear rate ̇, and thus represents viscosity. The two moduli G ′ and G ′′ are measured by rheologists as a function of frequency, , for a wide range of viscoelastic fluids. A typical example is shown in Fig. 3 . The inverse of the frequency where the two curves intersect, = 1∕ intersect , is described as the typical relaxation time of the fluid. This parameter allows us to nondimensionalise Eq. (3) [9] , writing the imposed shear as where the Deborah number, = , is the ratio of the relaxation time to the oscillation period, and the Weissenberg number, = 0 , is the ratio of viscous forces to elastic forces. However, since determining requires us to have already determined G ′ and G ′′ , we choose not to nondimensionalise the equations in this way.

System details for simulations
In this paper we perform oscillatory simulations on a sample of Newtonian fluid with dumbbells suspended in it. We implement the dumbbells in Stokesian Dynamics as pairs of small spheres (with radius = 1 in our choice of units), a variable distance apart, with a force between the two. We place our dumbbells in a periodic box of side length 150 units. The dumbbells are constrained to lie in a plane, 2 units deep, so that our simulations are carried out on a monolayer of particles. The method remains three-dimensional.
To stop the dumbbells contracting to zero length in an otherwise quiescent flow, the dumbbells are given a natural length of = 20 , with which they are initialised. We later vary the number of dumbbells in the sample, but unless otherwise stated, we use an area concentration of 10%. The suspension undergoes eight shear periods, with the number of timesteps per period ranging from 80 to 1200, depending on the frequency and amplitude of shear; see discussion in Section 3.2 . The final shear period is used for our analysis; and we use the ensemble average of three independent solutions.
The use of dumbbells in Stokesian Dynamics was first seen in Binous and Phillips [2] , where a relationship was formulated between the dumbbells' velocities and the forces exerted on them. In our implementation for small-amplitude oscillatory shear, we go further by letting beads interact with each other through lubrication, allowing us to examine concentrated suspensions. No non-hydrodynamic forces on the beads were found to be necessary. In our implementation for large-amplitude oscillatory shear, we do not include bead-bead interactions as the high amplitudes at insufficiently small large timesteps give the potential for numerical instabilities. This is discussed further in Section 4.6 .
We also develop the use of Stokesian Dynamics for rheology by introducing a background shear rate, ∞ , to our fluid. This is a convenient way to impose shear flow on the sea of dumbbells through the background shear, The beads are assumed to be torque-free, and so bead interactions are considered only up to the level of forces. Stresslets from the dumbbells are calculated from the forces on the beads.

Extraction of rheological measurements from simulations
To measure the moduli, G ′ and G ′′ , we need to be able to extract the fluid stress from the simulation. In laboratory rheometrical tests, the fluid stress is extracted from the force experienced by the walls of the rheometer. We can measure this stress indirectly by first considering the total force on a volume, V , of fluid, with boundary V , which is given by where is the stress in the fluid and n is the unit normal. We want to measure the force on an imaginary upper wall. Given our sample is a single layer of spheres of radius a , the depth of this upper wall is 2 a . Hence the force on the imaginary upper walls is where L W is the box width. The total stress (per unit monolayer volume), , in a Newtonian fluid with N suspended particles of area fraction c is [10] where we are summing over the stresslets, , on the fluid by the particles, indexed by . Note that c /2 a 3 N is the reciprocal of the fluid volume. The deviatoric stresslet on the fluid by a spherical particle, , at x with a force distribution, f , upon it is given by Treating the dumbbell as two point particles a vector distance Δx apart and with a force difference F , the stresslet on the fluid is given by Calculating the stress in the fluid this way removes the requirement for physical walls on which to measure the stress. We have also performed measurements with fixed walls, and with wall-imposed motion, and found that the results are all in good agreement. In the next section, we show the results of these tests for various spring laws to validate this method as a way to simulate viscoelastic background fluids.
Having measured the stress, we can extract G ′ and G ′′ in at least two ways. For small-amplitude oscillatory shear, the stress response is also sinusoidal, but delayed by a phase shift, , We can extract G ′ and G ′′ from this by using the sine sum identity: ( ) = 0 sin ( ) cos + 0 cos ( ) sin The phase shift, , and amplitude of the stress response, 0 , are easily measured by plotting the stress response against the strain of the system. These plots are Lissajous curves (sometimes Lissajous-Bowditch curves) [11] , and they provide a useful way of describing the viscoelasticity of a fluid. In a linearly elastic solid, = 0 . This corresponds, as can be seen in Fig  the form of a tilted ellipse. In a LAOS test, we once again apply an oscillatory shear, but this time we find that for these large amplitudes, 0 , the stress response is not a pure sine function; instead it is nonlinear, as in Fig. 5 a. With smallamplitude shear, we were able to define G ′ and G ′′ as rheological properties of the fluid we were sampling. The non-sinusoidal stress response also provides a fingerprint for the fluid, giving us not just equivalents to G ′ and G ′′ , but also additional information about the nonlinearity. LAOS tests can therefore be used to provide diagnostic information, allowing us to classify viscoelastic fluids [12] : see Section 4.1 .
We can test a sample to find the extent of the linear regime -to find how large 'large' has to be in LAOS -by performing a G ′ and G ′′ sweep on the shear amplitude for a fixed frequency. In the linear regime, these moduli do not depend on the amplitude of oscillation, but at larger amplitudes, they do. Fig. 6 confirms this for our Hookean dumbbell solution, tracking a fixed frequency over a range of amplitudes. For a more solid definition of 'large', some recent literature has defined LAOS for polymeric liquids as operating for 0 > 1 [9] . In the nondimensional terms introduced in Eq. (5) , the nonlinear regime is reached when the Weissenberg number, Wi , -the ratio of viscous forces to elastic forces -is not small and when the strain amplitude, Wi / De , is also not small [13] .
So far, however, G ′ and G ′′ have only been defined for SAOS, so we first have to find their LAOS analogues. We do this by presenting the stress response as a Fourier series [3] ,  where the first term is exactly the response from linear viscoelasticity, Eq. (12) . Note that the summation is over odd frequencies because the imposed stress is odd. We can split up this expression for ( t ) in the way we did before, to write In the literature, the first order terms ′ 1 and ′′ 1 are those which are typically given: these are the LAOS analogues of G ′ and G ′′ and are often called just that.
We can extract the ′ and ′′ terms from the measured stress by taking its discrete Fourier transform: a technique known as FT-rheology.
An important caveat we have so far neglected is that it is not given that the stress response to large-amplitude oscillatory shear should be periodic -sometimes described as the fluid reaching 'alternance' [9] . As collated by Khair [13] , experiments on polymer melts and colloidal gels have found stress responses to LAOS which are only quasiperiodic [14] or indeed completely aperiodic [15,16] . However, in our numerical experiments, we do find periodic stress responses.

Small-amplitude oscillatory shear
As we perform oscillatory simulations, we show here that different force laws on the dumbbells correspond to different fluid models: some by design, and some by observation.

Newtonian: No dumbbells
The trivial case we should check is the empty Newtonian fluid. As can be seen from Eq. (4) , for a Newtonian fluid, ′ = 0 and ′′ = .
Calculations of these moduli using the stresslet method will certainly agree with these figures, as it follows algebraically from Eqs. (9) and (11) .

Oldroyd-B model: Hooke's law
The simplest force law we can apply to our dumbbells is Hooke's law: for an extension x and spring constant k . The Oldroyd-B model [17] of a viscoelastic fluid is typically represented by a system of equations which govern the evolution of the constitutive equation in time. As it is, in fact, derived from a system of a Newtonian background fluid with bead-and-spring dumbbells, it should provide a theoretical validation of our rheological measurements.
It can be derived from the model [8] that under SAOS, its storage and loss moduli are Here, G is a constant proportional to the spring constant, = , where m is the number of dumbbells per unit volume, and is the relaxation time of the dumbbell, = 3 ∕2 . Since we are interested in the contribution from the dumbbells, rather than the solute, we choose to plot G ′′ without the viscosity component. Along with the shape of the graphs of G ′ and G ′′ , then, this also allows us to compare some extracted constants: G is the asymptote of G ′ as → ∞, and = 1∕ is the crossover point of G ′ and G ′′ .
These theoretical curves are plotted in Fig. 7 (dashed lines) with the constants G and extracted from the simulation (solid lines). There is excellent shape agreement throughout, which is as we would expect. However, to attain this agreement requires very fine timestep resolution in the simulation. At relatively large timesteps (40 timesteps per oscillation, dotted lines), we see good agreement for all G ′ and at low frequencies in G ′′ . However, high frequencies give higher readings of G ′′ than expected. To match the theoretical readings requires much smaller timesteps (1000 per oscillation in Fig. 7 , solid lines). We find measurements at our highest frequency in Fig. 7 converge to the theoretical value with an error of ( ) = 2 . 62 −0 . 866 , for N t timesteps per oscillation.
This sensitivity at high frequencies is because the sample acts almost entirely elastically in this region: the Lissajous plot is almost a (diagonal) straight line. Determining the phase difference between the maximum stress and maximum strain therefore requires very high resolution, even when using the Fourier transform method of finding the moduli.
Effect of concentration. Fig. 8 shows plots of G ′ and G ′′ for dumbbell area concentrations of 10%, 20% and 40%. Doubling the concentration leads to readings for G ′ and G ′′ which are, on average, twice as large. At low concentrations, this is in line with what we expect from Eq. (21) , since G is proportional to the concentration of dumbbells. At higher concentrations, we might expect that the dumbbells would interact in a way so as to not lead to a proportional increase in G ′ and G ′′ . This complication does not appear to happen, at least up to = 40% .  The storage and loss moduli, G ′ and G ′′ , for Hookean dumbbells with spring constants = 1 , 2 and 4. Arrows point in the direction of increasing k . At low frequencies, G ′′ is independent of k , converging to G (dotted line), and G ′ ∼ 1/ k ; but at high frequencies, G ′ , G ′′ ∼ k . Note that the crossover frequency 1/ (vertical dashed line) moves proportionally to the right as k increases. Right: Plotted against Deborah number and scaled by k , the different graphs of G ′ and G ′′ collapse onto each other.
Effect of spring constant. Fig. 9 shows plots of G ′ and G ′′ for spring constants of = 1 , 2 and 4. The dependence on k is a playoff between G ∼ k and ∼ 1/ k . Summarised: This is explained by the right-hand figure of Fig. 9 , where the plots of G ′ / k and G ′′ / k against the Deborah number, = , collapse onto each other. This collapse is to be expected from the definitions of G ′ and G ′′ for an Oldroyd-B fluid, Eq. (21) , nondimensionalised, (describing G ′′ without the viscosity component), recalling that G ∼ k . Furthermore, noticing that ∼ 1/ k , we see that for any constant c , This can be seen in the left-hand figure of Fig. 9 : with the logarithmic axes, doubling the spring constant is equivalent to a shift of the curve to the right by a factor of 2 and a shift upwards by a factor of 2. In other words, the effect of changing the spring constant is reproduced (with appropriate scaling factor) by simply changing the frequency of oscillation.
Effect of natural length of dumbbell. In the Oldroyd-B model, the idealised dumbbells have a natural length set by the strength of the Brownian motion, which we are not including. To stop the dumbbells contracting to zero length, then, we specify a natural length, L , for our dumbbells. Hooke's law then applies to the extension , x , from the natural length, for a dumbbell with end-to-end vector Δx . The impact of this imposed natural dumbbell length on G ′ and G ′′ is less clear. Fig. 10 shows the measurements of G ′ and G ′′ for three natural lengths: 10, 20, and 40 units, where lengths are nondimensionalised so that the dumbbell beads have radius 1. At high frequencies, both G ′ and G ′′ appear to scale as L 2 ; whereas at low frequencies, it is less clear. Some clarity is found by plotting G ′ and G ′′ against the Deborah number, = , where in the right-hand figure of Fig. 10 , they are seen to scale directly as L 2 . This is not altogether surprising: Eq. (11) shows that the dumbbell stresslet contribution, from which the total fluid stress is calculated, behaves as The dumbbell length Δx is dominated in the small-amplitude oscillatory regime by the natural dumbbell length L , and the force is proportional to the extension. The extension due to the background flow also scales as L , since dumbbells beads which are further apart feel a greater velocity gradient. Hence the stresslet behaves as L 2 , which is what we observe. The irregular scaling for low frequencies in the left-hand graph derives from how the relaxation time of the dumbbells, , depends on the natural length. Although this appears to decrease linearly (in this regime) for increasing natural length, the reason for this remains unclear.
Effect of amplitude. The plots of G ′ and G ′′ for oscillation amplitudes of 0.1, 0.2 and 0.4 are indistinguishable from each other. The moduli G ′ and G ′′ are shown to be independent of amplitude, which is what we expect from Eq. (21) . We will see in Section 4 how we can use this as a judge of when we are in the small-amplitude regime.  11. The force response to extension for Hookean, FENE and inverse Langevin (which FENE is an approximation to) dumbbells. SAOS experiments, such as we are doing here, stay in the small-amplitude regime where the difference is negligible.

FENE dumbbells
A limitation of using Hooke's law to model polymers is that real polymers have a maximum length. Furthermore, in extensional flow, Hookean dumbbell models can predict negative viscosities at sufficiently high strain rates [ 18 , p. 44]: a clearly unphysical result. Both these problems are often fixed by using a finitely extensible nonlinear elastic model, or FENE model. This force law is governed by where is the maximum extension of the dumbbell and = | |. More than just being a 'fix', it is a close approximation to a force law derived for polymer deformation from molecular arguments [ 19 , p. 428, eq. B-4]. This law, the inverse Langevin force law, has the form where ̂ is a unit vector in the extension direction, and the factor of 3 allows us to use the same spring constant as the Hookean and FENE models [ 20 , p. 14]. Fig. 11 shows the similarity. Although there are analytically derivable forms of G ′ and G ′′ for FENE dumbbells [21] in small-amplitude shear, as illustrated in Fig. 11 , FENE (and inverse Langevin) dumbbells in SAOS stay in the mostly Hookean regime. This leads to indistinguishable G ′ and G ′′ readings.

Large-amplitude oscillatory shear
Since in the small-amplitude cases, many spring laws reduce to Hooke's law, we continue our investigation by performing largeamplitude oscillatory shear (LAOS) experiments. In most industrial operations, deformation of fluids is large and fast, and therefore outside the linear regime.

Classifying fluids with LAOS
Hyun et al. [12] (with further discussion in Hyun et al. [3] ) outlined four types of behaviour for G ′ and G ′′ in a LAOS amplitude sweep, shown in Fig. 12 , which allow us to classify viscoelastic fluids. These classes are: • Type I (strain thinning): Both the elastic and viscous moduli decrease with increased amplitude. The mechanism behind this is attributed to entanglement of polymer chains. At low amplitudes, the chains remain entangled, but at high amplitudes, the chains untangle as they are stretched, aligning in the flow direction. This reduces the local drag, decreasing the viscous modulus. As the chains extend, they reach their maximum length, reducing the elastic modulus. Most polymer solutions and melts display this behaviour [12] . • Type II (strain hardening): Both the elastic and viscous moduli increase with increased amplitude. The mechanism behind this is attributed to strong interactions between segments of the fluid. As the amplitude of oscillation is increased, these segments interact in a way to increase resistance to flow. For example, polymer chains crossing can form microstructure which inhibits further deformation in the flow direction. Many biological gels -Hyun et al. [3] give the examples F-actin, fibrin and collagen -exhibit this behaviour, as well as PVA/Borax (polyvinyl alcohol/sodium borate) solutions [22] . • Type III (weak strain overshoot): Both the elastic and viscous moduli eventually decrease with increased amplitude, but the viscous modulus initially increases. The mechanism here is a combination of the previous two types: initially weak microstructure is formed which resists flow, but after a critical shear amplitude, the structure is destroyed and the chains align in the flow direction, aiding flow. This feature is common in soft, glassy materials: concentrated emulsions, suspensions and pastes among other examples [3] . • Type IV (strong strain overshoot): Both the elastic and viscous moduli eventually decrease with increased amplitude, although they both initially increase. The mechanism here is the same as Type III, but with stronger intermolecular interaction (but not as strong as in Type II). This behaviour is seen in associative polymer solutions [23] .
This classification of complex fluids looks solely at the behaviour of the first harmonic, which makes it a simple method for experimental use as it can be typically done with standard rheological equipment (using Fourier Transform analysis) without the need to collect full stress data. However, as we have seen, LAOS is typified by the non-sinusoidal stress response: to capture more information about the sample, other classifications exist which measure the contributions of higher harmonics. One such test we will use here is measuring the contribution of the third harmonic relative to the first harmonic, as first described by Hyun and Wilhelm [24] . (The factor of 1∕ 2 0 at the front is due to the tendency of the n th harmonic to be proportional to the n th power of the amplitude, at small amplitudes.) They found that the behaviour of Q ( 0 ) as the shear amplitude increases was different for different fluid samples: for monodisperse linear polystyrene melts, it decreases, but for comb polystyrene melts, it increases. Concluding that Q was somehow related to polymer topology, they went no further. Nonetheless, by its definition, higher values of Q correspond to increased nonlinear effects, and thus we will find it useful as a measurement of these.
A final measurement we can take, given that we have all the particle position data, is the degree of alignment of the dumbbells in the shear direction. Alignment of dumbbells in the flow direction can affect elasticity in that direction, and we can quantify how much by looking at the angle the dumbbells make with the shear direction: we call this the mean dumbbell pitch . There are two options to measure the angle from each dumbbell vector Δx , as shown in Fig. 13 . One option is to measure dumbbells leaning into the shear as distinct from dumbbells leaning away from the shear -we call this the signed mean pitch, where ̂ is the unit vector parallel to the shear, and this gives a pitch of 0°≤ < 180°. In a well-mixed solution, the mean signed dumbbell pitch should be 90°. The other option is measure pitch by how far each dumbbell is away from being horizontal -we call this the unsigned mean pitch, where taking the inverse cosine of the absolute value of the dot product gives the smallest angle between the two vectors, 0°≤ < 90°. In a wellmixed solution, the mean unsigned dumbbell pitch should be 45°. The two measurements, signed and unsigned, show us something different: for example, if we consider two dumbbells with signed pitches of 10°and 170°, their signed mean is 90°, whereas their unsigned mean is 10°. The signed measure shows that there is no left-right bias, whereas the unsigned measure shows that the dumbbells lie almost parallel to the shear direction. Fig. 14 shows the mean dumbbell pitch, signed and unsigned, for Hookean dumbbells over two cycles of shear at different amplitudes. For both measurements, the suspension starts well-mixed. As the suspension is sheared, the mean pitch decreases, as to be expected; larger amplitudes lead to larger decreases. Larger amplitudes also lead to flatter plots at the edge of the shear.  Both subfigures appear to show that the dumbbell pitches return to a well-mixed state at the middle and end of every shear cycle. This indicates that alignment is not permanently altered by the applied shear. This is quantified in Fig. 15 , where we describe alignment as a function of shear amplitude by measuring the peaks or troughs of the unsigned mean dumbbell pitch in Fig. 14 b, following the initial shear cycle. In this figure, we can use the peaks near a shear angle of 3 , when the system is sheared least, or the troughs near a shear angle of 7 /2, when the system is sheared most, to give a good measure for a given amplitude. We see from the slight decrease in the peaks in Fig. 15 that large amplitudes are, in fact, leading to some small alignment bias parallel to the shear direction.
As seen, LAOS tests are good indicators of the microstructure of the fluid, with particular sensitivity to molecular interaction and the shape and entanglement of polymers. Here we shall see the LAOS response to our bead-and-spring suspensions for varying parameters.

System details for LAOS simulations
Recall, from Section 2.2 , that we are considering a periodic sample of the dumbbell suspension, with a square periodic box of side length 150 (recall that in these dimensionless units, the dumbbell beads have radius 1) and 10% area concentration. The natural dumbbell length is set at = 20 . Although in SAOS we considered bead-bead interactions, for LAOS we neglect these, given the potential for numerical instabilities at large amplitudes with insufficiently small timesteps. We run these simulations at 1200 RK4 timesteps per shear period. We discuss this further in Section 4.6 .

Oldroyd-B model: Hooke's law
We start by examining dumbbells with their extension obeying Hooke's law, Eq. (20) . A true Oldroyd-B fluid should have a linear shear stress response at both small and large amplitudes.
Effect of spring constant. Fig. 16 shows the change in some rheometrical properties of our suspension after changing the spring constant, k . This is equivalent to changing the frequency of oscillation, , through the extension of the relationship we found in Section 3.2 for SAOS, Hence this is also equivalent to changing the Deborah number, = . Subfigure (c) shows G ′ and G ′′ , which are normalised in (a) and (b) on the smallest amplitude. Subfigure (e) shows the nonlinearity parameter Q , defined in Eq. (28) . Looking first at subfigure (e), we see that Q increases sharply as the amplitude decreases. This is in contrast to G ′ and G ′′ , which converge at this limit. However, we should not worry: as seen in the analysis of the theoretical fluid models in Section 4.5 , this exact phenomenon is an artefact of the timestep size. The measurement of Q is always worst at the smallest amplitudes, for as we head into the linear regime, the higher harmonics contribute less, and our measurements are therefore less precise. In particular, Q is inversely proportional to the amplitude squared, and this can lead to measurements for 0 ≲ 1 which are less reliable. The theoretical discussion in Section 4.5 shows that decreasing the timestep size by two decades leads to readings for Q which converge to the measurements already seen at the current timestep size for 0 ≈ 1. This observation of Q aside, we also see that increasing the amplitude decreases the nonlinearity, and it is mostly unchanged by changes in the spring constant. Subfigure (f) demonstrates this with Lissajous curves for the high amplitude of 0 = 25 . 6 : all the curves, despite having different 'widths', are propeller-shaped, with the same concavity. The different widths in the curves are physically explained by considering that stronger dumbbell springs require more applied shear to deform them.
The aforementioned concavity is described by Gilbert and Giacomin [25] as 'wrong convexity', which they found for De ≳ 1. Wrong convexity is not observed experimentally [26] .
The mean alignment of the dumbbells is measured in subfigure (d): changing the spring constant is shown to not have any great effect, although there is the slightest demonstration at the highest amplitudes that increasing the spring constant leads to stronger alignment at the moments of maximum shear rate (the peaks).
Looking at the normalised values of G ′ and G ′′ in subfigures (a) and (b), we see that both moduli ultimately decrease at high amplitudes, showing strain thinning, as particles align at the edge of the shear cycle. However, increasing the spring strength takes the suspension from a Type I, where we see mostly only strain thinning, to a hybrid Type I-Type IV, where G ′ sees an initial increase as the strain amplitude is increased, before eventually decreasing. This initial increase in G ′ , the elastic modulus, makes sense: as we increase the amplitude, the elastic forces in the system do more work.
We only observe a 'bump' in G ′ for spring constants k ≳ 1. Above this value, the size of the bump appears to be proportional to the spring constant. Below this value, we just see G ′ decreasing, although there is a turning point at k ≈ 0.25: for passive dumbbells ( = 0 ), we would expect flat G ′ ( 0 ) profiles, so for small enough spring constants, G ′ starts to decrease less. Thus the spring constant has two effects: to alter the size of the bump, and to alter the amount of strain thinning.
Interestingly, the transition stage of a bump in G ′ but not in G ′′ is not one of the four explicit types labelled by Hyun et al. [12] : the bump is discussed in Section 4.6 .
The raw measurements in subfigure Fig. 16 (c) show that the viscous modulus, G ′′ , increases for higher spring constants. However, at low amplitudes, the elastic modulus, G ′ , initially shifts upwards with increasing amplitude before shifting downwards again: readings for k are close to readings for 1/ k . This is only seen at low amplitudes, though: at high amplitudes, the overshoot in G ′ , which grows for increased spring constant, is large enough so that G ′ only shifts upwards for increased spring constant.
Effect of natural length. In Fig. 17 , we see G ′ and G ′′ for dumbbells of different natural lengths, L . Again as in SAOS, since the natural length of the dumbbell is a required physical parameter, but not part of the Oldroyd-B model, the effect of increasing it on the rheological parameters is quite complicated. If we first look at subfigure (d), though, the nonlinearity parameter Q has at least one point to make: at the largest amplitudes, there is not much difference between the samples. This makes sense: at the largest amplitudes, the differences between the natural lengths is comparatively small. For mid-sized amplitudes, shorter natural lengths display less nonlinearity. Looking at subfigure (a), as seen for SAOS, increasing the natural length leads to increases in the readings for G ′ and G ′′ . The normalised data in (b) shows a slightly larger bump in G ′ for higher dumbbell lengths: larger lengths have greater exposure to the shear, which could explain this. The viscous modulus, G ′′ , is unchanged.
The alignment of the three samples, in subfigure (c), are initially different. This difference is maintained at the peaks, but the plots converge at the troughs. The trend of slightly decreasing peak measurements, already identified in previous test cases, is also seen here.
Effect of concentration. We see in Fig. 18 an amplitude sweep for concentrations of = 10% and 20%. Large-amplitude effects begin to kick in at amplitudes of around 0 = 2 . Looking at subfigure (b), both suspensions again lie somewhere between Types I and IV, as we see a upward bump in G ′ before decreasing, yet G ′′ stays constant until decreasing.
As we saw in the frequency sweep, subfigure (a) shows that doubling the concentration leads to readings for G ′ and G ′′ which are, on average, twice as large. This remains true in the nonlinear regime as well, with the normalised plot, (b), showing similar readings for the higher concentration. This makes some sense: we are seeing simply the sum of forces from twice as many dumbbells. Similarly, the measure of dumbbell alignment in (c) shows very little difference between concen-trations, and the nonlinearity reading in (d) also shows little difference in nonlinearity between the two solutions.

FENE dumbbells
In small oscillations, as we have seen in Section 3.3 , dumbbells experiencing a finitely extensible nonlinear elastic force, Eq. (26) , are indistinguishable from Hookean dumbbells, as the applied shear does not cause the spring to extend outside of its mostly linear regime.
Effect of maximum extension. At larger shear amplitudes, we begin to see distinctions between the different FENE models and Hookean dumbbells, although they are slight. Fig. 19 shows the readings for G ′ and G ′′ for four samples: Hookean dumbbells (essentially FENE dumbbells with 'maximum extension' = ∞), and then FENE dumbbells with = 40 , 20 and 10. These values are chosen following preliminary simulations, in which we found that the maximum extension of Hookean springs with spring constant = 2 , undergoing oscillatory shear at the largest amplitude, is no larger than 50. In each case, the Deborah number is the same.
Until 0 ≳ 10, there is not a great difference in the G ′ and G ′′ plots. After that, we see that the shorter maximum extensions break away from the Hookean measurements, with G ′ decreasing faster for shorter maximum extensions, and G ′′ decreasing faster for larger maximum exten- sions. The nonlinearity parameter, Q , displays greater nonlinearity at the smaller maximum extensions, while the mean dumbbell pitch is mostly unchanged.

LAOS with fluid models
Several fluid models are derived from idealised dumbbell suspensions. We now compare the analytically-derived large-amplitude behaviour of these models with our LAOS simulations. For Hookean dumbbells, the appropriate model is the aforementioned Oldroyd-B model, and for FENE dumbbells we have the FENE-P and FENE-CR models.
The rheological behaviour of these models at large amplitudes can be derived from their governing equations, which for all three have the form -respectively, the continuity and momentum equations, together making up the Navier-Stokes equations -and the constitutive equation, where the evolution of the elastic or polymeric contribution to the stress tensor, , is given by The parameter is the relaxation time of the dumbbells, as introduced in Section 3.2 , and the functions f ( R ) and g ( R ) depend on which model we pick. In all cases, 2 = tr ( ) . The parameter G is proportional to the spring strength and dumbbell concentration. In the constitutive equation, is once again the symmetric part of the velocity gradient tensor, and is the identity tensor.
It is the extra contribution to the stress from the dumbbells, the polymeric stress = ( ) , which makes this different to the Newtonian fluid model: indeed, if ( ) = 0 , Eqs. (32a) to (32c) are the governing equations for a Newtonian fluid.
The evolution of this dumbbell stress contribution is given in Eq. (32d) , where the upper-convected time derivative of this tensor is given by Derived in Oldroyd [17] , this is the material derivative appropriate for a line element which rotate and stretches with the fluid. That is to say, for a infinitesimal line element d , advecting passively in the flow, the derived tensor d d has The dumbbell stress contribution comes from the failure of the dumbbells to deform with the fluid elements.
The functions f ( R ) and g ( R ) for each model are given by [27] Newtonian: Oldroyd-B: FENE-P: FENE-CR: where is the maximum dumbbell extension.
Suppose we now apply oscillatory xy -shear, Eqs. (1) and (3) . To find the rheological parameters, we need to find the xy term of the stress tensor, xy , which by Eq. (32c) depends on A xy and R ( = √ + ). The evolution equation of , Eq. (32d) , can be rearranged and written out as We therefore have to solve this coupled system of evolution equations. In the Oldroyd-B case, when ( ) = ( ) = 1 , we can solve this system explicitly [8] . Otherwise we can do so using a numerical ODE solver.
Taking the solutions for A xx , A xy , A yy , then, we construct the solution for xy and pass this stress through the same process described at the beginning of this section to find G ′ , G ′′ , their higher harmonics and Q , for a chosen frequency, as functions of amplitude.
Oldroyd-B fluid. The polymeric stress, , is derived from the forces on the dumbbells, which in the models come from Stokes drag, Brownian motion and the spring force, F . If the dumbbell length is given by R , the polymeric stress is given by where the angle brackets ⟨ · ⟩ represent some form of averaging over all the dumbbells. We can write the Hookean force law, Eq. (20) , as = . Therefore the polymeric stress is given by For a true Oldroyd-B fluid, then, G ′ ( 0 ) and G ′′ ( 0 ) should be constant while, Q ( 0 ) should be zero throughout, given that the higher harmonics ′ 3 and ′′ 3 are zero. Given that our simulations with Hookean dumbbells show nontrivial rheological parameters, these must be affected by something else.
FENE-P fluid. We can write the FENE force laws for the dumbbells, Eq. (26) , as The force law is no longer linear, making us unable to describe its evolution in a closed form. Approximations are therefore needed.
In the FENE-P model, developed by Bird et al. [28] after inspiration from Peterlin [29] , the force function f ( R ) is replaced by the force due  to the average dumbbell length, f ( ⟨R ⟩). Hence the force law is linearised as and the polymeric stress contribution becomes where 2 = tr ( ) . This evolves as in Eq. (32d) with ( ) = ∕ ( ) , and gives a shear-thinning viscosity.
The resultant rheometric plots are seen in Fig. 20 . Subfigure (a) shows both G ′ and G ′′ decreasing, in the style of a Type I fluid. A comparison with the simulation measurements can be seen in Fig. 21 a. The decrease in G ′′ matches well, similar to what we see in our simulations, but the behaviour of G ′ at moderate amplitudes is different: the model predicts no bump in this regime, although it does decrease at a similar rate following this.
The Lissajous curve in Fig. 20 b shows the evolution from a linear response at small amplitudes, to a flattened ellipse, to a flag-shape at the highest amplitudes.
The nonlinearity parameter, Q , in Fig. 20 c, decreases for increased amplitude, but we can see that the results at low amplitudes are very sensitive to timestep size. These are interesting when compared to the simulation measurements in Fig. 21 b. In the comparison figure, Q with the FENE-P model matches the simulation well at the highest ampli-tudes. However, the overall shape observed in the simulation, with the dip at moderate amplitudes, is similar to the overall shape for the FENE-P model at low resolutions. Increasing the resolution gives convergence to a flat profile for low amplitudes, which is what we would expect, given the convergence of G ′ and G ′′ here. This confirms that what we see for Q in the simulation at low amplitudes is noise, and nothing more substantial.
Changing the dumbbell concentration or spring law (and hence G ) simply scales the stress response, in agreement with the experiments in Fig. 18 , where we change the concentration, but in contrast to Fig. 16 , where we change the spring constant. In the latter case we saw higher spring constants leading to larger G ′ bumps, and the absence of this here is the largest difference between what we see in simulations and the FENE-P model.
Changing the frequency or maximum extension in Fig. 22 leads to more interesting regimes. Increasing the frequency leads to faster decreases in G ′ but the creation of a bump in G ′′ : the Type III regime. This is the opposite of what we see in the simulations, Fig. 16 , where increasing the frequency -equivalent to reducing the spring constant -reduces both G ′ and G ′′ . Increasing the maximum extension makes the dumbbells more Hookean, leading to a more linear response. For G ′′ , this agrees with our FENE suspensions, Fig. 19 , although again we do not see the bump in G ′ .

FENE-CR fluid.
The FENE-CR model, developed by Chilcott and Rallison [30] , makes an artificial change to the evolution of , Eq. (32d) , setting ( ) = 1 . This keeps the steady shear viscosity constant. Fig. 23 shows the rheometrical measurements which come from solving the governing equations. In subfigure (a), G ′ decreases with increased amplitude, whereas G ′′ rises slowly. The Lissajous plot in subfigure (b) shows the form of the nonlinear response at higher amplitudes, morphing from an ellipse to a rounded rectangle, before experiencing kinks towards the edge of shear. The nonlinearity parameter, Q , in subfigure (c) shows a familiar pattern, again decreasing at higher amplitudes.
Looking back at the comparison in Fig. 21 a, this fits our simulation measurement worse, given that G ′′ does not decrease. Of course, this is to be expected: the simulations show the suspension to be shear-thinning and this model is designed explicitly to have a constant viscosity.
The response of the FENE-CR model to changes in spring constant, concentration, frequency and maximum dumbbell extension is mostly similar to the FENE-P model: there is still a bump in G ′′ at higher frequencies.

Comparison of simulations and fluid models
Despite tentative matching of G ′′ between our simulations and the FENE-P model, none of the continuum models show the bump in G ′ that we see in our simulations. The bump is a robust feature at high spring constants which we also observe in fully 3D simulations of our suspensions.
Furthermore, the bump is also found in both the presence and absence of hydrodynamic interactions between the particles. In our LAOS simulations we ignore hydrodynamic interactions between beads since we are primarily interested in the dumbbell beads as markers, and also since insufficiently small timesteps can cause numerical instabilities in Stokesian Dynamics. Fig. 25 verifies the suitability of this approximation, demonstrating that G ′ and G ′′ amplitude sweeps where long-range Fig. 25. Amplitude sweep for G ′ and G ′′ for Hookean dumbbells with = 0 . 56 , where particles interact through far-field hydrodynamics using the Rotne-Prager-Yamakawa tensor. The G ′ bump is visible. Sweep undertaken at reduced resolution of 120 RK4 timesteps per oscillation, compared to the usual 1200, due to increased computational intensity.
hydrodynamic interactions are included, here using an Rotne-Prager-Yamakawa potential [31,32] when particles overlap, also display the G ′ bump.
We observed first in Fig. 16 that the Lissajous curves of the sheared suspension display 'wrong convexity' as described in Gilbert and Giacomin [25] , i.e., that there is an inflection point in the upper branch of the curve occupying the first quadrant. We see from Fig. 26 a that at high amplitudes, the Lissajous curve has a propellor shape, displaying stationary points near the neutral shear position which are not seen in real-world experiments [26] . In a shear cycle of the simulations, as the suspension is initially sheared, the dumbbells stretch but resist this stretch by pulling their endpoints together under the spring force. As the suspension is then returned to its neutral position, the dumbbells are compressed, and given their natural length, the dumbbells resist this compression by pushing their endpoints apart. The strength of this 'push' increases with spring constant, not just because of the increased strength of the spring force, but also because stronger spring dumbbells have a smaller extension at the extreme of the shear cycle, just before compression begins.
This pushing behaviour distinguishes our simulations from the fluid models. We can remove the stresslets derived from pushing dumbbells to produce the Lissajous curve in Fig. 26 b, However, this is unlikely to be the cause of the G ′ bump, since we find that at all amplitudes, the value of the shear for which the maximum stress is obtained, and the value of the stress when the maximum shear is obtained, both of which define G ′ (recall Eq. (15) ), remain mostly the same. Our Hookean springs are given a natural length since alternative Hookean spring models, for example giving the dumbbells a zero natural length, or keeping a nonzero natural length but specifying that below this natural length the dumbbell produces no internal force, produce dumbbells which contract to zero length quickly, even in strong shear flow. Fig. 27 shows the resultant stress in the fluid for zero-natural length dumbbells over two units of shear, where the dumbbells are placed in with extensions of Δ = 20 . Even at low spring strengths, we find the same pattern: dumbbells contract quickly to zero length, reducing the stress in the fluid to almost zero. In particular, the stress response is no longer periodic, and we cannot describe it using G ′ and G ′′ .
Again recalling Fig. 16 , wrong convexity is seen for all spring constants, whereas the G ′ bump is only seen for spring forces k ≳ 1: it is therefore unlikely that wrong convexity causes the G ′ bump. With smaller spring forces, the behaviour of G ′ and G ′′ matches the FENE-P model better. Simply, we may not be finding agreement with our analytical models because the spring forces we are applying may be unphysically large.

Conclusions
By placing simple bead-and-spring dumbbells, with various force laws, into a suspension, we are able to show that the resulting suspension behaves as a viscoelastic fluid. We have shown that we can tune the degree of viscoelasticity by tuning certain parameters in these force laws.
In Section 2 , we have developed Stokesian Dynamics to allow us to include these force laws easily, and then been able to measure the Fig. 27. Total system stress over the first unit of shear, where the Hookean dumbbells have zero natural length, starting with extensions of Δ = 20 . Here the spring force is = 2 . In particular, notice that the stress response is no longer periodic. rheological behaviour of the simulated suspensions. We have been able to do so using two techniques, neither of which require the addition of physical walls to be placed into the system.
Our small-amplitude tests in Section 3 confirmed that the dumbbell suspension behaved as an Oldroyd-B fluid, and we showed the effect of changing the parameters in the model. Our large-amplitude tests in Section 4 showed that in the absence of any interparticle forces, other than the spring laws, the dumbbell solution behaves as a strain thinning fluid. With stronger spring constants, the solution exhibits a strong strain overshoot. We were also able to quantify the degree of nonlinearity in our suspension, introducing a measure of dumbbell alignment, and how this depends on the parameters in the spring law. In particular, we have shown that in all cases, increasing the shear amplitude led to decreased nonlinearity. We also compared the measurements from our dumbbell suspensions with established non-Newtonian, constitutive laws, showing the parameter values at which we find agreement and disagreement.
This lays good groundwork for further studies using Stokesian Dynamics, confirming that if we added larger particles, we would be able to model a suspension in a viscoelastic fluid.