UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Taylor-Couette instability of a Bingham fluid 2003

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2004-0054.pdf
ubc_2004-0054.pdf
ubc_2004-0054.pdf [ 2.12MB ]
Metadata
JSON: 1.0080059.json
JSON-LD: 1.0080059+ld.json
RDF/XML (Pretty): 1.0080059.xml
RDF/JSON: 1.0080059+rdf.json
Turtle: 1.0080059+rdf-turtle.txt
N-Triples: 1.0080059+rdf-ntriples.txt
Citation
1.0080059.ris

Full Text

T A Y L O R - C O U E T T E INSTABILITY OF A B I N G H A M F L U I D by MARIA PAULINA LANDRY BMath (Applied Mathematics, Honours Co-op) University of Waterloo, 2001 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Mathematics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 2003 © Maria Paulina Landry, 2003 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer- ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Date De~o / b , 2JQ0 3 Abstract Stability of a Bingham fluid in Taylor-Couette flow is studied. Of interest is to determine the critical values of the inner and outer cylinder angular velocities which cause the linear instabil- ity in the system. We introduce a perturbation into the basic flow, and examine the linearized Navier-Stokes equations for the perturbation. Assuming a normal mode solution for the per- turbation, we obtain an eigenvalue problem which determines when the perturbation will grow or decay with time. In the case of a plug region at the outer cylinder wall, additional boundary conditions are applied at the yield surface. The marginal stability curve is generated numeri- cally using a finite difference method. In addition to numerical analysis, heuristic arguments and energy methods are applied to find a lower bound for stability. ii Table of Contents A b s t r a c t i i Tab le of Conten ts i i i L i s t of Tables v L i s t of F igures v i Acknowledgement v i i C h a p t e r 1. In t roduc t ion 1 1.1 Historical Perspective of the Taylor-Couette Problem 2 1.2 Bingham Fluids and the Taylor-Couette Problem 3 1.2.1 Applications 4 1.2.2 Previous Work 6 1.3 Thesis Outline 6 C h a p t e r 2. Gove rn ing Equat ions 8 2.1 Basic Solutions 8 2.2 Dimensionless Parameters 10 2.3 The Velocity Profile 11 2.4 Relation Between TJ and 0 14 2.4.1 Region I 14 2.4.2 Region II 15 2.5 Visualizing Regions I, II and III 16 2.6 Limit as B ->• 0 17 C h a p t e r 3. L inea r S tab i l i t y Ana lys i s 20 3.1 Linear Perturbation Equations 20 3.2 The Method of Normal Modes 25 3.3 Boundary Conditions for the Perturbation 27 3.3.1 Region I: No Plug 27 3.3.2 Region II: Plug on the Wall 28 C h a p t e r 4. A x i s y m m e t r i c D is tu rbance 31 4.1 Axisymmetric Equations 31 4.2 Numerical Methods 34 4.2.1 Basic Flow Equations 35 4.2.2 Finite Difference Method 35 4.2.3 Finding the Critical Rei 35 C h a p t e r 5. Resu l t s for an A x i s y m m e t r i c D is tu rbance 39 5.1 Numerical Results • 39 iii Table of Contents 5.1.1 Compaxison with Newtonian Results 39 5.1.2 Results for B > 0 39 5.2 Heuristic Results 43 5.3 A Lower Bound for Stability 47 5.3.1 Simplifying the Integrals 50 5.3.2 Bounding the Integrals 53 5.3.3 Applying the Poincare Inequality 56 5.3.4 Results 57 C h a p t e r 6 . C o n c l u s i o n s 5 9 B i b l i o g r a p h y 6 0 A p p e n d i x A . N u m e r i c a l A p p r o x i m a t i o n s o f t h e E i g e n v a l u e P r o b l e m 6 1 iv List of Tables 5.1 Approximate least squares values for a and q for B > 200 with N = 50 43 v List of Figures 1.1 Taylor-vortex flow of a Newtonian fluid [10] 2 1.2 Regions observed for the Taylor-Couette flow of a Newtonian with R1/R2 = 0.883, An- dereck et al [1] 3 2.1 Regions with rj = 0.883 for various values of B. Region III is the line Re\ = r]Re2- • • 18 2.2 Regions with B — 5 for various values of 77 18 4.1 Re\ is too small. Instability has not yet occurred 36 4.2 -Rei is too large. Instability has already occurred 37 4.3 itei = Relcrit 37 4.4 Marginal stability curve for B = 10, 77 = 0.883 38 4.5 Flowchart of marginal stability code 38 5.1 Comparison of our numerical calculations against data from [1] for m = 0 40 5.2 Marginal stability curves versus regions for 77 = 0.883 with N = 50 40 5.3 Marginal stability curves 77 = 0.883 and varying B 41 5.4 ReiCTit becomes progressively independent of Re2 as B increases 42 5-5 ReXcrit for B = 1000 with the asymptotic solution Relcrit = aB"; a = e 4 0 4 6 7 , q = 1.2169. 42 5.6 For a larger gap, Re\crit is not necessarily an increasing function of B 44 5.7 Marginal stability curves with 77 = 0.5 and varying B 44 5.8 Marginal stability curve for a Newtonian fluid along with the Rayleigh criterion. . . . 46 5.9 Solid lines: Rayleigh criterion. Dashed line: Rigid rotation Re\ = r\Re2\ n = 0.883. 48 5.10 Solid line: f(Rei,Re2) = 0. Dashed line: Rigid rotation Rei = r)Re2; r\ = 0.883 58 vi Acknowledgement This thesis was conducted under the supervision of Dr. Ian Frigaard and Dr. Mark Martinez. I would like to thank them for their guidance, insight and encouragement throughout. I would like to acknowledge the help of Dr. Randall Tagg from the University of Colorado who provided data for this thesis. For funding this research, I would also like to thank the Advanced Systems Institute of British Columbia and Metso Corporation. vii Chapter 1 Introduction We will consider the stability of a Bingham fluid between two infinitely long concentric cylinders with inner and outer radii R\ and R2 respectively. Both the inner and outer cylinders rotate independently at respective angular speeds Q\ and &2- The Taylor-Couette problem is the study of the flow behaviour between the cylinders. For sufficiently small values of the angular speeds, Couette flow is observed; a laminar azimuthal shear flow between these two concentric cylinders. For Newtonian flow, the angular speeds can be increased so that a critical flow state is ob- tained and the particles travel in helical paths, where the vortices are layered one on top of the other as seen in figure 1.1. These helical vortices are known as Taylor vortices and this primary instability is known as Taylor-vortex flow. Taylor-vortex flow is, however, a stable flow and further instabilities can be reached by increasing the angular speeds. Observe that this linear instability is axisymmetric and periodic. Furthermore, note that the Taylor-Couette flow demonstrates the non-uniqueness of the Navier-Stokes equations, since there are at least two flows which mathematically satisfy these equations with the same boundary conditions: the laminar flow and the perturbed flow. We would like to determine the critical values of &i and CI2 for instability of a Bingham fluid. In the sections to follow, we will describe in detail the properties of a Bingham fluid and how the study of Taylor-Couette flow of a Bingham fluid differs from the study of Taylor-Couette flow of a Newtonian fluid. 1 Chapter 1. Introduction A Q i Figure 1.1: Taylor-vortex flow of a Newtonian fluid [10]. 1.1 Historical Perspective of the Taylor-Couette Problem In 1923, G.I. Taylor [11] was the first to predict analytically, and support experimentally, the occurrence of Taylor vortices with a Newtonian fluid. He assumed an axisymmetric disturbance, and a narrow gap between the cylinders. Since then, there have been many papers written on the Taylor-Couette problem. The literature review by Tagg [10] is a plentiful source for references, methods of approach to the Taylor-Couette problem of a Newtonian fluid, and results obtained. An extensive study by Andereck, Liu and Swinney [1] shows regimes observed for various angular speeds as seen in figure 1.2. In this figure, Re\ and Re2 are the inner and outer cylinder Reynolds numbers respectively (see section 2.2), and the radius ratio is fixed at r? = R1/R2 = 0.883. The transition from Couette to Taylor-vortex flow is observed to be the primary instability. Flow regions beyond the primary instability are shown as well. Of interest in this thesis is to produce the linear instability curve in the (i?e2,.Rei)-plane when the fluid is a Bingham fluid. 2 Chapter 1. Introduction • oL I 1— « 1 _ i _ •^mo -ym -awo - t w o o> IOOQ ' ' ^ 2 Figure 1.2: Regions observed for the Taylor-Couette flow of a Newtonian with R1/R2 = 0.883, Andereck et al [1]. 1.2 Bingham Fluids and the Taylor-Couette Problem Unlike a Newtonian fluid, a viscoplastic fluid is a fluid with a yield stress. That is, the yield stress needs to be exceeded in order for the fluid to flow. Wherever the stress does not exceed the yield stress, there is no flow and this region is called the unyielded region. The region where there is flow is called the yielded region. Viscoplastic fluids can be modeled by one of the Bingham, Casson, or Herschel-Bulkley equations which all provide mathematical relations between the stress and strain-rate. The most widely used model, because of its simplicity, is the Bingham fluid model [2], where the stress/strain-rate relation is given by fij = I pP + ^ I jij f>fy 7 7 = 0 T <Ty where fy and jij are the deviatoric stress and rate of strain tensors respectively, 3 Chapter 1. Introduction T = ( ) 2 is the second invariant of the deviatoric stress tensor and 7 — I 2'^''^' ) is the second invariant of the rate of strain tensor. Here pp is the plastic viscosity and fy is the yield stress. In contrast, the stress/strain-rate relation for a Newtonian fluid with viscosity p is given by Note that once the yield stress is exceeded, the stress varies linearly with the strain rate as does that of a Newtonian fluid. In regions where the yield stress is not exceeded, the fluid moves as if it were a solid. The region where the fluid moves in rigid motion is called a plug region. In the case of a plug, additional boundary conditions are required at the yield surface. For the Taylor-Couette flow of a Bingham fluid, we find that there are three basic types of flow: I. No plug flow, II. Plug attached to the wall, and III. Rigid body rotation. Examples of fluids with a yield stress include: pulp slurries, sewage sludge, drilling muds, unhardened cement, toothpaste and paint. 1.2.1 App l i ca t ions P u l p suspensions: Pulp suspensions are known to exhibit a yield stress. To the lowest order of approximation, we can approximate a pulp suspension with a Bingham fluid model. Papermaking fibre suspensions are a mixture of wood fibre, water and clay as well as small quantities of inorganic salts and polymeric additives. This physical makeup leads to complex 4 Chapter 1. Introduction rheological behaviour with the suspensions showing non-linear stress-strain relationships. The major feature of the stress-strain relationship is that of a non-Newtonian fluid. It is widely known that pulp suspensions do not flow until a certain critical stress (or yield stress) is exceeded. Suspension yield stresses have been either measured directly using viscometers or inferred from dynamic and tensile measurements. For pulp suspensions, the yield stress has been typically correlated with the suspension mass concentration using equations of the form where cm is the mass concentration and cs is the gel concentration. The parameters a and b depend upon the test procedure employed, the pulp type, and the degree of treatment. Values of b have been reported to range between 1.25 to 3. With flowing papermaking suspensions, most experimental data relating the transition from a plug-like to a well-mixed state and then turbulent flow has been observed in pipe flows. The onset of drag reduction, observed to occur in the well-mixed state, has been characterized by criteria based upon power dissipation per unit volume. There has been very little work regarding transition in a rotating device, i.e. a pressure screen. The information generated in this thesis will help address this issue. Electrorheological (ER) fluids: ER fluids are suspensions of dielectric particles in a non- conducting liquid. The primary effect of an electric field is to create a yield stress in the fluid. Typically, this can increase up to 5 kilopascals in the order of a millisecond. Papers such as [5] use the Bingham model to model the flow of an ER fluid. There are many technological applications including shock absorbers and electronically controlled clutches. Often in these applications, the ER fluid serves as a hydraulic lubricant where the ER device involves rotational shear. The Taylor-Couette flow is an archetypical flow and can be used to model instabilities in these applications. For safety concerns, it is important to determine when such instability would occur. 5 Chapter 1. Introduction 1.2.2 P r e v i o u s W o r k Linear stability of Taylor-Couette flow of a Bingham fluid has been studied by Graebel [8] who generalized the methods of Chandrasekhar [3] by using a narrow gap approximation to analyze the stability of a Bingham fluid once the normal mode equations were obtained. In the case of a plug on the outer cylinder wall, the gap was considered to be the fluid region between the inner cylinder and the plug. Graebel concluded that the yield stress acts as a stabilizing agent in the flow. Recent papers on linear stability of Poiseuille flow of a Bingham fluid, such as [6] and [7], present the approach which has been applied to Taylor-Couette flow in this thesis. Methods from [6] include linear stability analysis and applications of boundary conditions at the plug, which we follow in chapter 3. Furthermore, [7] provides a method for finding a conservative estimate for linear stability, which we apply in section 5.3. 1.3 Thesis Outline Our objective is to determine when the linear instability of the Taylor-Couette flow of a Bing- ham fluid occurs. To do this, we will use linear stability analysis. In chapter 2, we first nondimensionalize all relevant parameters, and determine all the basic equations for the flow. Of importance are the dimensionless groups n, Re\, Re2 and B, the Bingham number (see section 2.2), which determine all characteristics for the flow. By finding the stress at any point, we can determine where a plug region occurs. From this point, we obtain three basic types of flow: I. No plug flow, II. Plug on the wall, III. Rigid body rotation. Furthermore, a velocity profile is calculated for each of these flow types. In chapter 3 we introduce a perturbation to the basic velocity profile, and study the motion of the perturbation by linearizing the Navier-stokes equations, and assuming a normal mode expansion for the perturbation. Furthermore, no-slip boundary conditions are applied at the cylinder walls and, in the case of type II flow, additional boundary conditions are applied at the yield surface. 6 Chapter 1. Introduction Since the linear instability for a Newtonian fluid has been observed to be axisymmetric, we focus our efforts on an axisymmetric disturbance in chapter 4. The normal mode approach leads to an eigenvalue problem for the growth rate of the perturbation. Depending on the sign of the real part of the eigenvalue, we can determine whether the perturbation grows or decays with time. We then describe the numerical methods used to solve this eigenvalue problem. Numerical results and analytical lower bounds are presented in chapter 5. 7 Chapter 2 Governing Equations 2.1 Basic Solutions Consider the flow between two infinitely long concentric cylinders with inner radius Ri and outer radius Ri- The cylinders rotate independently with inner angular speed £l\ and outer angular speed Cl2- Let r, 0 and z be the radial, azimuthal and axial directions of motion respectively. Assuming a velocity profile of the form U = V(f)eg gives the stress/strain-rate relation tfe = [frp + - r j Ife T>Ty (2.1) 7 = 0 f <fy (2.2) where p,p is the plastic viscosity, f = \T?e\ is the stress, 7 = {jf-el 1S the strain rate, and fy is the yield stress that needs to be exceeded in order for fluid to flow. A l l of the other stress components vanish. If the yield stress is not exceeded at a point, then the strain rate is zero at that point so that there is no relative motion between particles and the fluid moves as if it were a rigid body. For this reason, the region where f < fy is called a plug region. The momentum equations are given by y 2 0 _ dp df l_d_ f.2 Qf. (r 2rfe) (2.3) (2-4) 8 Chapter 2. Governing Equations There axe two possible boundary conditions that can be applied to the inner cylinder. We can apply either an inner stress Looking at the (9-momentum equation (2.4), we see that the former boundary condition is a more natural choice. Solving for T?O in the (9-momentum equation and applying the inner stress boundary condition gives We see that ffe decreases as r increases, therefore if a plug occurs in the flow it will be against the outer cylinder wall. If the stress at the outer cylinder is greater than the yield stress, then there will be no plug in the flow. If the stress at the inner cylinder is less than the yield stress, then the entire region between the cylinders will be a plug and the fluid will move in rigid body rotation. We obtain three basic solutions: Tfe(Ri) = n (2.5) or an inner velocity V(RJ = Mi. (2.6) III. Solid body rotation II. Partial plug I. No plug For each of regions I, II and III we need to apply the boundary condition V{R2) = CI2R2. (2.8) Denote the location of the yield surface, or plug, by f = Ry. In case of a plug, we apply a third boundary condition 9 Chapter 2. Governing Equations V(Ry) = Ci2Ry. (2.9) 2.2 Dimensionless Parameters To nondimensionalize the problem, we will follow [4] by setting TJ = - \ TJ, r — i f , z = ii, t = i£ , P = 4-P, r = 4-f where d = R2 - Ru f = and P 0 = ^PR}^\ By convention, we take fii > 0. This nondimensionalization gives the rescaled Navier-Stokes equation U t + i2ei(U- V)U = - V P + V - r (2.10) where is the inner cyhnder Reynolds number, and set Re2 = (2.12) as the outer cylinder Reynolds number. The stress/strain-rate relation is given by Tre = (l + ̂ jjre <^ T> B (2.13) 7 = 0 <*=>• r <B (2.14) where 7 = \ jrB\ and T = |-rrd|- Here, B = (2.15) p,vR\Vt\ 10 Chapter 2. Governing Equations is the Bingham number which is the ratio of the yield stress to the viscous stress. If B = 0, then the fluid is Newtonian since no stress needs to be exceeded in order for the fluid to flow. Let 77 = R1/R2. Then . R 1 and R2 get rescaled to B-i = ~~~ (2.16) 1 -v R2 = j^— (2-17) respectively. Note that with this rescaling, the gap width R2 — Ri will always be 1. At the inner cylinder, we can apply either the inner stress boundary condition rre{Ri) = n (2.18) or the inner velocity boundary condition V(Ri) = 1. (2.19) If there is no plug, the outer cylinder boundary condition is given by V(R2) = Q^ = - (2.20) RX T] where Cl = CI2/CI1. If we are in region II, then the boundary condition at the yield surface r = Ry is given by V(Ry) = n^. (2.21) 2.3 The Velocity Profile The rescaled stress is given by 11 Chapter 2. Governing Equations TrO = ^ (2-22) and the three basic solutions are given by I. No plug II. Partial plug III. Solid body rotation Izsl > i < t a i < ( i ) ! M < 1 B — If there is a plug (or unyielded region), the yield surface Ry occurs when \TTQ\ = B. Therefore we solve the equation R = B for Ry to obtain (2.23) Note that in the yielded region, r» > 0 TTQ > 0 + jrg > 0 jre > 0. Similarly, T j < 0 jre < 0, therefore, in the yielded region we have sgn(7rg) = sgn(r;). Using this result we obtain Tt0 = ire + 5sgn(7 r f i) = jr9 + Bsgn(Tj) . (2.24) We can rewrite irg as nR2 ire = Tre - Bsgn{ri) = - jBsgn(rj). (2.25) Using the above information and the fact that jrg = 757 — 7- = r-jjjp ( 7 ) , we obtain the differential equation ( 7 ) ( 2 - 2 6 ) which gives 12 Chapter 2. Governing Equations 1 T R2 V = - - - ^ - r51n(r)sgn(ri) + Cr. 2 r (2.27) R e g i o n I: Apply the boundary condition V(R.2) = fljft- which gives C = ^ - + ^ + flln(il2)8gn(7i) so that for Ri < r < R2 F=|r +r^(ii^)+ i J r l"(^ | S 6 n ( T , )' (2.28) R e g i o n II: Applying the boundary condition V(Ry) = ^ 7 ^ , and using the fact that 7 = 0 for Ry <r < R2 gives V=< ( n R '-r + l-TiR\r ^ L - ^ j + Br ln ( ^ j sgnfa) for Rx < r < Ry for Ry<r< R2 9. (2.29) R e g i o n I I I : We have that 7 = 0 everywhere, so for R\ <r < R2 we obtain V Ri (2.30) 13 Chapter 2. Governing Equations We can express the velocity profile more concisely as f SLr + lTiRlr( ^ - V)+5r ln(—W(T*) f o r R i < r < R o v = ) Ri 2 \R* r V V r ) ( 2 - 3 1 ) 1 Tr-r for RQ<r<R2 V i l l where R0 = mm{Ry, R2}. 2.4 Relation Between Tj and Q. Applying the boundary condition V(Ri) — 1 to the velocity profile in the yielded region gives a relation between TJ and fi. 2.4.1 R e g i o n I In region I we obtain the relation 1 = n + ^TiRi ( f | - l ) + B i * i l n 0 | ) s g n ( T i ) = n + ^TMn2-l) + BRx\n(^jsgn(Ti) (2.32) which can be rearranged as 1 _ fi = ^ (Htf _ 1 ) + ] n ( J _ ^ s g n ( T i ) . ( 2.33) Since we are in region I, we have ^ > therefore ^ - ' ) + " ( ? ) < ̂ +>-(?) = ta(?) - ? +1 <0 <2-34> for n < 1. Therefore, if sgn(Tj) = —1 then 1 — fi > 0 or T j < 0 fi < 1. Similarly, T j > 0 4$ fi > 1. Therefore sgn(Tj) = sgn(fi — 1), and can solve for TJ to obtain 14 Chapter 2. Governing Equations n = T?2 — 1 2.4.2 Region II (sgn(fi - 1)B ln(r?) + (1 - fi) ( ^ p ) ) ' ( 2" 3 5 ) It is not possible to solve for TJ explicitly when we are working in region II. Applying V(Ri) = 1 gives the relation (2.36) or n - 1 = -^r (ln (51)" l4 + 0sgn(Ti)- (2-37) If we are in region II, then 1 < ^ < ^ , therefore ln ( ^ p ) - ^ + 1 < 0. Again we obtain Ti < 0 fi < 1 and n > 0 fi > 1. Therefore, sgn(ri) = sgn(f i - l ) . (2.38) Let F (Tt) = V(iJi) - 1 = fi - 1 + fffajsgnfa) (2.39) where ^ ) - ^ ( b . ( M ) - M + I | . (2.40, We want to find a TJ such that F ( T J ) = 0 with i? < |TJ | < -B/V. Note that for r, > 0 15 Chapter 2. Governing Equations so that F is strictly decreasing on (B, oo) with F(B) = fi — 1 > 0. This result implies that we can solve for TJ numerically for example by using a bisection method on F. Furthermore, F(ri) is an odd function so that if TJ < 0, F'fa) is strictly decreasing on (—oo, —B) with F(—B) = fi — 1 < 0 so that we can also use the bisection method for TJ < . —B (see section 4.2.1). 2.5 Visualizing Regions I, II and III In this section, we would like to determine how to represent regions I, II and III on the (Rei, i?e2)-plane. Let us first begin by noting that we can express Re\ in terms of Re2 as itei = ^Re2. (2.42) Note that if we are in region III, then fi = 1 or Rei = qRe2. (2.43) If we are in region I with ^ > \ > 0, then we obtain the relation n - i = - ^ - I ) + ! • ( £ ) ) > G (-44) where G - ^ H ? H + I ) > 0 - (2-45) Since fi > 1 + G > 0, equation (2.44) can be rearranged to read 1/fi < 1/(1 + G). Therefore, there is no plug if 16 Chapter 2. Governing Equations Rei < T^QRe2- (2-46) Consider the other no plug case where g < — ̂  < 0. Equation (2.33) gives 1 _ n _ _ ^ f J j j ! ( i r - _ i ) + l n ( £ ) ) > e (2.47) which can be arranged to read Q < 1 — G < 1. It follows that Re2 = -Rex < -—-Rei (2.48) The no plug region with TJ < —B/rj2 will fall into one of the three cases below: 1 - G > 0 : Rei > ~~^Re2 (2.49) 1 — Cr 1 - G < 0 : i?ei < -r^Re2 (2.50) 1 — G 1 - G - 0 : i?e2 < 0. (2.51) Figure 2.1 shows where the regions fall on the (jRe2,i2ei)-plane for different values of B with T] = 0.883. It is clear that as the Bingham number increases, region II increases. Note that region III is a line. It is the rigid rotation line Rei = rjRe2 and it stays the same for all B. For a fixed value of B, figure 2.2 shows that region II decreases as rj increases. 2.6 Limit as B -> 0 In this section, our goal is to show continuity of the flow with respect to B. As B —> 0, Ry —> oo so that R0 = R2 which implies that we need to take the limit of the velocity profile in region I. Rewrite the velocity profile as V(r) = Air H Br ln(r)sgn(Tj) r 17 (2.52) Chapter 2. Governing Equations B=5TI=0.883 2000 1500 1000 500 / . / / I. No plug ' / Ml. Rigid Rotation / // \ / / I -4000 -3000 -2000 -1000 0 Re, 1000 B=5011=0.883 2000 1500 II. Plug on wall III ^ 1000 \ V 500 ^ - - ^ II n -4000 -3000 -2000 -1000 0 Re, 1000 B=10r|=0.883 2000 \ 1500 \ \ II. Plug \ on wall 1000 I. No plug \ / 500 III. Rigid Rotation n -4000 -3000 -2000 -1000 0 Re 2 B=100 n=0.883 1000 2000 1500 II. Plug on wall 1000 500 III ^ 7 " n i r ~r—-—^—— -4000 -3000 -2000 -1000 0 1000 Re, Figure 2.1: Regions with 77 = 0.883 for various values of B. Region III is the line Re\ = rjRe^- B=5 r|=0.6 B=5 n=0.7 -4000 -3000 -2000 -1000 Re, 0 1000 4000 -3000 -2000 -1000 Re, 1000 B=511=0.8 B=511=0.9 2000 1500 1000 500 0 V 1 1 II. Plug on wall \ '" V 1. No plug \ / " 2000 1500 * 1000 -4000 -3000 -2000 -1000 Re, 0 1000 500 0 1 1 i 1 "/ III i / I. No plug -4000 -3000 -2000 -1000 Re, 0 1000 Figure 2.2: Regions with B — 5 for various values of rj. 18 Chapter 2. Governing Equations where M = R\ + 2 / f + 5 1 n i ? 2 S g n ( r , ) (2.53) and 1 2 A2 = --TiRl (2.54) First, let us determine the relationship between Cl and TJ as B —>• 0. From equation (2.32) it follows that 1 - n = Urn { ^ R i f o 2 - 1) + BRX In sgnfa)} = | ( 7 - ^ ) ^ " ^ < 2- 5 5) Therefore as B —> 0, we obtain 2 v v J W - i . Taking the limits as B —> 0 of A\ and A2 and using equation (2.56) gives (2.56) Mm = (^ -4V -5- ) (2.58) so that as B —> 0 we obtain the velocity profile for < r < which is that of a Newtonian fluid in Couette flow [101. 1—7/ — — 1 — ?7 L J 19 Chapter 3 Linear Stability Analysis To determine for what values of Rei and Re2 the primary instability occurs, we will follow the method in [6] by adding a small perturbation to the basic flow, and determining whether this perturbation grows or decays with time. By assuming a normal mode solution for this perturbation, the linearized perturbation equations become a system of ordinary rather than partial differential equations, which leads to an eigenvalue problem. Furthermore, boundary conditions for the perturbation will be derived. In addition to the no slip boundary conditions on the cylinder walls, a boundary condition needs to be applied to the plug region in the case of type II flow. 3.1 Linear Perturbation Equations For now, consider some arbitrary basic, steady flow (U , P) where U is the basic velocity and P the basic pressure. . We would like to examine the linear stability of U by introducing a small disturbance in the velocity of the form eu' = (eu'(r, 9, z, t), ev'{r, 9, z, t),w'(r, 9, z, t)) (3.1) and in the pressure ep'(r,9,z,t) (3,2) where c « l . The perturbed flow is expressed as U + eu', P + ep'. Both the basic flow and the perturbed flow satisfy the Navier-Stokes and the continuity equation. That is, the perturbed 20 Chapter 3. Linear Stability Analysis flow satisfies V - ( U + eu') = 0 (3.3) eu't-f i ? e i ( U + eu') • V ( U + eu') = - V ( P + ep') + V • T ( U + eu') (3.4) and the basic flow satisfies V - U = 0 (3.5) i t e i ( U - V ) U = - V P + V - r ( U ) . (3.6) Subtracting equation (3.5) from (3.3) gives the perturbed continuity equation „ , du' u' 1 dv' dw' „ .„ V • u' = — + - + - — + — = 0. (3.7) or r r oo oz Subtracting equation (3.6) from (3.4) gives e(u't + Rei(u' • V ) U + ( U • V)u') + e2i?ei(u' • V)u' = -eVp + V • ( r ( U + eu') - r ( U ) ) . (3.8) Note that r is determined only in regions where T > B. Suppose we are in a region where both T ( U + eu) > B and r ( U ) > B. Then we have, Tij(U + eu') - 7 y ( U ) = Me(U + e u ' ) 7 I J ( U + eu') - / i C ( U ) 7 0 - ( U ) (3.9) where fie = 1 + ? • (3-10) 21 Chapter 3. Linear Stability Analysis Using the fact that jij(V + eu') = 7y (U) + 67^(11'), we linearize pe(U + eu') about the basic flow U to obtain M e ( U + eu') = / ie(U) + e V J 7 m n ( u ' ) % ^ + 0(e2). (3.11) mn ° l m n Note that 9 7 m n dj d-ymn &y \2 7 / 2 #7 Substituting this result into equation (3.11) gives pe(V + eu') = /J e(U) - e ^ — r g ^ 7 m n ( u ' ) 7 m n ( U ) + 0 ( e 2 ) = M e ( U ) - € 7 7 f j r 3 7 r e ( u ' ) 7 ^ ( U ) + 0(e2). (3.13) Substitute this result into equation (3.9) to obtain 22 Chapter 3. Linear Stability Analysis Tij(V + eu') - TijiU) = ( M e ( U ) - e^yrB(u')jre(V)j ( 7 0 -(U) + € 7 i j(u')) - Me(U) 7 i j (U) + 0(e 2) = 6 ̂ e ( U ) 7 y ( u ' ) - ^7y(U)7rtf(U) 7rfl(u')) + 0(e 2) Me(U)7ij(u'), ij^r9,6r = e < . , M ) ^ ( U ^ C u O - B ^ , . ij = re,Or A*c(U)7ij(u'), ij^rO,Or = e < + 0(e z) 7ij(u') ij = r6,0r e< 7(U) +0(e 2 ) 7ij(u') ij = r0,9r _ e < 7«K, + ™ « , ^ r * A + 0 ( £ 2 ) ( 3 1 4 ) [ 7ij(u') = r0,(9r where M;,- = T. . After some algebra, we find that 7 ( U ) , Y7TT . T T „ , v (du' \ ( ,dv ,v vdv'\ (vdw'\ Substituting this result, and our expression for T^(U + eu) - r^(U) into equation (3.8), and retaining terms of order e gives du' „ V (du' n , . <V 1 9 . . . ... 1 d . . . ... a .. , 1 . +3 ( - — (rMrr) + — {Mrz) - -Mr0 1 r ar az r (3.15) 23 Chapter 3. Linear Stability Analysis dv' „ / ,dV ,V Vdv' — + Reii u'— + u ' - + —— ot \ dr r r dv I dp' 1 0 , 2 . , 1 3 . . . d ,. , +B (1-^{Mee) + -fz(M6z) (3.16) aw' „ (Vdw'\ dp' 19. . . 1 9 , , a.. , +jB (^{rM-} + lii{Mze) + i { M z z ) ) • (3-17) Simplifying the above equations gives the following linearized perturbation equations dv' _ / ,_ Vdv' - + R e i (u'D*V + - - 1 dp' „•> , 2 6V r 90 +5 rI(I»2Hte2 + * S ^ l } (3.19) j [r d9 dz ) dw' „ V dw' dp' _2 , -^- + Re1-— = - - f + W at r cw where D* — 4- + \ and 7 = 7(U) where 7 = 17*1 = ^ - * (3.21) 24 Chapter 3. Linear Stability Analysis As required, 7 is nonnegative on [Ri R0] since it is decreasing on that interval with j{Ri) = \T{\ — B > 0 and j(Ry) = 0 or 7(^2) > 0. Note that if we are in region II, there is a singularity in equations (3.18)-(3.20) since in the denominator we will have j(Ry) = 0. This singularity, however, is remedied by the boundary conditions (3.60) that we will derive in section 3.3. These boundary conditions help to define the linear perturbation equations as r —>• Ry. 3.2 The Method of Normal Modes This system of PDEs can be reduced to a system of ODEs by applying the method of normal modes. That is, by assuming a solution of the form («', v', w',p') = ext+i(m9+kz\u{r), v(r),w(r),p(r)). (3.22) Let D = and substitute the above into each strain rate component 7 ^ to obtain ^rr(u') = 2 ^ = 2Duext+i^me+k^ (3.23) or ^(u'» - \{% ~"')+ % = (I(im"" ">+Dv) eA'+'M+") (3'24) 7r»M = ~ + ̂  = (iku + Dm) e*+«™»+*«> (3.25) oz or 7«(u ') = * (% + « ' ) = \ (imv + u) (3.26) ~M(U') = 2^- = 2ikwext+i(m6+kz\ (3.28) oz 25 Chapter 3. Linear Stability Analysis Let a = Xt + i(mO + kz) (3.29) and substitute all of the above into the Bingham terms of the linearized perturbation equations (3.18) - (3.20). The Bingham terms in equation (3.18) become 1 d (rirrW)\ + 1 fdjrzW) 7*«(u') r dr \ 7 / 7 \ dz r 1 d {2rDuea\ Id..., „ , „, 2 ,. , a r dr \ 7 / 702 rzi fid [2rDu\ -k2u + ikDw 2 ,. A a ,„ O A , = + : + e Q (3.30) \ r or \ 7 / 7 r̂ 7 / while for equation (3.19) the Bingham terms can be expressed as 1 fldjeoW) , dj6z(u') 7 \r d6 dz 1 d (2.. , A 1 d ( (1. ., \ a — 7 ^ -hmv + u)e + — — -zmw + z/cu e r^ydO \r ) j dz \\r J — -Trhmu — m v) k v e (6.61) 7 \ H r J and for equation (3.20) they can be rewritten as 1 d frizr(u')\ + 1 (ldjzoW + dizz(u') r dr y 7 J J \ f d6 ' dz ld_(r(iku + Dw)ca\ + J_j9 / / t m u ; + i j f c \ A + i^_t2ikwea) r dr \ 7 J rj dO \ \ r J J 7 dz na Miku + p^n + ± _ k \ +1 2 \ e „ \r or \ 7 / r 7 \ r / 7 / 1 3 / r ( » f c n + Z?ii;)\ _ 1 / m 2 ^ + fcrm; + 2 f c 2 \ \ g Q r dr \ 7 / 7 \ r 2 r / / Define 26 Chapter 3. Linear Stability Analysis <f>r(r,u,v,w,Du,Dv,Dw) = ±D (^p) + \ [ikDw - k2u - U)^j (3.33) „ „ ^ , 1 /2(imu — m2v) kmw , 9 \ , (f>e(r,u,v,w,Du,Dv,Dw) = - - i = '- k2v) 3.34 7 \ rz r J , , . 1 _ (r(iku + Dw)\ 1 (m2w kmv „ , 9 \ .„ „ , <j)z(r,u,v,w,Du,Dv,Dw) = -D f r - J - v ( - ^ + - " " - + 2k2w\ .(3.35) Then, substituting expressions (3.22) and (3.33)-(3.35) into the linearized perturbation equa- tions (3.7) and (3.18)-(3.20) gives V, , „ „ o Du m2u , 9 2imv u „ , .„ „„. Au + Rei — {imu - 2v) = -Dp + D2u + — ^ - - Art* —2 ^ + Bcpr (3.36) „ / ^ Tr V. \ im ~ Dv m2 ,2 \v + Rei [ uD*V-\ imv) = p + Dzv-\ Tv - k^v \ r ) r r r z 2im v „ , .„ + — 5 - u - - ~ + . B ^ (3.37) DID 77i^ Xw + i?ei —imw = —ikp + Z?2u; H TTW — k2w + Bct>z (3.38) THV D*u + — + ikw = 0. (3.39) r 3.3 Boundary Conditions for the Perturbation The boundary conditions for the perturbation are divided into two cases: no plug and plug on the wall. 3.3.1 Region I: No Plug From the no-slip conditions at the cylinder walls, it follows that 27 Chapter 3. Linear Stability Analysis u = v = w = 0 at r = Rur = R2. (3.40) Equation (3.39) implies that we can rewrite the boundary conditions as u = Du = v = 0 at r = Rur = R2. (3.41) 3.3.2 R e g i o n II: P l u g o n t h e W a l l In addition to the zero boundary conditions at r = R\ and r = R2, we must apply boundary conditions at the yield surface r = Ry. Since the velocity and pressure are being perturbed, then the stress is being perturbed, so the yield surface must also be perturbed. Since the perturbations are periodic, we assume that the yield surface has also been perturbed in a periodic manner, so we write the perturbed yield surface as r = Ry + eh{0, z,t) (3.42) where h = fi0ext+t^m6+kz^ and h0 is a constant. The yield criterion along with continuity of stress suggests that 7(U + eu') = 0 at the perturbed yield surface. That is, jij (U + eu') = 0, r = Ry + eh, i,j = r, 9, z. (3.43) Let tij = jij(V(Ry + eh) + eu'(Ry + eh, 6, z,t)). (3.44) Evaluating the perturbed strain rate at the perturbed yield surface, and linearizing about r — Ry , it follows from equations (3.23)-(3.28) that r r r = 2eDu(Ry + eh)ea = 2eDu(Ry)ea + 0(e2) 28 (3.45) Chapter 3. Linear Stability Analysis trz = e {iku(Ry + eh) + Dw(Ry + eh)) ea = e (iku(Ry) + Dw(Ry)) ea + 0(e2) (3.46) tgz = e f i m ^ + f * + t i few^ + e/*)) eQ = e f i m ^ ^ + ikv(Ry)) ea + 0(e2) (3.48) tzz = 2eikw(Ry + eh)ea = 2eikw{Ry)ea + 0(e2) (3.49) tr6 = e( im^^M + Dv(Ry) - ^ + h0D2V(Ry)) ea + 0(e2). (3.50) V Ky Ky J Note that in equation (3.50) we used the fact that j(Ry) = DV(Ry) — V ^ y ^ = 0. Using equation (3.43) we will set each = 0. Assuming k ^ 0 and retaining terms of order e gives 0 = Du(R-) (3.51) 0 = ikw(R-) (3.52) 0 = imw(Ry) + ikRyv(R~) = ikRyv{R~) (3.53) 0 = imv(R-) + u(R-)=u(Ry) (3.54) 0 = iku(Ry) + Dw{Ry) = Dw(Ry) (3.55) 0 = im^-^L + Dv(Ry) - + h0D2V(Ry) = Dv{Ry) + h0D2V{Ry). (3.56) Ry Ry Prom the above equations we see that u(R~) — 0. We will show that u(Ry) = 0. By continuity of velocity, we must have U(i?~) + eu'{R~,0,z,t) = U(i?+) + eu'(R+,0,z,t). Since the basic 29 Chapter 3. Linear Stability Analysis flow U is continuous everywhere between the cylinders, we have U(R~) = XJ(R+). It follows that u'(R~,0,z,t) = u'(R+,0,z,t) which implies that u(R~) = u(i?+). Furthermore, we can show that u = 0 throughout the plug. The boundary conditions at the yield surface are given by U(Ry) = v(Ry) = W(Ry) = 0 (3.57) Du(R-) = Dw(R-) = 0, Dv(R-) = -h0D2V(Ry). The above boundary conditions together with the fact that n Du u . (Dv v i . , , „ D2u + iz = -tm 5 - ikDw, (3.58) rp y2 \ rp rp2, which follows from the continuity equation, imply the boundary condition D2V(R-) D2u{Ry) = imh0 jj-2-L. Therefore, the boundary conditions in terms of u and v are given by i{Ry) = Du(R-) = 0, D2u(Ry) = imh0D2V(R-)/Ry v{Ry) = 0, Dv(Ry) = -h0D2V(Ry). (3.59) (3.60) 30 Chapter 4 A x i s y m m e t r i c Disturbance 4.1 Axisymmetric Equations If m = 0, equations (3.36)-(3.39) become An - 2Rex (^j v = -Dp + D2u + ^ ~ k 2 u - ^ + B<f>r (4.1) Xv + RelUD*V = D2v + — - k2v - - k2BV- (4.2) Xw = -ikp + D2w + ^ - k2w + B(pz (4.3) D*u + ikw = 0 (4.4) where 1 „ (2rDu\ 1 / , „ , 9 2u\ , , r D [—) + j \ % k D w ~ k u ~ 72 J (4-5) 1 _ (r{iku + Dw)\ 2k2 w ,„ „. -D : — . (4.6) r \ 7 7 7 ^ , and £)*D = j^i + \ J : , we can rewrite equations <f)r(r,u,v,w,Du,Dv,Dw) = (f>z(r,u,v,w,Du,Dv,Dw) = Using the fact that DD* = £z + ± £ - (4.1) through (4.4) as 31 Chapter 4. Axisymmetric Disturbance (DD* - k2 - X)u + 2Rei [-y)v = Dp- B<f)r (DD* —k— X)v = Rei(D*V)u + k2B- 7 (D*D — k2 — X)w = ikp - Bcbz D*u = —ikw. Substitute equation 4.10 into equation 4.9 to obtain (D*D — k2 — X)D*u = -ik(ikp - B<pz) = k2p + B$z where ( 1 1 $ = ikd)z [ u,v,——D*u,Du,Dv,——DD*u \ ik ik ik ( 1 D ( r ^ i k u + -rkDD*u)\ 2k2(-±D*u) V \ ^ / i . 1 D/r(k2u + DD*u)\ | 2k2D*u r \ 7 7 7 = a^Dzu + a2D2u + a\Du + OQU and the parameters 00-03 are given by 32 Chapter 4. Axisymmetric Disturbance k2 k2Dj Di 1 a0 = — + -TO 0 7 2 ~ " 3 ^ ( 4 - 1 2 ) ai = — + — + ^ (4.13) a2 = - 4 + ^ (4.14) a3 = - I . (4.15) 7 Differentiate equation 4.11 and solve for Dp to obtain Dp=^D(D*D-k2-X)D*u-^BD^z. (4.16) Substitute the above expression for Dp into equation 4.7 to obtain (DD* - k2 - \)u + 2Rei^v = ~D(D*D — k2 - X)D*u - ^BD^Z - BG>r (4.17) where $r = (j)r \u,v,—^-D*u,Du,Dv,—^-DD*u \ ik ik 1 n (2rDu\ 1 / ., 1 ^ , 9 2u\ = -D —.— + - -ik—DD*u -k2u--^\ r \ 7 / 7 \ ik rl J 1 „ (2rDu\ 1 / „ 2uN = — J - - DD.u + k2u+^T r V 7 / 7 V r , 7 / \ r 7 7^ / \ r 2 7 7 Note that ^D(D*D — k2 — X)D*u - (DD* — k2 - X)u = ^(DD* — k2 - X)(DD* - k2)u (4.18) so that we obtain the coupled differential equations 33 Chapter 4. Axisymmetric Disturbance (DD* - k2 — \)(DD* - k2)u = 2k2Re1 (^j v + B (D$z +{k2$r) (DD* — k2 - X)v = Rei(D*V)u + k2B^ (4.19) (4.20) which can be rewritten as the linear system (MNewtonian + BMBingham) = A DD* — k2 0 (4.21) where MNewtonian — ' (DD* - k2)2 -2k2Rex ^ ^ -Rei(D*V) DD* - k2 1 (4.22) and MBingham 0 D 0 k2 T / (4.23) Here, D is the operator which satisfies Du = D$z + k2$r. From equation (4.21) we see that the method of normal modes leads to an eigenvalue problem with eigenvalue A. 4.2 Numerical Methods In our numerical analysis, given the outer cylinder Reynolds number Re2, we would like to determine the inner cylinder Reynolds number Rei at which point the primary instability first occurs. If we are in a region where there is a plug on the wall, we will only consider the interval [Rl Ry]. 34 Chapter 4. Axisymmetric Disturbance 4.2.1 B a s i c F l o w E q u a t i o n s All of the basic flow equations can be fully determined by the set of parameters (77, B, Rei,Re2)- Before we calculate any of the flow equations, we first determine the region of interest (I or II) which is dictated by these four parameters as seen in section 2.5. If we are in region I, then T J is given by equation (2.35). Otherwise, if we are in region II then we apply a bisection method to equation (2.39) searching for a root T J such that F(T{) — 0. Since we are in region II, the root T{ must satisfy B < | T J | < B/rf. After having found T J , we can easily find Af(r\Ti, 77, B), Ro(i~i,rj,B) and V(r\Ti, 77, B, Rel, Re2). 4.2.2 F i n i t e D i f f e r e n c e M e t h o d Once all the basic properties of the flow have been determined, system (4.21) is fully defined. To approximate the derivatives, we will use a finite difference method over the yielded region [Ri R0] by taking N data points r\, r2, • • • rjv where ri = Ri + (i-l) N _ 1 • This mesh leads to an (JV — 2) x (N — 2) matrix approximation for each derivative (see Appendix A). Hence, system (4.21) becomes a (2N — 4) x (2iV — 4) eigenvalue problem. 4.2.3 F i n d i n g t h e C r i t i c a l f t e i To find the critical inner cylinder Reynolds number, Re\CTit, we employ a bisection method on Rei. Suppose we are working with 77 and B fixed. Then for each set of values (Rei,Re2, k) we obtain a spectrum of eigenvalues A. Let Am o x(A;) be the eigenvalue with the largest real part for a specific value of k. Suppose that \crit is the eigenvalue with the largest real part over all wave numbers k. That is, real(Amax(A;)) < real(Acrjt) for all k. Furthermore, suppose that real(Am a a ;) attains its maximum at k = kcru. Then Reicrit should satisfy rea\(\crit(Reicrit,Re2, kcrit)) = 0. (4.24) 35 Chapter 4. Axisymmetric Disturbance B=10, r|=0.883, N=50, (Re2,Re1)=(0,1080), real(Xcr|t)=-6.6379, kor| =2.6953, Region 2 10 - 0 —• k Figure 4.1: Re\ is too small. Instability has not yet occurred. If rea\(Xcrit(Rei, Re2, kcru)) < 0 then the perturbation will decay with time, so that Re\ is too small and the instability has not occurred yet in the system . Similarly, if real(Acrji(-Rei, Re2, kcrit)) 0 then the perturbation will grow with time, so that Rei is too large and the instability has already occurred for some smaller value of Re\ . For example, suppose we choose the set of values B = 10, rj = 0.883, Re2 = 0 with an initial guess of Re\ = 1080. Figure 4.1 shows that real(Acrj/,) < 0 so that we need to increase Re\. As a second attempt, we guess a larger value for Re\, say Re\ = 1180. Figure 4.2 shows that real(Acrjt) > 0 so that Re\ is too large. Finally, choosing Re\ = 1131.3351 results in a "small enough" value for real(AcrJi) as seen in figure 4.3. Finally, we plot the pair of points (Re2,Rei) = (0,1131.3351). By continuing in this manner for several values of Re2, we obtain the entire marginal stability curve, as seen in figure 4.4 where we use N = 50. The relative error for Re\crit between N = 50 and N = 75 is approximately 0.18 % while the relative error between N = 75 and JV = 100 is approximately 0.054 %. Figure 4.5 shows a flowchart for the marginal stability code. 36 Chapter 4. Axisymmetric Disturbance B=10, TI=0.883, N=50, (Re2,Re,)=(0,1180), real(X.rit)=7.1027, k,rit=3.1278, Region 2 1 1 1 1 I r I l I I 1 1 L 0 1 2 3 4 5 6 k Figure 4.2: Rei is too large. Instability has already occurred. B=10, ii=0.883, N=50, (Re2,Re1)=(0,1131.3351), real(Xor|1)=0.0060712, kcri=2.9371, Region 2 k Figure 4.3: Rei = ReiCTit 37 Chapter 4. Axisymmetric Disturbance Marginal stability curve: ti = 0.883, B = 10, N = 50 2000 500 Figure 4.4: Marginal stability curve for B = 10, 77 = 0.883 InputB, n, Rej (B, n.Rej) Bisect to find Rej (B, % Rej, Rei) Find J<nt (B, rj, Rej, Rei, k) Fmdregiqn (B,n, Rej, Rei, k, region) Find-tj, Ro, %y (B, i), Ke% Rei, k, % V, Y); Calculate the central difference i matrices and solve the eigenvalue problem. Figure 4.5: Flowchart of marginal stability code 38 Chapter 5 Resul ts for an A x i s y m m e t r i c Dis turbance 5.1 Numerical Results 5.1.1 C o m p a r i s o n w i t h N e w t o n i a n R e s u l t s Our first goal was to reproduce the marginal stability curve for the linear instability of a Newtonian fluid. We were able to compare our curve for m = 0, rj = 0.883, B = 0 with the linear instability curve in [9] courtesy of Dr. Randall Tagg from the University of Colorado in Denver who provided us with the numerical data. From figure 5.1 we can see that our numerical results are in agreement. 5.1.2 R e s u l t s for B > 0 We will first begin by looking at results for a narrow gap with rj = 0.883. Figure 5.2 shows marginal stability curves for B — 1, B = 10, B — 20 and B — 100 versus the region as dictated by 77, Rei, Re2 and B. Note that above the marginal stability curve, the partial plug and no plug regions may change due to the instability. We can see clearly in figure 5.3 that as B increases, the system becomes more stable so that the yield stress acts as a stabilizing agent in the system. This result is in agreement with that of Graebel [8]. In figure 5.4 we plot the critical values of Rei, denoted by Reicrit, versus B for different values of Re2. We can see that Reicrit is an increasing function of B for a fixed value of Re2. Furthermore, this figure shows that as B -> 00, ReiCTit becomes progressively independent of Re2. That is, 39 Chapter 5. Results for an Axisymmetric Disturbance Marginal Stability Curve for a Newtonian Fluid, T\ = 0.883 190 180 160 130 — Numerical results o Tagg -120 -100 -80 -60 -40 -20 20 40 Figure 5.1: Comparison of our numerical calculations against data from [1] for m B=1 ti=0.883 B=10ti=0.883 -4000-2000 0 2000 4000 6000 8000 Re„ Figure 5.2: Marginal stability curves versus regions for rj = 0.883 with N = 50 40 Chapter 5. Results for an Axisymmetric Disturbance Marginal stability curve: T) = 0.883, N= 50 2000 1500h 500 Figure 5.3: Marginal stability curves rj = 0.883 and varying B. Reicrit becomes nearly constant as B —)• oo. For large B, we may assume a solution of the form aBq (5.1) and use a least squares approximation to solve for a and q. Only considering values of B > 200 and with r? = 0.883 we get a « e 4 0 4 6 7 and q « 1.2169. Figure 5.5 shows the approximation (5.1) along with a plot for Re\CTit when B = 1000. The dependency of Re\CTit on Re2 is clearly much weaker than with respect to B. Table 5.1 shows approximate values of a and q for different values of n. We now wish to see if this result holds for flows with a wider gap. As seen in figure 5.6, for a gap width of rj = 0.5, Re\CTit still becomes progressively independent of Re2 for large B. However, we note that Re\cHi is not necessarily an increasing function of B. For example, the curve with Re2 = 1000 shows that Re\crii first decreases from B = 0 and then increases again at around B = 2. Therefore, for a wider gap, the presence of a yield stress is not necessarily 41 Chapter 5. Results for an Axisymmetric Disturbance Marginal stability values of Re with n, = 0.883, N=50 10 10 Re„=0 10 10 10' B 10 Figure 5.4: Re\CTit becomes progressively independent of Re2 as B increases. .<p e- 1.5 B=1000,11=0.883, N=50 R elcri, „ Re, = aBq -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 500 1000 Re, Figure 5.5: Re\crit for B = 1000 with the asymptotic solution Re\crit = aBq; a = e4 0 4 6 7 , q = 1.2169. 42 Chapter 5. Results for an Axisymmetric Disturbance a Q 0.3 g 3.8687 1.1619 0.5 g 3.7996 1.1946 0.75 g 3.9206 1.2076 0.883 g 4.0467 1.2169 Table 5.1: Approximate least squares values for a and q for B > 200 with N = 50. a stabilizing factor in the flow. In fact, we see that for certain values of B it may destabilize the system. Figure 5.7 shows curves for n = 0.5 and various values of B. We can see that for Re2 > 0, some values of B are more destabilizing than others, but there is no obvious pattern. However, the yield stress is still a stabilizing factor for Re2 < 0. In figure 5.7 we note also that the rigid rotation line Re\ = r]Re2 acts as a lower bound for stability. We wiU show in section 5.3 that this lower bound holds for all values of B and rj. Intuitively, when the entire region between the two cylinders moves in rigid rotation, we do not expect any instabilities to occur. 5.2 Heuristic Results In this section we will apply Rayleigh's criterion for an inviscid flow to our velocity profile for a Bingham fluid. Rayleigh's criterion gives a lower estimate for stability using physical arguments. Firstly, we will give a brief outline of Rayleigh's argument. Euler's equation for an inviscid flow in dimensional variables is given by p \oT + { t ' ^)tJ) = ~ P P (5"2) where tj = Uer + Veg + Wez. Assuming an azimuthal axisymmetric flow gives the 9 component 43 Chapter 5. Results for an Axisymmetric Disturbance Marginal stability values of Re with n, = 0.5, N=50 — Re2=0 -e- Re2=1000 ~o~ Re2=-1000 Figure 5.6: For a larger gap, Re\crit is not necessarily an increasing function of B. Marginal stability curve: TI = 0.5, N= 50 1 1 1 1 r~ 1500 1000 500 — B=0 - e - B=1 - A - B=2 B=5 — Rigid rotation Figure 5.7: Marginal stability curves with 77 = 0.5 and varying B. 44 Chapter 5. Results for an Axisymmetric Disturbance dv Tydv uv ~dv ^ + U— + — + W— = 0 (5.3) dt or r oz which can be rewritten as ^ 1 = 0 (5.4) Dt K ' where ^ is the material derivative. It follows from equation (5.4) that angular momentum is conserved for an invsicid fluid. Let us return to our dimensionless equations. Now, consider two small fluid elements, one located at r\ and one located at r2 with R\ < r\ < r2 < R0- Suppose that the fluid element at r\ moves up to position r2. Then, according to the conservation of angular momentum, it's velocity at r2 is given by V(r2) = r-^l. (5.5) r2 According to the r-momentum equation (2.3), the local pressure gradient at r = r2 is given by ™fA=Rei(m)i. ( M ) or r2 In order for the fluid element to return to its position n , we must have r2 r2 which can be expressed as (r2V(r2))2>(nV(ri))2. (5.8) That is, we will have stability when 45 Chapter 5. Results for an Axisymmetric Disturbance 1500 71 = 0.883 n = o.5 1000 500 -1000 -500 0 500 1000 Re, -1000 -500 0 500 1000 Re, Figure 5.8: Marginal stability curve for a Newtonian fluid along with the Rayleigh criterion. d{rVf dr > 0. (5.9) Equation (5.9) is known as Rayleigh's criterion for an inviscid fluid [10]. It follows from this criterion that for a Newtonian fluid in Couette flow, the flow will be stable when f2 > rf or Rei < Re2 V (5.10) This analysis is not rigorous, but from figure 5.8 we see that the marginal stability curve for a Newtonian fluid asymptotically approaches the Rayleigh criterion as Re2 -» oo. Let us now apply the Rayleigh criterion to our velocity profile for a Bingham fluid. Again, this method will not be rigorous, but should indicate when the velocity profile has a tendency to be stable, neglecting viscous effects. It follows from equation (5.9) that stability occurs if 46 Chapter 5. Results for an Axisymmetric Disturbance 2 ( r M r ) ) ^ > 0 (5.11) where u(r) = V/r is the angular speed at any point r in the flow. We will assume that 0,2 > 0. It follows that oj(r) > 0 so that stability occurs wherever  d^J^ > 0. Assuming that TJ < — B we obtain = 2Axr + Br(l + 2ln(r)) > 2Axr + Br{l + 2 ln(J2i)). (5.12) dr Therefore, our stability criterion becomes 2Ai + B ( l + 21n(i2i))>0 (5.13) where A ^ ^ + ^ - B H R o ) . . (5.14) In figure 5.9 we plot the above criterion and note that that the rigid rotation is always stable. In fact, Rayleigh's criterion for large B and the rigid rotation criterion almost coincide. This result also implies that whenever TJ > 0, the flow will be stable. When fl = 0 we obtain the Raleigh criterion for a Newtonian fluid. 5.3 A Lower Bound for Stability The preceding analysis of the Rayleigh criterion indicates that the rigid rotation curve is always stable, and this result is also confirmed by our numerical results. Furthermore, intuitively we expect that when the cylinders are locked (Cl = 1) an infinitesimal perturbation will be unable to break the unyielded fluid. In this section, we demonstrate that the above information is correct by deriving a rigorous bound for linear stability. We will show that real(A) < 0 if the inequality 47 Chapter 5. Results for an Axisymmetric Disturbance H = 0.883 15001 i 1 1 1 1 r Re2 Figure 5.9: Solid lines: Rayleigh criterion. Dashed line: Rigid rotation Rei — nRe2\ rj = 0.883. ^imaX 5 (Ro ~ Ri) ~ min | l , |v} ( l + JL} < 0 (5.15) holds. This inequality is a sufficient but not necessary condition for linear stability. It serves as a conservative approximation for linear stability. We will proceed as in [7]. Let u, v and w be the complex conjugates of u, v and w respectively. To begin the derivation of equation (5.15) we multiply equation (4.1) by ru, (4.2) by rv, (4.3) by rw and integrate from Ri to R0 to obtain 48 Chapter 5. Results for an Axisymmetric Disturbance A / rvv dr = —Re\ I (rDV + V)vu dr JRi JRI + j \rvD2v + vDv - k2rvv - — J (5.17) -k2B / — dr JRi 7 rRo rRo A / rwwdr — —ik / rwpdr JRi JRi rRo + / (rwD2w + wDw — k2rwvj) dr (5.18) rRo +B I rw<f>z dr. JRi Add equations (5.16) - (5.18) to obtain A||u|| — Re\Iinertia ^viscous BlBingham (5.19) where u|| 2 = fR°r(\u\2 + \v\2 + \w\2) dr JRi and (5.20) hnertiai = (tVvu - r (DV - ^ vu^j dr (5.21) iBingham = ~ jR {^u<t>r ~~ ^ ^ j " + d r ( 5 > 2 2) 49 Chapter 5. Results for an Axisymmetric Disturbance viscous (5.23) with rRo I (ruDp + ikrwp) dr JRi (5.24) h = - Ri rRo ruD2u + uDu — k2r\u rvD2v + vDv - k2r\v\2 - — (5.25) j (rwD2w + wDw — k2r\w\2) dr. Ri We will proceed to simplify the integrals Imertiah hiscous and I Bingham using integration by parts and applying the boundary conditions (3.60) with m = 0. Then we will proceed to bound these integrals using the Cauchy-Schwarz and triangle inequalities. Finally, we will apply the Poincare inequality to find our conservative estimate. 5.3.1 S i m p l i f y i n g t h e In tegra ls We will show that both of the integrals IViscous and Isingham are real and positive, while the integral Iinertia is complex. Of interest is the real part of equation (5.19) because we would like to determine the condition for which real(A) < 0. Note that I\ = 0 since by the continuity equation (4.4) and by integration by parts h I (rDu + u)p dr + r (Du H — H i JRi V r pdr = 0. (5.26) We can rewrite io. as 50 Chapter 5. Results for an Axisymmetric Disturbance I2 = k2 J*° r\u\2 dr + J*° (H 2 + H2) dr - I2l - 7 2 2 where and rRo I22 = / (uDu + vDv + wDw) dr. (5.29) JRi Using integration by parts, we see that hi — [{ru)Du + (rv)Dv + (rw)Dw]^ rRo - / (D(ru)Du + D{rv)Dv + D{rw)Dw) dr JRi rRo = - r{\Du\2+ \Dv\2+ \Dw\2)dr JRi rRo — I (uDu + vDv + wDw) dr JRi rn0 = - / r\Du\2dr-I22 (5.30) JR, so that integral IviSCOus becomes Let us rewrite the Bingham integral I Bingham- Integrating by parts we obtain (5.27) rRo J 2 1 = / (ruD2u + rvD2v + rwD2w) dr (5.28) JRi Iviscous = {r\Dxx\2 + A;2r|u|2 + M 2 + M 2 ^ d r ( 5 3 1 ) 51 Chapter 5. Results for an Axisymmetric Disturbance rRo / ru(/)r dr + u dr V 7 J \ R L JRX v 7 ; + JRi 7 V kDw k 2 + ^)u] dr. (5.32) Note that if R0 = Ry then it follows from the boundary conditions (3.60) with m = 0 that v 2rDu 2(rD2u + Du) 0 lim — — = hm — = , . = 0 (5.33) r-+Ry 7 r-*Ry Dj D^yRy) so that 2rDu L"( 7 = 0 (5.34) H i and we obtain J ru<j)r dr = - J 4 (2r\Du\2 + r (k2 + ̂  |u|2 - ikruDv^j dr. (5.35) Through similar methods, we obtain / ™ ^ JRi 2rk2 7 \w\2 1 dr = \w r r-.Ro 7 /•Ho Jill -{Dw(iku + Dw) + 2k 2\w\2) dr R1 Ri 7 j ~ (ikuDw + \Dw\2 + 2k2\w\2) dr. (5.36) The Bingham integral becomes 52 Chapter 5. Results for an Axisymmetric Disturbance iBingham = \ (2 (r\Du\2 + + rk2(\v\2 + 2\w\2) + r |x | 2 ) dr (5.37) where x = ku- iDw. (5.38) Note that both IviSCous a n d IBingham a r e r e a l a n d positive, which implies that the imaginary part of A can only come from the inertial integral Iinertia- Let us proceed by only considering the real part of equation (5.19). First, define the following quantities v — VR + ivi (5.39) u — UR + iui. (5.40) We can rewrite the integrand of Iinertia &s 2Vvu - r [ DV - — )vu = 2V(vu - vu) - (rDV - V)vu r = -rireiuRVR + mv!) +i(AV + r-yre){viuR - vRui) (5.41) so that rRo real(Jjn e r tio) = / r{-jrd)(uRVR + u / U / ) dr. (5.42) JRi 5.3.2 B o u n d i n g t h e In tegra ls We will begin by bounding real(ljnertia)- First note that 53 Chapter 5. Results for an Axisymmetric Disturbance Khnertio) < / r|7 r0|(|u f l||wfl| + |UJ | |UJ | ) dr rRo < imax r(\uR\\vR\ + |u/||v/|)dr. (5.43) JRi We can now apply a bound to the integrals of both r|?j#||vft| and r|ttj||uj| as follows. Rewrite the integral of r l n ^ H ^ ! as rRo rR0 I rr / r\uR\\vR\dr = / \vR\\ D(yuR{y))dy dr JRi JRi \JRi < JR °\VR\ (Jr \D{yuR{y))\dy^ dr = jR \vR\ ^ \kyu>[\dyj dr lywildy^J dr rRo rRo = k I \vR\dr / \rwr\ dr JRi JRi rRo / r \ \ rRo i ~ KJR \VR\drjR R%r~2\wi\dr (R0\^ fRo i fRo i = k J J r2\VR\dr J r*\u>i\dr. (5.44) Using the Cauchy-Schwarz inequality with the inner product <f,9>= [R° f(r)g(r)dr (5.45) JRi gives <l,rl\vR\> < y/< 1,1 >\f<r2\vR\,A\vR\ > <l , r5 | i o / |> < ^<l^T>y <A\wi\,r^\wi\ > (5.46) (5.47) so that 54 Chapter 5. Results for an Axisymmetric Disturbance jR ° r\uR\\vR\dr<k 2 (R0 - RJWVRWWW^. (5.48) Similarly, we obtain j*° r|uj |M dr <k ' (R0 - î lMIIKII- (5.49) 'Ri We can rewrite equation (5.43) as R0\ * rea.l(Iinertia) < 7maz ( "j^l (Ro-Rl)(k\\vR\\\\wi\\+k\\vi\\\\wR\ 1 |2 < * (R°Y2(R n ^ ^ l M ' + ll^ll2 , k2\\vi\\2 + \\wR - Jma*{Tj (R°-Rl>{ 2 + 2 = \imax(^~y (R0-Ri)(k^v\\2 + \\w\\2)^ (5.50) From inequality (5.50), it follows that we need to bound IBingham and IViSCous below with the integrals fc2||i>||2 and ||w|| 2. In-order to bound IBingham, we first set IBingham > / ° (2 (r\Du\2 + ̂ ) + rk2(\v\2 + 2\w\2) + r | x | 2 ) dr. (5.51) Umax JRl \ \ r / J Note that \x\2 = fc2|u|2 + ik(uDiB — uDw) + \Dw\2. Applying the continuity equation to the ikr(uDw — uDw) term in the above integral, and integrating it by parts gives rRo rRo / \u\2\ ik J ((ru)Dw - (ru)Dw) dr = - J 2 \r\Du\2 + J dr. (5.52) Therefore, 55 Chapter 5. Results for an Axisymmetric Disturbance B IBingham > ^— {k2r(\v\2 + 2\w\2) + k2r\u\2 + r\Dw\2) dr Umax J Ri B [R° (k2r\v\2+ r\Dw\2) dr. > Imax JRi (5.53) The final and simplest integral to bound is Iviscous and we set rRo Iviscous > / (k2r\v\2 + r\Dw\2) dr. (5.54) From the inequalities (5.50), (5.53) and (5.54) it follows that A ^ i i t x H 2 < r ( ^ 0 - ^ 1 ) ( f c 2 i b i i 2 + i i ^ n 2 ) - f i + b i ' ^ ' . - " ^ » - • • « * I * Imax" {k2\\v\\2 + \\Dw\\2). (5.55) 5.3.3 Applying the Poincare Inequality We now need to find a lower bound for | | J D U J | | 2 involving ||u;|| 2. It follows from the Poincare inequality that f1 ^ dy>*2f\W\2dy. Jo dy J0 Introduce a change of variable r = R\ + {R0 - Ri)y- Then we obtain f R ° \ d ^ 2 d r = 1 f1 JR1 \dr R0- Ri Jo dw dy dy n2 f1 ~ R0 — Ri Jo Wo \w\2 dy IT2 rR° (5.56) (5.57) 56 Chapter 5. Results for an Axisymmetric Disturbance Therefore where j r\Dw\2dr > Ri \Dw\2 dr (5.58) JRi JRX > — T r 2 / r\w\2dr (5.59) rRo > min{l,C} / r\w\2dr (5.60) JRi C = I-TT2.' (5.61) Substituting \\Dw\\2 + & 2 | |u| | 2 > min{l, C } ( | H | 2 + & 2 |M| 2 ) into equation (5.55) gives A*||u||2 < (|^) 5 (Ro - Ri) ~ min{l, C} ( l + JL^J (k2\\v\\2 + H | 2 ) (5.62) so that we will have stability when XR < 0 or whenever equation (5.15) holds. 5.3.4 R e s u l t s Let f(ReuRe2) = ^ 7 m a i ( ^ j ' (R0 - Rx) - min |l, | v j (l + JL) . (5.63) The plot of f(Rei, Re2) - 0 is shown in figure 5.10 for B = 10 and B - 100. We can see that the rigid rotation line Rei = rjRe2 is contained inside this contour, so that analytically we have proven that the rigid rotation is stable. 57 Chapter 5. Results for an Axisymmetric Disturbance 58 Chapter 6 Conclusions Using finite difference methods we were able to determine marginal stability curves for various values of B and 77. In the case of a narrow gap, our results agree with [8] in the sense that the yield stress acts as a stabilizing agent in the flow. As the Bingham number increases, the system becomes more stable. However, in the case of a wider gap, the yield stress can destabilize the flow when Re% > 0. Regardless of the gap width, we found that as B —» oo, the critical inner cylinder Reynolds number Reicrit becomes weakly dependent on Re2- That is, Re\ depends primarily on B and 77. Using numerical methods, heuristic arguments and analytical methods, we were able to show that the rigid rotation flow is always stable and acts as a lower bound for stability for any values of B and 77. This result implies that whenever the stress at the inner cylinder is positive or CI2 > fii5 the flow will be stable. Future considerations for this problem could include analysis with the m > 0 case and deter- mination of the partial plug region after the onset of instability. 59 Bibliography [1] C D . Andereck, S.S. Liu and Harry L. Swinney. Flow regimes in a circular Couette system with independently rotating cylinders. J. Fluid Mech., 164, pages 155-183, 1985. [2] R. Byron-Bird, G.C. Dai and Barbara J. Yarusso. The Rheology and flow of viscoplastic Materials. Reviews in Chemical Engineering, 1, pages 1-70, 1982. [3] S. Chandrasekhar. Hydrodynamic and hydromagnetic stability. Oxford, Clarendon Press, 1961. [4] P. Chossat and G. Iooss. The Couette-Taylor problem. Springer-Verlag, 1994. [5] B. Engelmann, R. Hiptmair, R.H.W. Hoppe, G. Mazurkevitch. Numerical simulation of electrorheological fluids based on an extended Bingham model. Computing and Visualiza- tion in Science, 2, pages 211-219, 2000. [6] I. Frigaard and C. Nouar. On the stability of Poiseuille flow of a Bingham fluid. J. Fluid Mech., 263, pages 133-150, 1994. [7] I. Frigaard and C. Nouar. On three-dimensional linear stability of Poiseuille flow of Bing- ham fluids. Physics of Fluids, 15, pages 2843-2851, 2003. [8] W.P. Graebel. The hydrodynamic stability of a Bingham fluid in Couette flow. Interna- tional Symposium on Second Order Effects in Elasticity, Plasticity and Fluid Dynamics, Haifa, Pergamon, pages 636-649, 1964. [9] W. F. Langford, Randall Tagg, Eric J. Kostelich, Harry L. Swinney, and Martin Golubitsky. Primary instabilities and bicriticality in flow between counter-rotating cylinders. Physics of Fluids, 31, pages 776-785, 1988. [10] R. Tagg. The Couette-Taylor problem. Nonlinear Science Today, 4, pages 2-25, 1994. [11] G.I. Taylor. Stability of a viscous liquid contained between two rotating cylinders. Phil. Trans. Roy. Soc. London A 223, pages 289-343, 1923. 60 Appendix A Numerical Approximations of the Eigenvalue Problem We will approximate the derivatives of our eigenfunctions u and v by using a finite difference method. Let N > 0 be an integer and divide the interval [Ri,R0] into N — 1 subintervals each of length h = — ^ — ^ and set = Ri + (i — l)h for i = 1,2,..., N. The derivatives of u (or at any point will be approximated by the central differences Du* = — ^ (A.1) ^ = ( A 2 ) D \ , = ^ - ^ 1 + ( A 3 ) Ar D m = (A.4) where ui « ufo) and Du,- « Du(rj). The boundary conditions u = i> = 0 a t r = i? 0 satisfied by the eigenfunctions can be rewritten as u(r\) = « i = 0 (A.5) U(TN) — UN = 0 (A.6) v(n) = «i = 0 (A.7) «(rjv) = «iv = 0 (A.8) while the boundary conditions Du = 0 at r = R0 can be approximated as 61 Appendix A. Numerical Approximations of the Eigenvalue Problem Du(r i ) Du(rN) u2 - u0 = 0 2h 2h (A.9) (A.10) That is uo = u2 UN+I — UN-I (A.U) (A.12) Each of the estimated derivatives for u (A.1)-(A.4) along with the boundary conditions for u (A.5), (A.6), (A. l l ) and (A.12) can be expressed as the (N - 2) x (N - 2) matrices D D 1 u2 ^ «4 UN-2 UN-l U2 U3 U4 UN-2 UN-l ( 1 h 0 0 1 0 0 0 0 0 1 2 0 0 0 0 1 0 -2 1 0 0 0 \ \ 0 1 - 2 1 0 0 1 - 2 U2 U3 U4 UN-2 UN-l U 3 U4 UN-2 UN-l (A.13) (A.14) 62 Appendix A. Numerical Approximations of the Eigenvalue Problem D4 U2 tt 3 U4 U5 UN-2 y uN-i j U3 U4 U5 UN-2 \ UN-1 ) 1 i_ 1 2 -1 1 2 0 0 0 0 1 0 -1 1 2 0 0 0 1 2 1 0 -1 1 2 0 0 0 1 2 1 0 -1 1 2 0 0 0 1 2 1 0 -1 1 2 0 0 0 1 2 1 0 -1 0 0 0 0 1 2 1 1 2 7 - 4 1 0 0 0 0 -4 6 - 4 1 0 0 0 1 -4 6 -4 1 0 0 0 1 - 4 6 -4 1 0 0 0 1 -4 6 -4 1 0 0 0 1 - 4 6 - 4 0 0 0 0 1 -4 7 J U2 U3 Ui U5 UN-2 \ UN-1 ) U2 Uz Ui U5 UN-2 \ UN-1 J Note that for the eigenfunction v we will only require the first and second derivative matrices, and according to the boundary conditions (A.5) - (A.8) these will be the same ones as those that are used for u. 63

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 4 0
Russia 3 0
Canada 1 0
United States 1 0
City Views Downloads
Unknown 3 2
Beijing 3 0
Shenzhen 1 0
Sherbrooke 1 0
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items