Interaction between a free-falling sphere and structure dynamics in a heterogeneous thixotropic fluid

The impact of thixotropy on the settling behavior of a solid sphere is investigated utilizing a ﬁnite element-computational ﬂuid dynamics simulation. Flow behavior is evaluated by coupling the Navier–Stokes equations with the dynamic evolution of an initially heterogeneous ﬂu-id’s microstructure. Studying the structure dynamics around the settling sphere allows us to identify a variety of irregular and linear settling regimes. Settling regimes are varied by the degree of structuring, the degree of associated heterogeneity, the local morphology of the heterogeneous microstructure, and the stress induced by the sphere. In addition, the settling velocity proﬁle of the relatively light spheres temporarily ﬂuctuates in a case where the settling time of the sphere is long enough to capture the local heterogeneity. Ultimately, we compare the results of the simulation of dropping spheres with those of the numerical simulation of different rheological tests. This illustrates that the competition between kernels of orthokinetic and perikinetic build-up and shear-induced break-down of the microstructure indeed allows an understanding of the connection between the ﬂuids’ ﬂow curve and the settling behaviors. Furthermore, settling regimes are characterized based on the rates of build-up and break-down of the microstructure. Moreover, the loss of fore-aft symmetry is observed in the ﬂow ﬁeld around the settling sphere as a result of a viscosity gradient behind and ahead of the sphere.


I. INTRODUCTION
Thixotropic fluids are ubiquitous, ranging from biological substances 1,2 including human blood, tissues, and the cytoplasm of cells, to some natural materials such as waxy crude oils, 3 some varieties of honey, 4 and some types of clay, 5 to several synthetic products such as fumed silica, 6 Laponite, 7 and cellulose nano-fiber (CNF) suspensions. 8,9They consist of an internal microstructure immersed in a liquid medium. 104][15][16][17][18] However, even after a century of studies on thixotropy, it still remains a very challenging topic to understand. 12,13,17xperimentally, depending on their internal structure, thixotropic fluids show different levels of elasticity varying from thixotropic yield-stress fluids (TYSFs), where thixotropy is more pronounced, to thixotropic-elastoviscoplastic fluids (TEVPs), where elasticity is more important. 13,19The principal difference is the volume fraction of the discrete phase.Typically, at high discrete phase volume fractions, the elasticity is enhanced, while the thixotropic time-scales are relatively short (e.g., Carbopol gel) compared to the systems with low discrete phase volume fractions (e.g., Laponite).Hence, the definition of thixotropy is already complicated, since a wide range of materials show some relatively weak degree of recoverable viscosity. 12,17,206][27] The complexity of thixotropic fluids is reflected by their modeling, which is difficult due to the lack of comprehensive datasets including a variety of rheological protocols. 17hese are required to address the multiple length/time scales in the evolution of the microstructure. 12,28Such complexity renders the topic an appealing subject of research.
The microstructure of a thixotropic fluid dynamically evolves due to the inter-particle network being formed by the flow-independent Brownian motion and flow-induced growth and breakage of flowdriven inter-particle network. 13Depending on the imposed agitation, the competition between these mechanisms divides the material into yielded (fluid-like) and un-yielded (solid-like) regions, 22,23,29 where, at each region, the associated rates of these processes are actually a function of the microscopic level properties at the time of observation. 13,16,30ccordingly, the transition point from the yielded to the unyielded region (yield stress), elastic storage modulus G 0 , and viscous loss modulus G 00 are dependent on the state of the microstructure. 16,17,22ondensing the microstructure dynamics into a single constitutive equation considering all rheological features applicable to all thixotropic fluids is known to be very complicated. 13,16,17urthermore, even if such an equation would be available, developing a suitable simulation framework to examine this detailed model against some experimental data, such as dropping solid objects in a thixotropic fluid, would be numerically challenging.6][47][48] In the SKM, the kinetics of forming-breaking of "bonds" in the continuous phase (fluid) at a macroscopic level is modeled.In the PBM, the dynamics of the "aggregation of clusters" of the discrete phase forming the microstructure at a microscopic level and considering the relevant physical processes are formulated. 16In the FCM, the dynamics of the total stress considering both elastic and viscous components is described, where a fluidity parameter stands to determine the influence of the state of the microstructure. 45,48mmersing solid objects is a well-known method of probing the thixotropic fluids.Particularly, this is an approach to investigate the thixotropic phenomena at the local rather than global level.Numerous experimental and numerical studies have highlighted the different settling behaviors observed in thixotropic fluids, as we briefly review in the following paragraphs.
The gravity-induced stress causes a confined sheared area around the solid object. 29,49,50The size of this fluidized region varies with the elastic and thixotropic relaxation times 50,51 and is associated with the drag force. 50In addition, the size of the sheared region ahead of the solid object is very thin for thixotropic yield stress fluids in comparison to simple yield stress fluids. 49Behind the solid object, however, the sheared region characteristics are dependent on the fluid's shear thinning properties. 52haracterizing the flow field around the settling solid object 48,50,53 in TEVP fluids indeed reveals the complexity of the drag force and yield stress. 51,53,54The onset of fore-aft asymmetry of the velocity field and the associated negative wake is believed to be a fundamental feature of elasticity, not thixotropy. 51,53Hence, considering thixotropy, the size of the negative wake increases with aging time. 53oreover, fore-aft asymmetry is reported in numerical simulations that used a structural kinetic model for an inelastic thixotropic fluid. 18,50In a specific case of settling of double spherical solid objects, the fore-aft asymmetry decreases by the distance of the two spheres 55 and by the Peclet number. 18he drag force experienced by the sphere in thixotropic yield stress fluids is computed utilizing the surface integral of the total stress in a direction normal to the yield surface around the solid object. 49hus, the drag force varies with the geometry of the yielded region around the solid object, the smoothness of the solid-fluid interface, and the boundary condition at the walls of the fluid container. 51,56,57n practice, the smoothness is determined in comparison to the size of the discrete phase particles in the fluid, whereas in numerical simulation, it is based on the degree of orthogonality of the mesh lines. 51,56ccordingly, the rough interface implies that fluid velocity has a tangential component equal to the velocity of the solid object and a normal component equal to zero, i.e., no-slip boundary condition. 51In contrast, the smooth interface implies that the velocity of the fluid is proportional to the normal derivative of the velocity at the interface wall where the slip length controls the proportionality, i.e., Navier-slip boundary conditions. 51In addition, the drag force increases due to the wall effect of the container.To avoid this, the size of the width of the container needs to be greater than five times the size of the solid object, 49,57 and in the numerical simulation, no-slip boundary condition needs to be applied on the walls.
The trajectory of the solid object in thixotropic fluids follows unexpected settling velocity profiles.These settling regimes are characterized based on the gravity-induced stress imposed by the solid object and the aging time of the fluid. 9,29,48,58,59Among them, the irregular semi-oscillatory and the continuously decreasing velocity regimes are recurrent.A continuous decrease in the solid object velocity, only reported for thixotropic fluids, is believed to occur due to a continuous restructuring in front of the solid object. 9,29In addition, the irregular settling regime in the presence of thixotropy is believed to be due to shear-banding, 59 developing heterogeneous flocs, 9 and other instabilities and/or rheochaos in the flow field behind the solid object. 48In a porous non-thixotropic yield stress fluid, this regime is believed to be caused by the solid object capturing the porosity of the fluid. 60hixotropy also influences the occurrence of other settling regimes such as steady-settling and no-motion regimes.In this matter, the rate of microstructure build-up leading to the increase in the "dynamic yield stress" determines whether the settling velocity reaches its terminal value or ends up in stoppage. 61on-homogeneous flow, where the fluid viscosity varies spatially, is another characteristic feature of thixotropic fluids that could affect their interaction with falling solid objects.Such spatial heterogeneity occurs due to different reasons.Indeed, in practical applications where flow geometries induce a spatial variation of the local shear stress, the fluid experiences locally a distinct shear history. 18Since thixotropic suspensions are known to constitute flocs, aggregates, and colloidal particles, consequently such discrete (porous) elements surrounded by fluid make them naturally prone to flow instabilities.. 6,13,16 Under shear, unstable colloidal suspensions exhibit the development of a locally heterogeneous microstructure due to the continuous breaking and re-aggregating of particle networks. 62In addition, heterogeneity by aging is also observed.4][65] Apart from flow-driven or aging-related heterogeneity, in some thixotropic gels such as in CNF gels, the particle size distribution is extremely wide, which causes spatially non-homogeneous thixotropic evolution that eventually leads to the development of flocs. 8,9,66,67ollowing our previous experimental work on settling steel spheres in TEMPO-CNF gels, 9 this work targets at a numerical investigation of the impact of thixotropy in the settling behavior of solid objects.There, through rheological measurements, a non-trivial thixotropic, an age-dependent yielding point, and elastic behavior for TEMPO-CNF gels were disclosed. 68In addition, in the dropping of the solid spherical objects, we observed different types of irregularities (depending on the aging time of the sample) and logarithmic behaviors (mostly in rejuvenated samples). 68Considering the wide size distribution of CNFs, we concluded that microstructure dynamics of TEMPO-CNF gels are most likely different from other low-density gels, and consequently, the aging effect that occurs due to the restructuring of fibrous elements by Brownian diffusion does not develop homogeneously across the whole sample. 9,69Finally, we offered that such specific structural dynamics eventually could lead to the observation of irregular and continuously decreasing settling (logarithmic behavior) regimes. 68However, to add evidence to this hypothesis, a well-defined simulation framework and a realistic rheological model constitute a possible approach.
Accordingly, our first aim in this research is to predict the condition of different experimentally reported settling regimes in TEMPO-CNF gels.We consider thixotropy isolated from elasticity to simplify the problem.This is justified since in low-density gels such as CNF gels or Laponite, the elastic effects are limited to relatively small deformations, and the viscous flow is more prevailing nearby the sphere. 59,68,702][73][74] In addition, to investigate the impact of possible spatial heterogeneity, simulation is specifically performed considering a degree of non-homogeneity associated with the initial microstructure.Furthermore, even though it is well-known that rheological properties dictate the settling behavior, finding an accurate link between settling and rheological data has proven to be difficult. 53Therefore, our second aim is to find a relationship between gravity-induced flow and rheological tests, such as flow curves, i.e., explaining how the material's rheology manifests in the sphere drop experiment.Moreover, the breaking of fore-aft symmetry of the flow-field could be due to either elasticity 51,55 or thixotropy; 18,50 hence, the difference is not fully understood.In this regard, our third aim is to characterize the flow field around the falling sphere and predict the condition that thixotropy changes the geometry of the flow field.
The rest of the article is organized as follows: In Sec.II, the numerical method and rheological model are explained.In Sec.III, our observations are reported and discussed, which constitutes three subsections related to the targets of this manuscript.Finally, concluding remarks are provided in Sec.IV.

II. MODEL AND NUMERICAL METHOD
A computational fluid dynamics (CFD) simulation utilizing commercial software COMSOL Multiphysics 6 is used to study the settling of solid spheres in an inelastic thixotropic fluid.To decrease the computational burden, a single-phase fluid flow in a two-dimensional axisymmetric domain is simulated.Figure 1 shows the domain X and a rectangle of height H and width W=2.A half-circular spherical object of diameter d s is removed from this domain, which places the sphere initially at a distance H 0 from the bottom.Table I displays domain dimensions and fluid and sphere parameters.
We used the no-slip boundary condition on the bottom (z ¼ 0) and right (r ¼ W=2) walls of the domain and on the boundary of the FIG. 1.An axisymmetric numerical domain X is shown.The rectangular domain of height H and width W =2 contains half hole of diameter d s , which represents the solid sphere.This is initially placed at a distance H 0 from the bottom, and the total force exerted on it is computed on its surface @X s (red boundary).A refined grid is used around the sphere over an extension of 3d s (region A).

Physics of
half-circle describing the spherical object @X s .The axial symmetry boundary condition is applied on the left wall of the domain (r The open boundary condition is used on the top wall (z ¼ H) of the domain.
The finite element method (FEM) is employed to numerically solve Navier-Stokes equations coupled with partial differential equations (PDEs) describing the fluid-sphere interactions.These sets of PDEs include the constitutive equations of motion of the fluid and sphere combined with the rheological model for the evolving microstructure of the fluid.
In the FEM framework, the (P 2 þ P 1 ) scheme is applied for the discretization of the fluid's velocity components with second-order elements and the pressure field with linear elements.A second-order (quadratic) Lagrangian function discretizes the PDE of the rheological model.Additionally, the backward differentiation formula (BDF) method is used for time stepping.Furthermore, the PDE system is combined with a grid generation scheme that re-meshes when the sphere interface undergoes large deformations.The automatic remeshing is applied to the whole domain based on the degree of distortion of the elements.To adequately solve the PDE system, the mesh is refined around the moving solid-fluid interface as Fig. 1 shows.The validation of the numerical method is detailed in the supplementary material.
Comparisons between the flow curve and the CFD-based settling simulations are enabled by performing numerical experiments on the rheological model using a Python-based script for its integration.This allows for computing the model's responses under different rheological tests, such as flow curves using different point times (times probed at a constant shear rate or stress), and performing transient simulations under constant shear rate or stress conditions.

A. Solid-fluid equations of motion
We employed the incompressible Navier-Stokes equations in cylindrical coordinates (r, z) axisymmetric around r ¼ 0. The continuity and momentum equations are as follows: (1) where U is the fluid velocity, and T ¼ 2gð/Þ_ c is the viscous stress tensor.Here, _ c is the rate of deformation tensor defined as _ c ¼ ðrU þ rU T Þ=2.The half-sphere boundary @X s (red boundary on Fig. 1) follows Newton's second law for a rigid body motion with the net force given by where v s is the sphere velocity, and the second term on the right hand is the drag force F d . 49part from the equations described above, to efficiently account for the movement of the sphere through the fluid, we use an arbitrary Lagrangian-Eulerian (ALE) method with a mesh updating scheme.The movement of the sphere interface is then defined by the kinematic condition

B. Rheological model
6][77] Coupled with the Navier-Stokes equations, it follows in which, the fluid's microstructure (/) evolves by three contending kernels that are 1.K B ðgð/ÞÞ, the build-up of the microstructure due to Brownian motion, which reads 2. K sha ð_ cÞ, the flow-induced aggregation kernel, which reads 3. K shf ð/; _ cÞ, the break-down of the microstructure due to flowdriven fragmentation kernel that reads Here, A s and B s are constant relaxation rates of restructuring and disintegration, respectively.n a , n f , and _ c 0 determine the volume fraction's sensitivity toward shearing.A b is a constant rate determining shearindependent restructuring (due to Brownian motion).The Brownian motion is assumed to be a power-law function of the viscosity gð/Þ in which n b determines the slow-down rate of the motion with the increase of viscosity.This kernel, therefore, models the development of the microstructure during aging.
The viscosity follows Quemada's constitutive equation, 10,78 which reads where / max is the maximum effective volume fraction.g 0 is the limit viscosity when / approaches the minimum effective volume fraction at high shear rates (that is zero here), corresponding to the limit of vanishing aggregates, and n g is a fitting parameter.To analyze the model, we use the thixotropic number N thix , which evaluates the strength of the thixotropic response of the imposed shear stress or shear rate. 13,32oreover, to consider a non-homogeneous fluid structure, a random Gaussian distribution of initial effective volume fractions of / 0 is defined.Depending on the state of the initial microstructure to be examined, N f flocs with / 0 < / max are randomly allocated.To consider a heterogeneous microstructure, we assumed that / 0 is between / 0Min and / 0Max < / max .Having N f distinguish values of / 0 resembles the discretization of the microstructure.In addition, the location of each floc is defined by two-dimensional coordinate (r, z), where 0 < r < W 2 and 0 < z < H.Then, this distribution, which is initially produced utilizing a Python script, is included in the rectangular domain X (Fig. 1) in the COMSOL code.This allows us to create an initial spatially heterogeneous microstructure characterized by three parameters: (1) the initial degree of structuring, / 0 / max that is the average distribution of the initial effective volume fraction in the whole domain scaled to the maximum volume fraction where The geometrical configuration of the initial microstructure in the whole domain that is dependent on the random generator seed value, N f the number of allocated / 0 with coordination (r, z), and the ratio of the length scale and the size of the sphere ' h ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi . In addition, to more accurately examine the impact of the geometrical configuration, we defined a regular distribution of flocs only in front of the sphere where N f is the number of regular flocs and the dimensionless variable ' h ¼ H N f ds .Finally, the scalar version of the structural model described in Eqs.(5) (without U Á r/), (6), ( 7), (8), and ( 9) is numerically solved (using the Python-based script) for constant shear rate or stress to mimic experimental rheological tests.This implementation involves another set of parameters / 0 ; _ c, and s representing the initial effective volume fraction and the scalar values of the shear stresses and shear rates.The selection of the model parameters used in this study is detailed in the supplementary material.

III. RESULTS AND DISCUSSION A. Settling velocity regimes
We ran 45 CFD cases utilizing the set of parameters in Table II.They have the same geometrical configuration for the spatial distribution of the initial microstructure.In each case, / 0 / max and Stdð/ 0 Þ are varied.Then, different spheres with size d s and density difference Dq with zero initial velocity are dropped into the fluid.In addition to these cases, we conducted others in which only the geometric configuration of the initial spatial heterogeneity changed.
As a result, for high enough gravity-induced stress, s s ¼ gd s Dq, where the sphere is able to penetrate the fluid, various irregular or linear regimes are observed.Depending on the combination of the fluid and sphere characters, the overall velocity magnitude of the sphere either reaches a terminal value after a short time or continuously reduces or even increases over time.In contrast, at low s s , the sphere fails to break down the microstructure and penetrate the fluid.Accordingly, settling regimes are classified as The no-motion regime is defined as a case where the sphere travels a distance of ' z < 0:2H, while the settling time for this displacement is relatively long, t s > 800 s.Then it reaches stoppage.The steady settling regime implies a case where the sphere moves with a terminal velocity after it travels a distance of ' z < 0:2H.In the irregular steady settling regime, terminal velocity is fluctuating around an average value, while, in the linear steady settling, the terminal velocity is almost constant.In the continuously decreasing velocity regime, the velocity of the sphere becomes logarithmic and slowed down over time.The variation of the velocity can be accompanied by some fluctuations or be linear.In this regime, the velocity of the sphere finally reaches either stoppage or a terminal velocity.The last regime is the rapid fall, where the velocity of the sphere continuously increases over time.Similar to other regimes, the velocity of the sphere either varies linearly or irregularly.In addition, it is also possible that the settling velocity profile reaches a terminal value.

Gravity-induced stress
Gravity-induced stress s s is a crucial factor governing the settling behavior.Due to microstructure, imposed stresses less than a threshold value prevent lighter spheres from penetrating the fluid.This is due to the yield stress associated with the microstructure where the gravityinduced stress needs to overcome the local yield stress over time.At the far extreme, a heavy sphere forces its way through the fluid, while its trajectory does not vary by the local spatial heterogeneity.Then, the velocity profile of a sphere that induces the fluid s s between these two extremes is in steady settling or continuously decreasing velocity regimes.In this range, the spatial heterogeneity of the microstructure could be captured by the sphere; thus, the settling velocity would exhibit some fluctuation.Figure 2 shows example cases of different settling regimes obtained from imposing different values of s s .In Fig. 2(a), a quick continuously decreasing velocity regime (leading to stoppage) without any fluctuation is observed.This is a result of a sphere of d s ¼ 3:5 mm and Dq ¼ 143 kg=m 3 , imposing the stress of

Degree of structuring and spatial heterogeneity
In general, the velocity magnitude of the sphere increases with reducing the degree of structuring / 0 .Apart from this, a no-motion regime turns into either a steady settling or a continuously decreasing velocity regime by going away from / max .For example, a sphere with d s ¼ 3:5 mm and Dq ¼ 6893 kg=m 3  regime [Fig.3(f)].These results are comparable the variation of the settling velocity and regimes by aging time. 9he manner in which the degree of spatial heterogeneity influences the settling behavior, however, is dependent on the settling time of the sphere t s .In addition, the setting time increases by increasing the degree of structuring while keeping the sphere parameters constant.Figures 3(e), 3(f), and 3(d) show three cases with roughly the same t s , each of them having a distinct degree of heterogeneity Stdð/ 0 Þ of 0.03, 0.012, and 0.0006, respectively.These show that increasing the degree of spatial heterogeneity increases the fluctuations in the settling velocity profile.In contrast, cases with similar parameters, having longer settling time, such as the one in Fig. 3(a), where t s ¼ 1367 s and Stdð/ 0 Þ ¼ 0:0008 [compare to Fig. 3(f) with t s ¼ 0:94 s and Stdð/ 0 Þ ¼ 0:012], offer more temporal fluctuations in the settling velocity.

Geometrical configuration of the spatially heterogeneous microstructure
For a given degree of structuring and spatial heterogeneity, we described different geometrical configurations of the initial microstructure.Figures 4(a)-4(d) show different geometries, where / 0 / max and Stdð/ 0 Þ are constant for all of them.In Fig. 4(f), the impact of each geometrical configuration on the settling velocity of the sphere is shown.As a result, the overall settling behavior follows a continuously decreasing velocity path, and the range of magnitude settling velocity is the same.However, the local variations in each geometrical configuration appear only as fluctuations in the temporal settling velocity profile.Indeed, the velocity profile is more irregular by decreasing the ratio of the length scale and the size of the sphere ' h .To more accurately characterize this, we consider regularly distributed flocs only in front of the sphere [similar to Fig. 4(d)].As can be seen in Fig. 4(f), increasing the number of flocs N f and reducing ' h induces more linear behavior.These results also imply that the sphere only interacts with the spatial heterogeneity in the regions it travels, which is in agreement with the experimental observation on a confined shear region near the sphere. 29

B. Discussion on settling regimes
We estimated the drag coefficient C s 18,51,71 and Reynolds number Re as follows: Here, kF d k and kv s k are the temporally averaged magnitude of drag force [Eq.( 3)] and velocity of the sphere, respectively.Average fields of _ c and gð/Þ are computed on the A around the sphere (Fig. 1) and over N t settling time steps (the average of any flow-related field such as X is computed as . Consequently, C s reduces by Re and a drag correlation of C s ¼ 24Re À0:79 þ 5000Re À0:2 describes the variation adequately (Fig. 5). 28The low range of Re in Fig. 5 signifies the creeping fluid flow around the settling sphere. 79In addition, a high range of C s implies a sluggish motion for the sphere in such a thixotropic fluid.
Considering the kinematic Eq. ( 5), it is expected that far away from the sphere, where the shear rate goes to zero, the Brownian kernel is dominant and builds up the microstructure.On the other hand, near the sphere, shear-dependent kernels play a bigger role.Hence, depending on s s and the duration time of applying this stress, the state of the fluid is expected to evolve locally at different rates determined by the dominant mechanism (kernel in the model context) and ultimately evolve individually.Finally, this offers a velocity profile for the settling sphere complying with this structural dynamics.In Fig. 6, the diagram of settling regimes over all ranges of s s and / 0 / max is disclosed.In general, for a given microstructure, no-motion, continuously decreasing velocity, steady settling, and rapid fall regimes occur in this particular order by increasing s s .Considering that by aging the degree of structuring increases, such a diagram is comparable with the experimental results. 9,59o correlate settling behavior with rheology, we numerically solved a stress-controlled flow curve test.Then the consequent flow curve is compared to the plot of s s -_ over the surface A in Fig. 1 for each settling case.As a result, we obtain a "settling flow curve," which probes the materials' deformation rate under different stress loading.It is worth mentioning, in the CFD simulation, the settling time for a light sphere in a fluid with .Thus, to understand the dynamics of settling based on a rheological test, the difference in the probing times naturally occurring in settling needs to be addressed.This can be done by varying the number of probing points or changing the measuring time of the points.In Fig. 7(a), where / 0 / max ¼ 0:7, and in Fig. 7(b), where / 0 / max ¼ 0:99, an increasing range of 5 < s < 10 4 Pa is applied for one second to 3000 and 20 consecutive probing points, respectively.Appropriately considering the settling time of CFD simulation, the stress-controlled flow cure that is measured using a rheometer type of approach precisely fits the flow curve obtained from the settling.Accordingly, the no-motion regime is comparable to the yield stress transition point of the flow curve where the shear rate goes to zero.Then, by increasing s, the settling velocity increases, and two regimes are observed over the non-monotonic part of the flow curve (continuously decreasing velocity and steady settling).Finally, the rapid fall regime appears at high-stress values, where the flow turns to Newtonian and the viscosity plateaus.Calculating the thixotropic number N thix [Figs.7(c) and 7(d)] also indicates that by increasing the shear rate, where N thix ! 1, rapid fall regime is the most likely settling velocity profile.Beyond that, by decreasing N thix , the continuously decreasing velocity regime is often observed.Here, _ c approaches lower values; therefore, the rate of Brownian kernel K B increases as expected [Figs.7(c) and 7(d)].As a direct result of this, and particularly at the lower degree of structuring, the settling sphere faces a more structured layer of fluid in front of itself.Consequently, its velocity continuously decreases over time. 9,29At a high degree of structuring, a relatively small rate of Brownian kernel is compensated by having a longer settling time and thus the same phenomena lead to the continuously decreasing velocity regime.
In addition, we numerically solved a rheological creep test.Figures 8(a) and 8(b) show how the viscosity relaxes by imposing different stress.The final state of viscosity and associated relaxation times are independent of the / 0 .However, in the beginning, at higher / 0 , a halting time is required for the fluid to deform adequately.This halting time is 1 s when s ¼ 10 000 Pa and increases to 10 4 s when s ¼ 5 Pa.In contrast, at lower / 0 , deformation is undertaken almost immediately after the stress is imposed, and no halting time is obtained.In comparison with the settling simulation, which is essentially a creeping flow-FIG.7. The imposed stress vs average shear rate near the sphere (surface A in Fig. 1).Different settling regimes show the settling flow curve that is compared to a stress-controlled flow curve obtained from the numerical solution of Eq. ( 5) (gray curve).In (a), / 0 / max ¼ 0:  like scenario, observing longer settling times at a higher degree of structuring implies a local halting time.Thus, due to the fact that the sphere has more time to interact with the local microstructure, 9 by approaching / max , this halting time is expected to add more irregular behavior in the velocity profile of the settling sphere [Fig.3(a)].Such a halting time also was predicted in the experiment of dropping steel spheres in aged samples of TEMPO-CNF gels. 9In addition to the halting time, a viscosity bifurcation 13,22 is obtained for Judging by Figs. 7 and 8, settling regimes are a direct consequence of how the microstructure evolves over time, due to the stress caused by the gravity acting on the sphere.As mentioned earlier, this dynamic evolution controls the rates of shear-dependent and shearindependent kernels of thixotropy.In Fig. 9, the magnitudes of the kernels describing the thixotropy are compared together.As expected, the Brownian kernel is dominant in pursuing a continuously decreasing velocity regime, while in rapid fall regime, the shear-dependent fragmentation and aggregation kernels are dominant, particularly by increasing / 0 / max .Since distinguishing the steady settling from the continuously decreasing velocity regime in this regard is difficult, we placed them in the same panel [Fig.9(b)].Still, in a steady settling regime, we expect a smaller rate for the Brownian kernel.It is worth mentioning that the set of parameters in Table II, in general, suggests a higher value for K B .Therefore, to study this, we carried out simulations at varying n b .As Fig. 9(c) indicates, by reducing n b , the steady settling regime becomes a continuously decreasing velocity regime.This results from the increasing impact of the Brownian kernel on the structure dynamics.

C. Fore-aft asymmetry
We observed the onset of fore-aft asymmetry in the flow field as a result of thixotropy.Here, there is no negative wake or velocity overshoot, which makes it different from fore-aft asymmetry in a viscoelastic fluid. 28,51In Fig. 10, two fore-aft asymmetry and symmetry cases are compared.Figures 10(d q and streamline of U, respectively.As Fig. 10(a) shows, in this case, a low viscosity tail is drawn behind the sphere, while ahead of the sphere on its surface at 30 < h < 90 , viscosity reaches its largest value in the domain (K B follows similar behavior with opposite values).In addition, as Fig. 10(b) shows, the value of _ c peaks in front and right next to the sphere (the variation of K sha and K shf follow that of the shear rate  II.field).Therefore, as Fig. 10(c) shows, the onset of asymmetry in the ratio between breaking down and up of the microstructure, indicated by the N thix , leads to flow field asymmetry around the sphere.In contrast, when the flow field is fore-aft symmetric [Figs.10(i) and 10(j)], variations of the gð/Þ and _ c ahead and behind the sphere exhibit the same symmetry [Figs.10(f) and 10(g)].Therefore, thixotropic kernels and consequently N thix are also symmetric.In conclusion, the onset of asymmetry in an inelastic thixotropic fluid is mainly a viscosity-originated phenomenon, where the viscosity gradient causes a faster reduction of the flow field ahead of the sphere than its growth rate behind it.Furthermore, as shown in Fig. 11, variations of gð/Þ; _ c, u z , and s zz behind and ahead of the sphere are examined for a sphere of the same size, but different densities dropped into a given fluid.Similar to Fig. 10(a), gð/Þ under the sphere on its surface is maximized and minimized behind it.Also, _ c on the surface ahead of the sphere is larger than behind it.The normal stress shows no overshoot s zz and varies in a similar manner, relaxing far away from the sphere.While the u z in front of the sphere quickly reaches zero, behind the sphere, it slowly decays causing the symmetry breaking.This effect becomes more pronounced with both increasing sphere density over time.
Comparing Figs. 6 and 12(a) reveals that the onset of asymmetry occurs for all settling regimes.Additionally, it is most likely to occur at a reduced degree of structuring and elevated gravity-induced stresses [Fig.13(a)].In such circumstances, asymmetry is developed at a high enough averaged value of thixotropic number N thix .
To further study the relaxation of the viscous stress in the rheological model, we numerically solved a rheological shear startup test.Similar to the CFD settling simulation, no stress overshoot is observed.In addition, a quite long relaxation time from 10 2 to 10 6 is obtained (Fig. 13), which is independent of the degree of structuring, and comparable to the relaxation time of the creep test in Figs.8(a

IV. CONCLUSIONS
An inelastic thixotropic model is examined to study the correlation between a solid sphere's settling and the fluid's dynamical evolution.The model proposes that evolving the thixotropic fluid during aging occurs through Brownian interactions, which is defined as a shear-independent function of viscosity.Eventually, depending on the local viscosity, this could lead to the development of spatial heterogeneity over time.To study the consequences of this, we included an initial non-homogeneous microstructure in the CFD simulation framework to mimic the development of (1) a heterogeneous microstructure by aging, (2) floc characteristic of some low-density gels such as cellulose nano-fib (CNF) gels, and (3) a microscopically heterogeneous flow, the typical case in practical applications of fluids.Then, this initial microstructure evolves locally due to structural evolution heterogeneity combined with the impact of gravity-induced stress.The key findings of this research are as follows: • In settling solid spheres in a thixotropic fluid modeled by the current scheme, we observe no-motion, continuously decreasing velocity, steady settling, and rapid fall regimes, in the order of increasing gravity-induced stress, respectively.• Fitting a flow curve test to a "settling flow curve" demonstrates that the no-motion regime is equivalent to the yield stress region of the flow curve.Then, continuously decreasing velocity and steady settling are equivalent to the non-monotonic region in the middle shear stress range of the flow curve.Finally, the rapid fall regime is equivalent to the Newtonian plateau region in the high shear stress range of the flow curve.• of the probing method (dropping sphere or rheometry), competition between the thixotropic kernels controls the deformation of a thixotropic fluid.There is a direct connection between the shape of the flow curve and the settling regime.• Deformation time scale of a thixotropic fluid is dependent on the initial degree of structuring.It decreases when the maximum effective volume fraction is approached.This needs to be considered when comparing the solid sphere drop experiment and a flow curve test.• The continuously decreasing velocity regime is a result of the continuous forming of the microstructure in front of the sphere.
Here, the Brownian kernel is dominant in evolving the microstructure, as predicted based on the experimental findings. 9,29• The rapid fall regime is a consequence of the shear-dependent kernels of thixotropy and is not influenced by the Brownian motion.
• Temporal fluctuations in the settling velocity of the sphere only appear when the motion of the sphere is slow enough to catch the spatial heterogeneity.This occurs, for instance, when dropping relatively light spheres in a fluid with a high degree of structuring.
• Following the interpretation based on experimental observations, 9 the model shows that by increasing the degree of structuring, the sphere momentarily halts to overcome the local yield stress.This increases the settling time and leads to an irregular settling regime.• We confirm the previous observations that thixotropy, isolated from elastic effects, causes an onset of flow field fore-aft asymmetry. 18,51This occurs as a consequence of developing a gradient between a low viscous tail along the symmetry axis behind the sphere and a high viscous region near the surface of the sphere in front of it.Such asymmetry is developed whenever a thixotropic fluid with a relatively low degree of structuring is loaded by a high enough gravity-induced stress.
This study provides a framework for a general understanding of the interaction of solid spheres in inelastic thixotropic fluids.However, it would also be of great interest to complement it by experimentally characterizing the heterogeneity of the fluid under shear and by aging.This specificity could target estimating the heterogeneous length scale of a low-density gel, such as TEMPO-CNF or Laponite gels, in comparison to the size and speed of the settling sphere.Additionally, the settling simulation could also be improved by utilizing the finite volume method (FVM) and a yield stress TEVP rheological model. 48his could be helpful in (1) obtaining a more accurate flow field around the sphere, (2) reducing singularity points associated with solving the viscosity, 80,81 and (3) examining other possible material-tomaterial related mechanisms leading to an irregular regime. 48,59xtending the rheological model for each specific low-density gel and parameter determination procedure could be done by considering the practical datasets of rheological measurements. 6

SUPPLEMENTARY MATERIAL
See the supplementary material A for the procedure of determining the rheological model's parameters and supplementary material B for numerical simulation validation.

Fig. 8 (
b) when / 0 / max ¼ 0:99.In this case, viscosity increases when s ! 15 Pa and decreases to a lower value when s 15 Pa.
) and 10(e) indicate the asymmetry of the flow field for the relative velocity magnitude ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 r þ ðu z À kv s kÞ 2

FIG. 10 .
FIG. 10.Comparison between the onset of fore-aft asymmetry and symmetry: dropping a sphere with d s ¼ 3:5 mm and Dq ¼ 6893 kg=m 3 in a fluid with / 0 / max ¼ 0:55 exhibits fore-aft asymmetry.Here, plot (a) is a contour plot of the viscosity field gð/Þ, (b) is the contour plot of the rate of deformation _ c, (c) is the contour plot of the thixotropic number N thix , (d) is the contour plot of the relative velocity ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 r þ ðu z À kv s kÞ 2 q , where u z and u r , respectively, are axial velocities of the fluid in the axis of symmetry z and in the axis r, and (e) shows the streamlines of the fluid's velocity field U.The same sphere in a fluid with / 0 / max ¼ 0:99 exhibits fore-aft symmetry.Here, plot (f) is the contour plot of viscosity field gð/Þ, (g) is the contour plot of rate of deformation _ c, (h) is the contour plot of the thixotropic number N thix , (i) is the contour plot of the relative velocity ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 r þ ðu z À kv s kÞ 2 q , and (j) shows the streamlines of the fluids velocity field U. Contour plots show variations of each field in the same region close to the sphere when it is located 3 cm above the bottom of the domain.Streamlines are shown in the whole domain.

TABLE I .
Computational domain dimensions, sphere, and fluid parameters.

TABLE II .
Parameters of the rheological model.