A phenomenological theory-based viscosity model for shear thickening fluids

A viscosity model for shear thickening fluids (STFs) based on phenomenological theory is proposed. The model considers three characteristic regions of the typical material properties of STFs: a shear thinning region at low shear rates, followed by a sharp increase in viscosity above the critical shear rate, and subsequently a significant failure region at high shear rates. The typical S-shaped characteristic of the STF viscosity curve is represented using the logistic function, and suitable constraints are applied to satisfy the continuity of the viscosity model. Then, the Levenberg–Marquardt algorithm is introduced to fit the constitutive model parameters based on experimental data. Verification against experimental data shows that the model can predict the viscosity behavior of STF systems composed of different materials with different mass concentrations and temperatures. The proposed viscosity model provides a calculation basis for the engineering applications of STFs (e.g., in increasing impact resistance and reducing vibration).


Introduction
Granular matter systems are complex systems that exist in nature and are encountered in everyday life. Research on the rheology of granular matter systems has long been of high scientific significance and industrial application value (Waitukaitis and Jaeger 2012). A shear thickening fluid (STF), a type of granular matter system, is a concentrated particle suspension formed by dispersing nano-or microscale particles in a polar medium. Under appropriate shear conditions, the viscosity of an STF increases rapidly and significantly (typically by several orders of magnitude) with increasing shear rate. This process is reversible; that is, when the stress is removed, the fluid returns to its original state. Under a high-speed impact, an STF may even change from a liquid state to a solid state (Barnes 1989, Boersma et al 1990.
STFs are widely used in the fields of personal protection, impact resistance, and structural vibration reduction due to their unique nonlinear mechanical response and good energy absorption performance.  studied STF-Kevlar composites and demonstrated that STF significantly enhances the ballistic performance of fabrics. Gong et al (2013) studied the impact resistance of STF-treated Kevlar fibers and found that STF significantly improves the impact resistance of Kevlar fibers. Park et al (2015) studied STF composite fiber fabrics under high-speed impact through experiments and finite-element models and showed that the frictional behavior of the STF particles is the main energy absorption mechanism. Fischer and colleagues (Fischer et al 2006, Fischer et al 2010 have applied STF to sandwich beams and found that when the vibration amplitude exceeds the critical vibration amplitude, the relative motion of the sandwich beam causes the shear thickening of the STF in the beam, thereby improving the structural damping performance. Zhang et al (2008a) designed a single-bar STF damper, studied its dynamic mechanical properties, and showed that such a damper can achieve dual control over stiffness and vibration damping under different frequency dynamic loading conditions. Lin and colleagues (Lin et al 2018, Wei et al 2019, Lin et al 2020 have designed a full-scale STF damper and experimentally examined its energy dissipation and dynamic properties; their damper exhibits both viscous and frictional behavior, which substantially benefits vibration control under varying loads.
With the development of experimental technologies and theoretical insights, research on STF has made significant progress in terms of its performance improvement, theoretical simulation, microstructure analysis, and mechanism exploration (Barnes 1989, Hoffman 1998, Yang Li et al 2001, Xu et al 2010, Katz 2014, Jiang et al 2015, Peters et al 2016. In practical applications, STFs are often used in complex environmental conditions, and their shear thickening properties need to be simulated and predicted to provide a basis for material preparation and performance optimization. However, research on viscosity models for STFs is in its infancy (Wei et al 2018, Manica et al 2004.
In this paper, a phenomenological constitutive model function based on an STF is proposed to simulate the viscosity-shear rate relationship curves of STF systems composed of different materials with varying mass concentrations and temperatures. The Levenberg-Marquardt algorithm is introduced to process specific experimental data in order to verify the fitting capability of the function and thereby provide a basis for the material preparation and performance optimization of STFs for engineering applications.

Overview of available viscosity models
The most commonly used constitutive viscosity model of shear thickening materials is the power-law model (Manica and De Bortoli 2004). This model broadly defines the nonlinear growth trend of fluid viscosity as an exponential power growth. In the power-law model, the viscosity model is expressed as follows: where h is the viscosity,  g is the shear rate, K is a material-dependent constant, and n is the power-law index; n is greater than 1, indicating that viscosity increases with increasing shear rate. This function can describe the trend of the instantaneous and exponentially increasing shear thickening effect of STFs. However, beyond the maximum shear rate, the STF shear thinning behavior cannot be fitted. In addition, the shear thinning behavior of STFs is difficult to model at low shear rates.
Addressing this problem, Galindo-Rosales and colleagues  have proposed a viscosity model that is more suitable for describing the rheological behavior of STFs: where h i is the viscosity function that fits the zone i of the general viscosity curve (for i=I, II, III),  g is the shear rate,  g c and  g max are the critical value and maximum value of shear rate. The viscosity model presented in equation (2) represents the three characteristic regions of a typical STF viscosity curve. The modeling involves three expressions that apply to different ranges of shear rates. However, this model has 11 independent parameters, and before the iterative fitting procedure, 6 parameters must be obtained from the experimental data to fit the STF viscosity curve. These aspects of the model make it difficult to find fitting parameters.

Proposed viscosity model
The typical viscosity curve of an STF is shown in figure 1. As shown in the figure, an STF exhibits three rheological characteristic regions under different steady-state shear rates: a shear thinning region at low shear, a shear thickening region at medium shear, and a failure region at high shear. Regions I and III show downward trends for the power function and exponential function, respectively. The shear thickening region shows an S-shaped growth trend. Therefore, based on phenomenological theory, we use a power function (with exponent<0) and an exponential function to fit the viscosity model of shear thinning and failure regions, and use the logistic growth model to regress the viscosity model of the shear thickening region, as shown in equation (3):  Figure 2 shows the graphs of the two logistic growth functions in equations (4) and (5). When the other parameters are held constant, changing parameter a 2 causes the curve to become more 'plump' or 'thin,' whereas parameter b 2 determines the left and right relative positions of the S-shaped curve.
In practical applications of STFs, the thickening ratio (TR) and thickening period (TP) are often used to characterize the rheological properties of shear thickening (Gürgen et al 2016). These measures are defined as follows: where h max is the maximum viscosity of the suspension beyond the thickening point (peak viscosity) and h c is the viscosity of the suspension at the critical shear rate. As evident from equations (6) and (7), the values of TR and TP are related to the critical shear rate, peak shear rate, critical viscosity, and peak viscosity. To improve the fitting capability of the phenomenological model such that it is more suitable for engineering applications, we add the critical viscosity (minimum viscosity) point and the peak viscosity (maximum viscosity) point (i.e., both ends of the shear thickening region) as constraints to the phenomenological model: By substituting equations (8) and (9) into equation (3), the phenomenological model of viscosity is obtained by simplification:

Validation of the proposed viscosity model
In this section, we prove the fitting capability of the proposed viscosity model by fitting experimental viscosity curves of STF systems composed of different materials with different mass concentrations and temperatures. The values of  g , c  g , max h , c and h max must be obtained from experimental data prior to starting the iterative fitting processes. The algorithm used in all of the fitting processes is the Levenberg-Marquardt algorithm. This algorithm combines the Gauss-Newton method and the steepest descent method and seeks the optimal parameter values by iterative chi-square minimization. When the difference between the values obtained in two consecutive iterations is less than the specified tolerance (we set this as 10 −10 during the fitting process), the fitting is assumed to have converged. In addition, the corresponding branch of equation (10) is used to fit each characteristic region of the curve. To simplify the discussion, the correlation coefficient R 2 is 0.9 in all cases in the following comparisons of fitting (in most cases, R 2 0.95).

Fitting capability for different material compositions
A suspension of cornstarch in water is a classic example of shear thickening behavior. The initial viscosity of such suspensions is generally large, and the critical shear rate is small. Figure 3 shows the ability of equation (10) to fit the experimental data set corresponding to a 55 wt% cornstarch suspension in water (Bischoff White et al 2010). In this case, because the curve does not show the failure region at high shear rates, only Region I (the shear thinning region at low shear rates) and Region II (the shear thickening region) of equation (10) are considered for the fitting procedure. As proven by the results in figure 3, the proposed viscosity model offers a satisfactory fitting.
Suspensions of relatively monodisperse polymethyl methacrylate (PMMA) particles in glycerine-water mixtures also exhibit shear thickening behavior. In this case (figure 4), the curve does not show the shear thinning region. Hence, the fitting procedure is performed by considering only Regions II and III of equation (10). The fitting parameters are shown in tables 1 and 2. As illustrated in figure 4, the proposed viscosity model provides a satisfactory fitting with the experimental data (Xu et al 2010) in this case as well.
Shear thickening behavior can also be realized (and the shear thickening performance improved) by adding certain dispersed phase particles to an STF and applying an external field; an example of such a case is magnetorheological STFs (MRSTFs). The left panel of figure 5 presents the fitting of viscosity model with the test data by Zhang et al (2008a and2008b). The right panel of figure 5 shows the viscosity curve of an MRSTF prepared from an STF (suspension of fumed silica in ethylene glycol) and carbonyl iron particles under different magnetic fields (0-100 kA m −1 ) (Zhang et al 2008b). The fitting parameters are listed in tables 3 and 4. Similar to the previous cases, the fitted curve agrees well with the experimental data (Zhang et al 2008a and2008b). Note that the overall viscosity curve shifts upward after the external magnetic field is applied, indicating a further increase in the viscosity. Figure 6 presents the viscosity curves for different solutions of a w , -Mg carboxylato-polyisoprene (Mn=21,000) in decahydronaphthalene with three polymer concentrations at 25°C. The fitting data in figure 6 were from tests by Maus et al (1995). As the shear thinning region is absent in figure 6, only Region II (the   shear thickening region) and Region III (the failure region at high shear rates) in equation (10) are considered for fitting. The figure shows that the proposed viscosity model has good fitting capabilities under different mass concentrations. The fitting parameters and experimental data (tables 5 and 6) change monotonically with  (2008a and 2008b)). Viscosity as a function of shear rate and magnetic field for an MRSTF composed of a 30 wt% suspension of fumed silica in ethylene glycol and carbonyl iron particles. The weight ratio of the iron particles and STF is 5:100 (rightpanel test data from Zhang et al (2008b)).  Figure 6. Phenomenological function simulating the viscosity of an STF (solutions of a w , -Mg carboxylato-polyisoprene (Mn=21,000 g mol −1 ) in decahydronaphthalene at 25°C and concentrations of 2.51, 2.54, and 2.56 g ml −1 , with test data from Maus et al (1995)).  In equation (19), the viscosity model contains only one parameter, mass concentration w, which means that it can be used to predict the rheological behavior of STFs with different mass concentrations.

Fitting capability for different temperatures
At present, silica-ethylene glycol (or silica-polyethylene glycol) suspension is the most common and widely used STF. The initial viscosity of the system is low, and the critical shear rate for shear thickening behavior is moderate. The experimental data (Tian et al 2013) and the model fitting for silica-ethylene glycol suspensions at 20°C, 40°C, 60°C, and 80°C are compared in figure 7. As evident, the model simulation perfectly fits the experimental data at all temperatures. The suitable parameters (for i=I, II, III) are summarized in table 7. In table 7, most parameters show a monotonically increasing or decreasing trend with increasing temperature; accordingly, the following equations are proposed to capture these trends between the five parameters and temperature T.    In equation (29), the proposed viscosity model contains only one parameter, temperature T, which means that it can be used to predict the rheological behavior of STFs at different temperatures. The structure diagram of an STF damper is shown in figure 8. The diameters of the piston and the piston rod are represented by D p and D r The lengths of the piston and the cylinder are expressed by L p and L. The gap between the cylinder and the piston is denoted by h. The effective fluid annular gap includes the entire annular space between the outer diameter of the piston and the interior of the damper cylinder. When the piston moves forward or backward, the STF is pushed through the annular gap. When the force acts on the damper, the piston moves in the cylinder at a speed of V p , and a pressure difference ΔP develops behind and ahead of the piston. The effective cross-sectional area of the piston is represented by A p . Accordingly, the damping force of the STF damper can be expressed as F d =ΔP×A p . Because h represents the height of the gap, the inner diameter of the damper cylinder housing can be written as (D p +2h). The flow of the STF in the gap between the cylinder and the piston meets the following assumptions: (1) The STF inside the damper is a sufficient, continuous, and incompressible fluid; (2) STF is a homogeneous fluid in the flow area of length L p and gap h; (3) The STF in the damper exhibits laminar flow; and (4) The flow gap h is much smaller than the piston diameter D p , and the flow of the STF in the gap can be considered as the flow between the parallel plates. Figure 9 shows the different operating modes of the STF damper. According to the stress state and flow form of the STF within the gap height of h, the working mode of the STF damper can be divided into three modes: flow mode, shear mode, and squeeze mode. The flow mode is shown in figure 9(a): the upper and lower plates are fixed, and the STF flows through the gap between the two plates under the effect of the pressure difference ΔP. If the effective cross-sectional area of the fluid is A L , the damping force of the STF damper can be expressed as F d =ΔP×A L . The shear mode is shown in figure 9(b): the upper and lower plates move relatively, and the STF in the cylinder is affected by the shearing force. The different shear strain rates at different heights cause the viscosity of the STF to change, which in turn generates damping forces. The squeeze mode is shown in figure 9(c): the upper and lower plates move toward each other, and the STF is squeezed to flow under the action of pressure, which in turn is subject to shearing action and generates damping force.
Based on the foregoing assumptions and analysis, the stress and flow velocity distribution of the STF can be analyzed by a combined model of shear mode and flow mode, as shown in figure 10. The width of the compound model can be expressed in terms of the average perimeter b of the annular flow path of the damper, where The control equation for the STF flow model with parallel plates with a gap height in the range of h can be expressed as follows: Consider that the volume flow Q has the following form: where A P is the effective area of the piston. Substituting equations (36)-(38) into equation (39), the pressure gradient can be expressed as follows: The damping force arising from the difference in piston pressure can be expressed as follows: By substituting equation (40) into equation (41), the damping force of the STF damper can be expressed as The experimental data from  depicting the shear thickening performance of a 20 wt% STF at different temperatures is shown in figure 11. Because the shear thickening region of the STF is affected by various factors, such as the dispersed phase particles, dispersion medium, additives, applied field, and temperature, this region is relatively uncontrollable. The structural parameters of the STF damper are usually optimized based on the experimentally measured STF shear thickening region and the relationship between the shear strain rate and the working speed of the STF damper such that the shear strain corresponds to the working speed of the STF damper. The shear strain rate lies in the shear thickening interval of the STF. Therefore, in the fitting of the STF viscosity formula, only the viscosity performance in the shear thickening interval is considered. The fitting parameters are given in table 8.The parameter data and test data of the shear thickening section in figure 11 corresponding to different temperatures are listed in table 8. Using the Levenberg-Marquardt algorithm, the following equations (43) and (44) are obtained by fitting to show the relationship between the parameter and the temperature T. Similarly, the following formula is obtained by fitting to obtain the relationship between the experimental data and the temperature T.   +  -+  -+  +  +  -+  +   -+  +  --+  +  -+  +  --+  +   ---------T  T  T  e  T  T  e  T   e  T  e  T  T  e  T  T  e  T  T   13  127  13  1  In summary, the STF viscosity η at different temperatures can be determined by equation (49). The working parameter of the damper is the moving speed of the piston in the damper, not the shear rate of the STF. Therefore, the relationship between the shear rate and the moving speed of the piston needs to be determined. The moving speed of the damper piston is V p , and the flow speed of the fluid is V l . According to the continuity principle of fluid mechanics, within any given period of time, the flow through the gap between the piston and the cylinder equals the flow from one side of the piston to the other. Therefore, the following equation can be obtained: where A p is the effective area of the piston, and A l is the area between the piston and the cylinder. Then, the velocity of the fluid can be obtained as follows: The shear rate can be determined from the piston speed, the diameters of the piston and piston rod, and the height of the flow gap. Substituting equation (53) into equation (49), the viscosity in equation (42) can be determined. Finally, the relationship between the damping force of the STF damper considering the effect of the temperature and the working speed of the STF damper is as follows:

Conclusions
By applying phenomenological theory, we propose a continuous viscosity model for STFs through the characteristics of the S-curve shown by the STF viscosity curve. The proposed function not only represents three characteristic regions in the typical viscosity-shear rate curve associated with STFs, but also predicts the viscosity of STFs with different material compositions at different temperatures and mass concentrations. In addition, through fitting of different experimental data, we find that the proposed model can satisfy the specific case wherein the failure region does not exist.