UNSTEADY FLOW OF THIXOTROPIC COLLAGEN SUBSTANCE IN PIPES

Unsteady flow of thixotropic liquid in pipes is solved by 1D and 2D numerical methods using the same constitutive equation — the only difference is in the radial diffusion of the structural parameter. Comparison shows that the neglected diffusion of structural parameter implicates a much stronger effect of thixotropy. The models are applied for analysis of the observed hysteresis of hydraulic characteristic of collagen.


Introduction
Thixotropy is a property of fluids having viscous characteristics that depend on their deformation history, see review articles by Mewis [1] and Barnes [2].Approximate solutions of steady flow of thixotropic liquids in a pipe based upon a constant power-law radial velocity profiles were presented by Sestak and Zitny [3] and by Kemblowski and Petera [4].A steady flow and variable radial velocity profile were considered in the FEM solution, Sestak et al. [5], using Houska's constitutive model [6] that takes into account yield stresses.Solutions describing unsteady flows of thixotropic fluids are not so frequent, and are concerned mostly with the specific problem of restarting waxy crude oil in pipelines, starting from the crude approximation by Zitny [7] or simulations presented by de Oliveira et al. [8] up to the most complex CFD methods published in a series of papers by Wachs et al. [9] and Negrao et al. [10].These CFD solutions of compressible fluids described by Houska's model of thixotropy are based upon techniques developed by Vinay et al. [11,12], for compressible Bingham fluids.The methods can be characterized as follows: 2D orthogonal staggered grid, augmented Lagrangian multiplier method for pressures and Uzawa algorithm applied to the resulting saddle point problem.All the previous models assume purely convective transport of structural parameter λ.The exception is a generalization of the Moore's model of thixotropy with a molecular diffusion term suggested by Billingham and Fergusson [13].The 1D and 2D models presented in this paper introduce a generalization of Houska's model by adding a diffusion term into the transport equation for λ.The 1D model assumes infinite radial diffusion and uniform radial profile of λ, while the 2D model assumes zero diffusion and a strong radial variation of λ.The aim of the following analysis is to compare predictions of pressure drops corresponding to the 1D and 2D extremes.

Motivation
The problem of unsteady flow of thixotropic liquid in a pipe was initiated by the industry and by the demand for an on-line recording of rheological properties of collagenous material used for the extrusion of collagen casings.The time independent Herschel Bulkley rheological model was capable of describing rheograms of collagenous material for a wide range of deformation rates; however it was not able to describe a slight observed hysteresis of hydraulic characteristic (pressure drop versus flowrate) at transient flows.The thixotropic properties of collagen were identified as a possible reason.The suggested thixotropic models can be used not only for modeling of the hysteresis, but also for simulation of production lines (collagenous materials used for extrusion of vascular grafts in biomedicine (Kumar [14]), or as sausage casings in food industry, Deiber et al. [15]).Nevertheless, information on rheological properties of collagen is restricted due to problems with the application of rotational rheometers.Usually only rheograms measured by capillary rheometers are available, giving no information about thixotropy.The only exception is a paper by Deiber et al. [15], reporting thixotropy of collagen suspensions at low concentration (0.5-3.7 %), using a rotational rheometer and a cone and plate configuration.

Methods
This paper is devoted to numerical modelling of unsteady flow of thixotropic liquids in circular pipes.As a possible application of the suggested methods the analysis of hysteresis of the hydraulic characteristics observed at a production line of collagen casings will be presented.

Experimental setup and procedure
The experiments were realized at the production extrusion line DEVRO Jilemnice shown schematically in Fig. 1.The processed material (collagenous matter: approximately 7 % mass fraction of bovine collagen and water) is delivered from a storage tank by an AQM 57-20/S003 positive displacement pump, through a relatively long pipe (L = 4.287 m, R = 0.0106 m) to an extruder.The flow rate is controlled by varying the pumping speed and the pressure drops (changes) are recorded by DMP 331P (BD Sensors) pressure transducers.Special experiments characterized by a gradual increase and decrease of mass flow rate at constant temperature and constant composition of tested matter were carried out during a break of production in the extrusion processing line with the goal of identifying parameters for the Herschel-Bulkley rheological model.Most details concerning the technical realization and the composition of the processed collagenous material are confidential, and only limited data files (time; flow rate; pressure drop) were released.Nevertheless these experiments serve as illustration and as a test of the possibilities of the suggested methods in real applications.
The processed material, at relatively high concentrations of bovine collagen (7 %), is a viscoelastic paste that looks like "silly putty".Our preliminary experiments in the laboratory confirmed linear viscoelasticity at small deformations (measured by the cone and plate oscillating rheometer Haake giving a module of G of 5 kPa, G of 2 kPa).Experiments with rotational rheometers at constant deformation rates failed due to very high consistency.Therefore, no independent evaluation of thixotropic time constants was realized.Laboratory tests of several collagen samples, carried out by using size exclusion chromatography and UV detection, identified the existence of very long collagen fibres (19 % of the longest fraction 780 kDa).

Constitutive equation
The Houska's thixotropic model [6] is a generalization of the routinely used Herschel Bulkley constitutive equation, simplified for the special case of simple shear flow to the scalar equation for shear stress τ : ( The structural parameter λ = 1 describes a fully recovered internal structure, while λ = 0 corresponds to a completely destroyed structure.Time changes for structural parameter λ are described by the transport equation: The diffusion term on the right side is usually neglected and only the terms for regeneration (a) and structure destruction (b) are considered.

Radially independent structural parameter (1D model)
The assumption that the structural parameter λ depends only on the axial coordinate and on time λ(t, x) would represent a great simplification.This case is probably more relevant for polymeric chains with clusters of a size comparable with the pipe radius (the tested collagen is characterized by extremely long chains with molecular mass 780 kDa).

Flow of a Herschel-Bulkley fluid in a pipe
If λ is independent from the radial coordinate, then the effective consistency K and the yield stress τ y depend only upon time and the axial coordinate: (3) Therefore, it is possible to use the RMW (Rabinowitsch, Mooney, Weissenberg) equation for volumetric flow rate expressed as a function of wall shear stress (or gradient of pressure): V Because we need to calculate the wall shear stress for a given flow rate, it is necessary to iterate the inverse relationship: After just a few iterations, the process converges for arbitrary model parameters.Equations ( 3)-( 6) were presented recently by Zitny et al. [16] in a simplified 1D model.

Flow of a Herschel-Bulkley fluid in a pipe
The transport equation ( 2) can be integrated across the cross section of the pipe and the last term describing the averaged decomposition of the structure can be expressed in an analytical form, because the radial velocity profile is known (the mentioned model in Zitny et al. [16] assumed a constant shear rate corresponding to Newtonian velocity profile): For D λ = 0, Equation ( 2) reduces to a hyperbolic equation that can be integrated analytically along its characteristic dx = u dt giving: where λ(t) is the value of the structural parameter of a fluid particle having a value of λ 0 at time t 0 , assuming that the particle is under the influence of a constant shear rate (constant flow rate) within the time interval (t 0 , t).In cases with variable flow rates, it is necessary to apply (8) using shorter time steps during the integration.The 1D numerical solution by method of characteristics calculates nodal values of structural parameters at a new time t (k+1) from the values of the structural parameter λ 2 , . . .at equidistant nodes (x 1 = 0, x 2 = ∆x, x 3 = 2∆x, . ..) and previous time t (k) .The time step ∆t (k) = t (k+1) − t (k) is not a constant and is determined by the volumetric flow rate so that the fluid particle moves a distance of ∆x: (Courant Friedrichs Lewy criterion CFL u∆t/∆x = 1).
The new value of λ (k+1) i is calculated by integration of a fluid particle along the characteristic from the old value λ (k) i−1 using (8).After the new values for structural parameters are updated, wall shear stress, pressure gradient and overall pressure profile can be calculated from (6).

2D model of velocities and structural parameter profiles (finite differences)
Results that are probably more realistic can be obtained by taking into account non-uniform structural parameter profiles.Values of axial u ij and radial v ij velocities and structural parameter λ ij are calculated at grid points x i , r j of orthogonal nonuniform collocated mesh using the finite differences method.Radial profiles of axial velocities are evaluated from axial momentum balance (assuming linear radial shear stress profile and a constant value of structural parameter in the annular sections r j , r j+1 for j = 1, 2, . . .n R−1 ).Radial velocities follow from the continuity equation and the profiles of structural parameters are calculated from discretized transport equation ( 2) assuming purely convective transport (therefore D λ = 0).More details of this numerical solution are presented in the Appendix.

Results and discussion
The numerical methods described in Section 3 were implemented in several Matlab programs and tested with the goal of evaluating the differences between the predictions made using the 1D and 2D model.Operational parameters were selected according to the geometry of the test setup (pipe L = 4.2 m, R = 0.01 m).The flow rate was progressively increased from 5 × 10 −7 to 2 × 10 −5 m 3 s −1 and then progressively decreased as shown in Fig. 2.
A comparison of the predicted hydraulic characteristic using the 1D and 2D model for selected rheological parameters is shown in Fig. 3.The 1D model, with a constant radial profile of structural parameter, predicts pressure drops approximately 10 to 20 % higher than the 2D model with a varying radial profile of structure.
As expected, the 2D model predicts a more intensive breakup of the thixotropic structure; particularly in the vicinity of the wall (the uniform distribution of the structure decay obviously underestimates the significant effect of the wall layer in the 1D model).The simulations were carried out for different meshes (spatial discretization using 21, 51, 101, and 451 nodes in the axial direction, and 11, 21, and 31 nodes in radial direction when testing the 2D method).The times steps were fully determined by spatial steps in the method of characteristics (1D), while in the 2D method, time steps were not restricted.Nevertheless, in the simulations the time steps were selected to be identical with the 1D method.Tests confirmed that both the 1D and 2D method were stable and consistent in all tested cases.Fig. 4 presents the numerical prediction of the 2D model and the data obtained experimentally with the collagen material used in the experimental setup described in paragraph 3.1.The volumetric flow rate was gradually increased from 7 × 10 −7 m 3 s −1 to 5.1 × 10 −5 m 3 s −1 .The rate of volumetric change was much slower than in previous tests (more than 20 minutes up and 20 minutes down).Such a slow variation of flow rate increases the number of time steps (4100 time steps for 151 grid points in the axial direction) and increases the run-time necessary for preliminary identification of thixotropic model parameters.The parameters (a, b, n, m, K, ∆K, τ y , ∆τ y ) were identified (approximately) by linear search using the least square criterion (the sum of the squares of the deviations between measured and predicted pressures during up and down phases).It is obvious that the models of thixotropy (regardless of whether it was the 1D or 2D model) were not able to describe the hysteresis at high flow rates, a trend that was also observed in previous simulations (Fig. 3).The fact that the hydraulic hysteresis of thixotropic liquid is significant only at low flow rates is closely related to the residence time t RT D = L/ū of fluid particles.The relative change of structural parameter (and hysteresis) can be quantified as For the analyzed case (and for the model parameters shown in Fig. 4) the ∆λ/λ ratio ≈ 100 at the minimum flow rate, and only ≈ 0.01 at the maximum flow rate.

Conclusions
The 1D method of characteristic and the 2D finite differences method were designed to simulate transient flow of thixotropic incompressible liquids in pipes using Houska's model with variable yield stress.The 1D model assumes that the radial diffusion mitigates non-uniformities and large gradients of structural parameter at the wall.The 2D model assumes only convective transport of structure and predicts large changes of structure at the wall where local residence times go up to infinity.The two suggested methods represent, in fact, the lower and upper bound of the pressure drop (the 1D method underpredicts and the 2D method overpredicts the effect of thixotropic structure decay).The developed models were used to evaluate rheological properties of bovine collagenous matter.We observed a phenomenon resembling hydraulic hysteresis in a production line facility and this paper is an attempt to explain this behaviour by thixotropy.However, the hysteresis can only be explained by thixotropy at low flow rates and the hysteresis observed at high flow rates would need a different explanation.

Appendix
Axial velocities u ij , radial velocities v ij and structural parameter λ ij are calculated in each time step in nodes of 2D rectangular mesh (i=index in the axial, j in the radial direction).The radial profiles of axial velocities are evaluated from the axial momentum balance (assuming a linear radial shear stress profile and a constant value of structural parameter in the annular sections r j , r j+1 for j = 1, 2, . . .n R−1 ), giving The nodal velocities were evaluated recursively from (12) starting from zero velocity at the wall.The corresponding flow rates in annular sections between the wall and radius r j are The velocities and flow rates in ( 12), ( 13) are expressed as analytical functions of wall shear stress τ w and the wall shear stress is determined by the prescribed flowrate.This constraint is in fact the algebraic equation V (t) = V (t, x i , r j ) , which has to be solved iteratively at each time step and at all cross sections x i .The radial velocity can be calculated from the mass balance expressed in terms of the radial profiles of flow rates Let us assume that the nodal values of structural parameter λ (k) ij are known at time t (k) and the values λ (k+1) ij at time t (k+1) and node x i , r j we want to calculate.It is first necessary to calculate the trajectory of the particles, or more precisely, to calculate the initial position of the fluid particles (coordinates ξ, η shown in Fig. 4) that ends at node x i , r j .Assuming that the change in the next time step ∆t (k) = t (k+1) − t (k) is small, the trajectory remains inside the cell formed by cross sections x i−1 and x i , and the increments ξ, η are given by velocities evaluated at node i, j from ( 12) and ( 14): The value of the structural parameters of a fluid particle at a "previous" time t (k) at position ξ, η can be interpolated from the nodal values (using bilinear interpolation): i−1,j−1 ξη for ξ > 0, η > 0; ( 16) The values of the structural parameters at a "new" time are obtained by integration of transport equation (2) along the particle trajectory (i.e., along the characteristic): • λ (k) (ξ, η) e −(a+b( γ(k) ij ) m )∆t (k)   .(18) This relationship can be applied to all internal points; the value of structural parameters at the wall can be calculated directly from equilibrium The method described is explicit because the velocities in (15), the associated trajectories and rate of deformation are evaluated at the "previous" time.i Index of cross-section (axial coordinate) j Index of node in radial direction w Subscript referring to the wall

Figure 2 .
Figure 2. Time course of flow rate for tested cases.

Figure 5 .
Figure 5. 2D mesh and interpolation of structural parameter at a "previous" time step.

anx
Coefficient of structure regeneration [s −1 ] b Coefficient of structure breakdown [s m−1 ] D λ Diffusion coefficient [m 2 s −1 ] K Consistency coefficient [Pa s n ] ∆K Increment of consistency coefficient [Pa s n ] L Length of pipe [m] m Breakdown index [-] Flow behaviour index [-] p Pressure [Pa] r Radial coordinate [m] R Radius of outer pipe [m] t Time [s] tRT D Mean residence time [s] u Axial velocity [m s −1 ] ū Mean axial velocity [m s −1 ] v Radial velocity [m s −1 ] V Volumetric flow rate [m 3 s −1 ] Axial coordinate [m] γ Shear rate [s −1 ] λ Structural parameter [-]τwWall shear stress [Pa] τy Yield stress [Pa] ∆τy Increment of yield stress [Pa]