50 Years of the K-BKZ Constitutive Relation for Polymers

Abstract

The K-BKZ constitutive model is now 50 years old. The paper reviews the connections of the model and its variants with continuum mechanics and experiment, presenting an up-to-date recap of research and major findings in the open literature. In the Introduction a historical perspective is given on developments in the last 50 years of the K-BKZ model. Then a section follows on mathematical modeling of polymer flows, including governing equations of flow, rheological constitutive equations (with emphasis on viscoelastic integral constitutive equations of the K-BKZ type), dimensionless numbers, and boundary conditions. The Method of Solution section reviews the major developments of techniques necessary for particle tracking and calculation of the integrals for the viscoelastic stresses in flow problems. Finally, selected examples are given of successful application of the K-BKZ model in polymer flows relevant to rheology.

1. Introduction

1.1. Rheology

Rheology is defined as the “study of deformation and flow of matter” [1, 2]. The term has ancient Greek roots that refer back to the 6th BC, when the Greek philosopher Heraclitus realized the relative change of all elements in his well-known motto “𝜏𝛼𝜋𝛼́𝜈𝜏𝛼𝜌𝜀𝜄́” or “everything flows” [3]. In our days the term “rheology” was first used in 1920 by the American chemistry professor Eugene Bingham in Lafayette College, Indiana, USA. Bingham consulted with colleagues in the Department of Classical Studies in his effort to explain the peculiar behavior of various colloidal solutions [4]. The term “rheology” and its above definition were accepted by the (American) Society of Rheology (SOR), founded in 1929 with its first president being Prof. Bingham. Many other national rheological societies have since come to being, with the European Society of Rheology (ESR) established in 1996 and encompassing many individual European societies. The various rheology societies celebrate every four years the advancements in rheology at the International Congress on Rheology. The last one took place in Lisbon, Portugal, in August 2012 [5].

According to the Heraclitian definition, the term “rheology” could be used for all materials, including the classical limit cases of Newtonian fluids, such as water, and elastic Hookean solids, such as metals. However, these limiting cases are often considered outside the scope of rheology, which deals mainly with materials characterized by complex behavior. As an example, the 1st annual meeting of SOR in 1929 in the USA included presentations on asphalt, lubricants, paint, plastics, and rubber, which gives an idea of the broad scope of the subject matter and the various scientific disciplines implicated in rheology. Today the subject of rheology has seen tremendous growth and includes polymer rheology (plastics and rubber), colloidal rheology, suspension rheology, biorheology, food rheology, and so forth [1].

Since rheology deals with flow and deformation, it has been customary to take as predominant the term “flow” (to flow means “𝜌^‵𝜀𝜔́” or “rheo” in ancient Greek) and study the flow (and deformation) of complex fluids, which are then called “non-Newtonian” fluids. The Journal of Non-Newtonian Fluid Mechanics (JNNFM) was established in 1976 and has been published ever since uninterruptedly by Elsevier [6].

Most non-Newtonian fluids under study in rheology exhibit highly nonlinear behavior under flow and deformation. This nonlinearity is known for polymeric materials as shear thinning (or pseudoplasticity) and governs the viscosity as a function of the shear rate. Other materials exhibit time-dependent phenomena of shear thinning and are called thixotropic. The presence of macromolecules in polymers is responsible for elastic phenomena, and this gives rise to viscoelasticity, which includes the relaxation time of the material.

When time of response of the materials comes into play, the traditional idea of solids and fluids may no longer be valid, and a different classification need be applied. This is taken into account when considering the rheological behavior of the materials, and the classification is more appropriate according to the rheology of the materials. The time constant present in the rheology of time-dependent materials gives rise to the dimensionless Deborah number (De). When De=0, the material behaves as a Newtonian liquid, while when De→∞, the material behaves as a solid. In-between values of De describe most of real materials.

Since rheology deals with the flow and deformation of matter in the most general case of time-dependent and three-dimensional space, it is necessary to describe such behavior by introducing differential calculus with tensor and vector notation in time and space. This makes the subject rather difficult and needs a good grasp of mathematics and physics. Some simple cases (usually one-dimensional) have analytical solutions. However, most cases require numerical solution, and this leads to the need for a good knowledge and experience with numerical analysis and computational methods. The presence of powerful computers and workstations has allowed the solution of many important problems in rheology. Still many more remain unsolved and are the subject matter of intensive research efforts.

The main current subjects of rheology deal (i) with rheometry (experimental methods for measuring the rheological properties of the materials), (ii) with constitutive equations (mathematical models for describing the rheology of the materials), and (iii) with complex flows (for visualizing and understanding the flow behavior of materials).

A major part of rheological research relates to polymers. Since the tremendous growth of the plastics industry, especially after the 2nd World War, it was realized that the making of plastic products involves the flow and deformation (hence rheology) of macromolecules (polymers). These were found to behave as highly non-linear non-Newtonian fluids, producing flow phenomena in sharp contrast to Newtonian fluids. The rheological community was then hard pressed (i) to find ways to measure the properties of these materials (hence rheometry), (ii) to model them with appropriate rheological models, and (iii) to predict their flow behavior in industrial operations, collectively known as “polymer processing.” This proved to be a daunting task, and many researchers around the world have contributed to a level of sophistication in all three areas above commensurate with the industrial importance of the subject.

An intellectually challenging endeavor was to formulate rheological equations of state (constitutive equations) able to describe well the unusual behavior of macromolecules (polymers). Many efforts were initiated, mainly at universities, and a plethora of constitutive equations saw the light with varying degrees of success. The present work refers to one such example of constitutive equations that saw the light 50 years ago (in 1962) and has since captured the imagination and research efforts of various scientific groups around the world. This is the K-BKZ integral constitutive equation [7, 8], which, arguably, has shown most promise in mimicking and predicting the behavior of polymer solutions and melts in flows of scientific and industrial importance. The present paper draws heavily on Tanner’s work [9] which celebrated the 25th Silver Jubilee of the K-BKZ model. It brings up-to-date information and more examples of the successful application of the K-BKZ model to polymer flows to help celebrate the 50th Golden Jubilee of this remarkable model.

1.2. The K-BKZ Original Papers

In 1963 the Transactions of the Society of Rheology (now Journal of Rheology) published the article known from the initials of its authors as the BKZ paper [7] (see Figure 1 regarding the three authors [2]). The article was presented to the Society of Rheology (SOR) at its annual meeting in October 1962; Alan Kaye’s note [8] expounding essentially the same ideas shares this date. Tanner in his review paper [9], celebrating 25 years of the model, charts the citations of these two articles from its birth. He depicted data taken from the Science Citation Index, covering the period 1963–April 1987. The two papers were closely linked and the rates of total citations were roughly in the proportion 140 (BKZ) to 50 (Kaye), each author seeming to attract about a quarter of the attention, which appears a fair result.

Figure 1 

The three authors of the BKZ article at the Rochester Society of Rheology meeting in 1991. From left: Barry Bernstein, Elliot Kearsley, and Louis Zapas [2].

These numbers have now been extended to the present date covering 50 years of the model. Taken from the same source (Science Citation Index, ISI) [10], the data up to 2012 give 501 citations for BKZ and 150 for Kaye, keeping roughly the 1/4 ratio for each author. The search for “BKZ” mentioned in an article brings out about 2000 papers (without self-citations by the authors) as shown in Figure 2. By all means the BKZ paper remains a citation classic, a status enjoyed by a mere fraction of published papers (in the rarefied world of 0.002%) [9].

Figure 2 

The Science Citation Index results for articles mentioning the BKZ model (ISI, [10]).

The pattern of citation is also quite abnormal. Generally, most citations occur within five or so years of publication and thereafter papers are quietly forgotten. Following the results presented by Tanner [9], the BKZ model got more attention in 2010 alone (220 citations) than in any other year since its presentation 50 years ago. Obviously, there must be good reasons for this, and in this paper more evidence will be given for the model’s use and its successes.

2. The K-BKZ Constitutive Model

From the inception of the basic ideas of a mathematical model to make it suitable for engineering applications, a long and arduous way has to be followed, along which a number of key people have contributed to make this happen. The following constitutes in the simplest possible terms a guide along that pathway that led from the original ideas to the current use of the model and its variants in engineering applications.

2.1. One-Dimensional Model

Before describing the general theory let us consider the one-dimensional (1D) shearing case, following Tanner [9]. Application of a step shear strain of magnitude Δ𝛾_𝑙 to an unstressed sample at time 𝑡₁ will result in a shear stress response 𝜏; if the response is linear in Δ𝛾_𝑙 we can write [1] 𝜏(𝑡)=𝐺𝑡−𝑡₁Δ𝛾₁, (1) where 𝐺(𝑠) is the relaxation modulus and 𝑠=𝑡−𝑡₁ is the time interval. Application of a series of steps at various times 𝑡₁,𝑡₂,…,𝑡_𝑁 results in the linear response 𝜏(𝑡)=^𝑁_𝑎=1𝐺𝑡−𝑡_𝑎Δ𝛾_𝑎 (2) and in the limit of a continuous strain history 𝜏(𝑡)=^𝑡_−∞𝐺𝑡−𝑡^′𝑑𝛾𝑡^′. (3)

We can integrate (3) by parts to find 𝜏(𝑡)=𝐺𝑡−𝑡^′𝑑𝛾𝑡^′𝑡_−∞+^𝑡_−∞𝛾𝑡^′𝑚𝑡−𝑡^′𝑑𝑡^′, (4) where 𝑚(𝑠)=−𝑑𝐺/𝑑𝑠 and 𝐺(∞)→0. The function 𝑚(𝑠) is usually referred to as a memory function as it depends on the time interval 𝑠=𝑡−𝑡^′. If we now measure 𝛾(𝑡^′) from the present configuration as reference, as is convenient in fluids, then (4) becomes 𝜏(𝑡)=^𝑡_−∞𝑚𝑡−𝑡^′𝛾𝑡^′𝑑𝑡^′. (5) One can look at (5) as a superposition of linear, noninterfering effects as before, known also as Boltzmann’s superposition principle. If we now permit 𝑚 to depend on the strain level, becoming a noninterfering superposition of nonlinear responses, then 𝑚=𝑚𝛾𝑡^′;𝑡−𝑡^′ (6) or letting 𝑚=𝜕𝑈,𝜕𝛾 (7) where 𝑈=𝑈𝛾𝑡^′;𝑡−𝑡^′ (8) then (5) becomes 𝜏(𝑡)=^𝑡_−∞𝜕𝑈𝛾𝑡𝜕𝛾^′;𝑡−𝑡^′𝛾𝑡^′𝑑𝑡^′ (9) which is the one-dimensional K-BKZ model.

2.2. Three-Dimensional Model

In the fully three-dimensional (3-D) case, for an incompressible fluid, the equation as developed in the original BKZ article [7] reads 𝜎^𝑖𝑗(𝑡)=−𝑝𝛿^𝑖𝑗+^𝑡_−∞𝑥^𝑖_𝑘𝑡,𝑡^′𝑥^𝑗_𝑙𝑡,𝑡^′𝜕𝑈𝜕𝐸_𝑘𝑙𝑡,𝑡^′𝑑𝑡^′, (10) where 𝜎^𝑖𝑗 is the total stress tensor, 𝑝 is the pressure, and 𝛿^𝑖𝑗 is the unit tensor. The notation for 𝑥^𝑖_𝑘 meant 𝑥^𝑖_𝑘𝑡,𝑡^′=𝜕𝑥^𝑖(𝑡)𝜕𝑥^𝑘𝑡^′. (11)

The notation introduced in (10) is not simple, and further useful definitions are required to make it easier to understand. Figure 3 shows a relevant notation of the position vector 𝐱, velocity vector 𝐯, and times 𝑡 (present) and 𝑡^′(past).

Figure 3 

Fluid particle trajectory (pathline) and relevant notation (in steady-state flows of viscoelastic fluids, the present state of the stress 𝝉(𝐱) is a function of the present kinematics [𝐯(𝐱),∇𝐯(𝐱)] as well as of those of the past states).

The right Cauchy-Green strain tensor is defined as 𝐶_𝑘𝑙𝑡,𝑡^′=𝛿_𝑖𝑗𝑥^𝑖_𝑘𝑡,𝑡^′𝑥^𝑗_𝑙𝑡,𝑡^′ (12) and the strain tensor 𝐸_𝑘𝑙 is defined as 𝐸_𝑘𝑙𝑡,𝑡^′=12𝐶_𝑘𝑙𝑡,𝑡^′−𝛿_𝑘𝑙. (13)

In (10) 𝑈 is a “potential function” and is assumed to be a function of 𝐸_𝑘𝑙 and (𝑡−𝑡^′). Furthermore, for the isotropic fluids considered here, 𝑈 is a function of the invariants 𝐼₁ and 𝐼₂ of 𝐄, where 𝐼₁𝐼=tr𝐄,₂=tr𝐄² (14) and 𝐄 is the strain tensor defined in (13). The third invariant of 𝐄 (𝐼₃=tr𝐄³) need not be an argument of 𝑈 for incompressible materials (being 1).

In (11) and (12) one notes the dependence of the strains on the 𝑥_𝑖 coordinates at two times, 𝑡 and 𝑡^′. We shall use the present (𝑥_𝑖(𝑡)) coordinates as a reference frame and write the equivalent of (11) as (Cartesian tensor notation, (𝑥,𝑦,𝑧)=(𝑥₁,𝑥₂,𝑥₃)) 𝐹_𝑖𝑗𝑡^′=𝜕𝑥_𝑖𝑡^′𝜕𝑥_𝑗(𝑡), (15) where 𝐹_𝑖𝑗(𝑡^′) is the relative deformation gradient tensor. With this definition the tensor 𝐶_𝑘𝑙(𝑡,𝑡^′) is equivalent to (𝐹⁻¹)_𝑘𝑚(𝐹^−𝑇)_𝑚𝑙 or the inverse of the Cauchy-Green tensor relative to the coordinates 𝐱(𝑡) as reference. We shall refer to this strain as the Finger strain tensor, 𝐂⁻¹(𝑡^′) [1].

The tensor 𝐄 is then just (1/2)[𝐂⁻¹(𝑡^′)−𝜹]. With these changes of notation, (10) becomes 𝝈=−𝑝𝜹+^𝑡_−∞𝐅⁻¹𝜕𝑈𝐅𝜕𝐄^−𝑇. (16)

This is the form of the original BKZ equation [7].

2.3. Kaye’s Contribution

Kaye [8] went straight for a generalization of the rubber-like liquid of Lodge [11]. With the notation used above, he wrote, for the rubber-like liquid, 𝝈=−𝑝𝜹+^𝑡_−∞2𝑆𝑡−𝑡^′𝐂⁻¹𝑡^′𝑑𝑡^′, (17) where 𝑆 is another memory function; he recognized that this is related to the simplest form of the theory of rubber elasticity [12] 𝝈=−𝑝𝜹+2𝜕𝑊𝜕𝐼₁𝐂𝑡,𝑡^′, (18) where 𝐂 is defined in (12). The simplest 𝑊 function is 1𝑊=2𝑁𝑘𝑇𝐼₁−3, (19) where 𝑁𝑘𝑇 is constant. In this case 𝝈=−𝑝𝜹+𝑁𝑘𝑇𝐂𝑡,𝑡^′. (20)

Equation (20) translates into (17) for a rubber-like liquid. The most general incompressible rubber-elastic theory [13, 14] was known by Kaye [8] to give 𝝈=−𝑝𝜹+2𝜕𝑊𝜕𝐼₁𝐂+𝜕𝑊𝜕𝐼₂𝐂𝐼₁−𝐂² (21) (here 𝐂 is measured relative to the undeformed configuration), where the corresponding rubber-elastic liquid was given as 𝝈=−𝑝𝜹+2^𝑡_−∞𝜕𝑈𝜕𝐼₁𝐂⁻¹+𝜕𝑈𝜕𝐼₂𝐂⁻¹𝐼₁−𝐂⁻²𝑑𝑡^′, (22) and 𝐼₁=tr𝐂⁻¹, 𝐼₂=(1/2)[(tr𝐂⁻¹)²−tr𝐂⁻²].

Kaye [8] also went on to write (22) as (in our notation) 𝝈=−𝑝𝜹+2^𝑡_−∞𝜕𝑈𝜕𝐼₁𝐂⁻¹𝑡^′−𝜕𝑈𝜕𝐼₂𝐂𝑡^′𝑑𝑡^′, (23) where 𝑝 is also a pressure, slightly different from that in (22).

Equation (23) appears to be the most convenient form of the constitutive equation [1]. Here we have 𝐼₁=tr𝐂⁻¹(𝑡^′) and 𝐼₂=tr𝐂(𝑡^′), and 𝑈 is a function of 𝐼₁, 𝐼₂, and 𝑡−𝑡^′. The transformation of the original form of (16) into the form of (23) did not appear in the 1965 BKZ article [15]; it appeared however in the article by Zapas and Craft [16] of 1965. Further developments for the nonisothermal case were given in 1964 [17]. Thus from 1965 on, the theory was essentially complete in the form now generally used. We now look at the variants of the model proposed by other researchers in an effort to make it tractable for numerical simulations.

2.4. Other Contributions

Tanner [9] describes how it is possible to trace the elements of K-BKZ theory in the continuum mechanics literature. It begins from the formulation [18, 19] which gives the stress tensor as a functional of C for a simple fluid over all past times. Pipkin [20] gives 𝝈=−𝑝𝜹+^𝑡_−∞𝑀₁𝑡−𝑡^′𝐂𝑡^′𝑑𝑡^′, (24) where 𝑀₁(𝑡−𝑡^′) is another memory function. For reasons beyond the scope of the present work, (24) with 𝐂 as a strain measure does not provide a realistic constitutive equation [1]. A better description is obtained if one substitutes 𝐂 with a function 𝐇(𝐂) to get 𝝈=−𝑝𝜹+^𝑡_−∞𝑚𝑡−𝑡^′𝐇𝐂𝑡^′𝑑𝑡^′. (25) Here the memory function 𝑚 is the same as that in (5) provided we scale 𝐇 appropriately; it is also equal to 2𝑆 in (17). The general form of 𝐇 is known to be, from rubber elasticity [1] (to within an isotropic tensor), 𝐇=𝑓₁𝐼₁,𝐼₂𝐂⁻¹+𝑓₂𝐼₁,𝐼₂𝐂, (26) where 𝐼₁=tr𝐂⁻¹, 𝐼₂=tr𝐂, and 𝑓₁ and 𝑓₂ are functions of these invariants. Hence (25) becomes 𝝈=−𝑝𝜹+^𝑡_−∞𝑚𝑡−𝑡^′𝑓₁𝐂⁻¹+𝑓₂𝐂𝑑𝑡^′, (27) where 𝐂 is evaluated at 𝑡^′ relative to the reference state at 𝑡. Choosing a strain measure 𝐇 optimal in some sense in the above equation is an idea with obvious appeal.

One can go further and note that (27) is the separable case (in the sense that the memory function m is independent of the strain measure 𝐇) of the more general rule 𝝈=−𝑝𝜹+^𝑡_−∞𝜙₁𝐂⁻¹𝑡^′+𝜙₂𝐂𝑡^′𝑑𝑡^′, (28) where 𝜙_1,2 are functions of 𝐼₁, 𝐼₂, and 𝑡−𝑡^′. This rule was first given by White and Tokita [21] and then later by Rivlin and Sawyers [22]. This generic integral model includes many other well-known constitutive equations, whose rheological response is studied in the monographs by Bird et al. [23], Larson [24], and Tanner [1]. For example, 𝜙₁=1 and 𝜙₂=0, and assuming separability for the memory function 𝑚(𝑡−𝑡^′) yields the Lodge rubber-like liquid model [25].

The factorized K-BKZ model proposed by White and Tokita [21] can be written for the extra stress tensor 𝝉 as 𝝉=^𝑡_−∞𝑚𝑡−𝑡^′𝜕𝑉𝐼₁,𝐼₂𝜕𝐼₁𝐂⁻¹𝑡^′+𝜕𝑉𝐼₁,𝐼₂𝜕𝐼₂𝐂𝑡^′𝑑𝑡^′, (29) where there is separability between the time-dependent function 𝑚 and the strain-dependent function 𝑉.

The Rivlin-Sawyers variant of the K-BKZ model is written as 𝝉=^𝑡_−∞𝑊₁𝑡−𝑡^′,𝐼₁,𝐼₂𝐂⁻¹𝑡^′+𝑊₂𝑡−𝑡^′,𝐼₁,𝐼₂𝐂𝑡^′𝑑𝑡^′, (30) where 𝑊 is not a potential function. Its factorized version is then written as 𝝉=^𝑡_−∞𝑚𝑡−𝑡^′𝐻₁𝐼₂,𝐼₂𝐂⁻¹𝑡^′+𝐻₂𝐼₁,𝐼₂𝐂𝑡^′𝑑𝑡^′, (31) where 𝐻 is a strain-dependent function.

Tanner [9] argues that, strictly speaking, the above equations are not the K-BKZ equation, and he offers an example of very fast deformations to see the difference. Also, the above equations do not include a potential function 𝑈, as the original BKZ model. However, as Tanner [9] continues to conclude, “there seem to be few cases where any problems have arisen with a lack of potential, and many applications proceed using versions of (28), despite the possible dangers that might arise if fast deformations were to be investigated.”

At this point and starting in 1976, there were serious efforts to make the above models tractable to fit rheological experimental data and make them available for simulations. The impetus was given by Meissner’s IUPAC Working Party data on low-density polyethylene (LDPE) published in 1975 [26]. The data was extensive, and a comprehensive set of rheological data for an important industrial polymer at last became available to people in the field.

Wagner [27] began to describe the IUPAC-LDPE data by K-BKZ forms. Generally he used the equation in the form 𝝉=^𝑡_−∞𝑚𝑡−𝑡^′ℎ𝐼₁,𝐼₂𝐂⁻¹𝑡^′𝑑𝑡^′, (32) where ℎ was termed the “damping function,” which is a function of 𝐼₁ and 𝐼₂; these are the traces of 𝐂⁻¹ and 𝐂, respectively, and 𝑚 is the time-dependent or memory function.

The memory function can be found from small-strain data according to the linear viscoelastic theory; it is usually given as a sum of exponential functions involving the relaxation times 𝜆_𝑘 and relaxation moduli 𝑎_𝑘 for 𝑁 relaxation modes 𝑚𝑡−𝑡^′=^𝑁_𝑘=1𝑎_𝑘𝜆_𝑘exp−𝑡−𝑡^′𝜆_𝑘. (33) The number of modes depends on the polymer melt, but usually 4–8 modes are enough to capture the linear viscoelastic data [27], while for polymer solutions less modes (2-3) are enough [1].

The damping function is not a unique function and many similar functions have been proposed to match experimental data. Wagner [27] for general flows used ℎ=exp−𝑛𝐼₂−3, (34) where 𝑛 is a fitted parameter. Later, Wagner et al. [28] used a strain invariant 𝐼=𝛽𝐼₁+1+𝛽𝐼₂ (35) to reconcile shear and elongational data; here 𝛽 is a fitted parameter to be determined from elongational data. Thus, the Wagner model uses 𝑚𝑡−𝑡^′=^𝑁_𝑘=1𝑎_𝑘𝜆_𝑘exp−𝑡−𝑡^′𝜆_𝑘,ℎ=exp−𝑛𝛽𝐼₁+1−𝛽𝐼₂−3. (36)

Around 1980 a number of damping functions were available, so that Dan Joseph in his graduate course on Rheology referred to his list of these functions as a “Drugstore of Rheology” [29]. Papanastasiou et al. [30], having taken Joseph’s course on rheology at the University of Minnesota, proposed the following fractional damping function: 𝛼ℎ=(𝛼−3)+𝛽𝐼₁+1−𝛽𝐼₂. (37)

The resulting model has been known in the literature as the Papanastasiou-Scriven-Macosko or K-BKZ/PSM model, and it has been widely used in numerical simulations as will be detailed later. The constants 𝛼 and 𝛽 in the damping function can be obtained independently from shear and elongational data, respectively.

A more complete form of the model was proposed by Luo and Tanner [31] to account for a nonzero second normal stress difference and for a better fitting of elongational data. The K-BKZ/L-T model is written as 1𝝉=1−𝜃^𝑡𝑁_−∞𝑘=1𝑎_𝑘𝜆_𝑘exp−𝑡−𝑡^′𝜆_𝑘×𝛼𝛼−3+𝛽_𝑘𝐼₁+1−𝛽_𝑘𝐼₂×𝐂⁻¹𝑡^′+𝜃𝐂𝑡^′𝑑𝑡^′, (38) where 𝑁₂𝑁₁=𝜃.1−𝜃 (39) In the above, 𝑁₁=𝜏₁₁−𝜏₂₂ and 𝑁₂=𝜏₂₂−𝜏₃₃ are the first and second normal stress differences, respectively. The constant 𝜃 is a negative number and usually in the range of −0.1 to −0.25 [1]. The multiple 𝛽_𝑘 values correspond to different elongational behavior for each relaxation mode.

The above more tractable integral equations were further rectified. Wagner [32] introduced the “irreversibility” argument to solve the excessively elastic recoil problem exhibited by the above models. In this model, the damping function ℎ becomes “inelastic” or “irreversible” in the following way: when the strain changes between time 𝑡^′ and 𝑡, one selects the lowest value of ℎ that exists over this time segment. Another rectification was given by Olley [33], who proposed a modification of the PSM damping function so that different responses could be obtained in uniaxial and in planar extension (note that with the PSM damping function the planar extensional viscosity follows the shear viscosity). The Olley damping function reads 𝛼ℎ=𝛼+𝛽𝐼₁−3+1−𝛽𝐼₁−3 (40) and can give independent planar extensional behavior from its shear counterpart without compromising the uniaxial elongational behavior.

Tanner [9] goes on to relate how the K-BKZ model is also connected with microstructural theories. Under this line of thought the reptation model of Doi and Edwards [34–36] for entangled linear polymers, further extended by Curtiss and Bird [37], can also be written in the form of (28) using Currie’s [38] approximation of the related potential function. Then the Doi-Edwards model is given by 𝜙₁=5𝐼₁+2𝐼₂,𝜙+3.25−1₂𝜙=−¹𝐼₂,𝑚+3.25𝑡−𝑡^′=96𝜇𝜋⁴𝜆^2∞_𝑘=0exp−(2𝑘+1)²𝑡−𝑡^′𝜆. (41) It contains only two adjustable parameters, namely, the zero-shear-rate viscosity 𝜇 and the time constant 𝜆.

Integral constitutive equations proposed more recently on the basis of molecular theories for linear and branched polymers have a more complicated mathematical structure (integrodifferential forms) and are beyond the scope of the present paper. The interested reader is directed towards the developments made by Mead et al. [39], McLeish and Larson [40], Ianniruberto and Marrucci [41], Wagner et al. [42], Hassager et al. [43], Hassager and Hansen [44], and so forth.

3. Fitting Experimental Data with the K-BKZ Model

The idea behind any model is to fit well experimental data in an effort to use physics and mathematics to mimic nature. The K-BKZ model and its more tractable variants (such as the PSM model [30] or the L-T model [31]) have enough degree of complexity as to guarantee good fit of many experimental data obtained from rheometry. In the following, examples are given of how to use the model to obtain useful fitting of experimental data for two main polymer melts of industrial importance, namely, a low-density polyethylene (LDPE) and a high-density polyethylene (HDPE).

3.1. Data Fitting

In the review article by Tanner [9] celebrating 25 years of the K-BKZ model, an example was given of how to fit the rheological data for the IUPAC-LDPE melt. It was shown that the time-memory function 𝑚(𝑡−𝑡^′) consisted of 𝑁=8 relaxation modes, and for each mode there were values of 𝑎_𝑘, 𝜆_𝑘, and 𝛽_𝑘, 𝑘=1,8. The fitting was completed by setting the other two constants, namely, 𝛼 and 𝜃. It was thus necessary to have 24+2=26 constants to fit the rheological data.

One may argue that these are just too many constants to fit. However, and upon closer inspection, this is not such a daunting task. The required relaxation spectrum (𝑁,𝑎_𝑘,𝜆_𝑘) can be easily obtained from a single experiment of collecting data for the linear viscoelastic storage and loss moduli, 𝐺^′ and 𝐺^′′, respectively. These in turn are given as a function of frequency 𝜔 as follows (by assuming a Maxwell model for each mode) [1]: 𝐺^′(𝜔)=^𝑁_𝑘=1𝑎_𝑘𝜔𝜆_𝑘²1+𝜔𝜆_𝑘²,𝐺^′′(𝜔)=^𝑁_𝑘=1𝑎_𝑘𝜔𝜆_𝑘1+𝜔𝜆_𝑘². (42) These functions are independent of the strain-memory function (hence separability of the functions or factorized K-KBZ model), and only 𝑎_𝑘 and 𝜆_𝑘 can be determined from dynamic data of the viscoelastic moduli. As an example we give in Figure 4 and Table 1 the fitting of experimental 𝐺^′and 𝐺^′′ data for an LDPE melt [45] and an HDPE melt [46], where it is shown that a spectrum of 8 relaxation times, ranging from 10⁻⁴ s to 10⁺³ s, is able to capture well the data at all frequencies. It should be noted that the fitting is not very difficult and can be done even with an MS Excel software, provided that one uses a reasonable number of relaxation times (8 for polymer melts or 3 for polymer solutions) with initial guesses for the relaxation times submultiples or multiples of 10 and appropriate initial guesses for the relaxation moduli [30, 47]. A commercial code IRIS [48] is also available based on the parsimonious model of Baumgaertel and Winter [49], which finds the minimum number 𝑁 of relaxation modes to fit the data. Thus, the time-memory function 𝑚(𝑡−𝑡^′) presents no problem in fitting well experimental data for 𝐺^′ and 𝐺^′′ and finding the related parameters.

Table 1 

Relaxation spectra and material constants for two polymer melts obeying the K-BKZ model: LDPE at 160°C [45] and HDPE at 190°C [46].

(a)

(b)

Figure 4 

Experimental data (symbols) and model predictions (lines) of storage (𝐺′) and loss (𝐺′′) moduli for two polymer melts: (a) LDPE melt at 160°C [45] and (b) HDPE at 190°C [46]. The relaxation spectrum is listed in Table 1.

The strain-memory or damping function ℎ(𝐼₁,𝐼₂) is derived from the first and second invariants of the Finger strain tensor. For simple shear flow, the strain-memory function is given as ℎ𝐼₁,𝐼₂=𝛼𝛼+𝛾², (43) where 𝛾 is the shear strain. The strain-memory function in simple shear flow is dependent on 𝛼 but not on 𝛽. This is expected since 𝛼 is viewed as a shear parameter, while 𝛽 is viewed as an elongational parameter. Thus, data on ℎ versus 𝛾 can be fitted to give the value of parameter 𝛼. However, an easier way is to have shear viscosity (𝜂_𝑆) data as a function of shear rate, something which is easily available with today’s rheometers. Again, the shear viscosity data is used to find the best value of parameter 𝛼. This is shown in Figure 5 for the two melts, namely, LDPE and HDPE.

(a)

(b)

Figure 5 

Experimental data (symbols) and model predictions (lines) of shear viscosity, 𝜂𝑆, first normal stress difference, 𝑁1, and elongational viscosity, 𝜂𝐸, for two polymer melts: (a) LDPE melt at 160°C [45] and (b) HDPE melt at 190°C [46] using the K-BKZ/PSM model with the parameters listed in Table 1.

The elongational parameter 𝛽 (or 𝛽_𝑘 if multiple parameters are needed) is determined from elongational data, either steady state or transient. In the earlier days of the 1960s and 1970s, elongational viscosity data were very hard to come by, and Meissner and his group [26] were virtually the only ones to produce reliable data on elongational viscosity. This situation has now been much improved because of two important developments, the filament-stretching device of Sridhar and McKinley [50, 51], and the Sentmanat Elongational Rheometer (SER) [52]. Thus, fitting elongational data obtained by these devices, the parameter 𝛽 of the model is readily obtained, as shown in Figure 6 for the LDPE and HDPE melts. The value of 𝜃 is usually set to be a small negative number (around −0.1) according to experimental evidence [31].

(a)

(b)

Figure 6 

Experimental data (symbols) and model predictions (lines) of transient elongational viscosity, 𝜂𝐸, for two polymer melts: (a) LDPE melt at 160°C [45] and (b) HDPE melt at 190°C [46] using the K-BKZ/PSM model with the parameters listed in Table 1.

With the determination of the fitting parameters, it is thus possible to obtain a good fit for steady-state rheometric functions, such as the shear viscosity, 𝜂_𝑆, first normal stress difference, 𝑁₁, and uniaxial elongational viscosity, 𝜂_𝐸. Then the other extensional viscosities in planar extension, 𝜂_𝑃, and in biaxial extension, 𝜂_𝐵, are predicted by the model. These predictions can also be extended to transient effects for all rheological functions at different times [30, 47].

4. Mathematical Modeling

4.1. Governing Conservation Equations

In order to study polymer flows in rheology and inside processing equipment, it is essential to consider first the governing flow equations. The flow of incompressible fluids (such as polymer solutions and melts, at least in situations where they are considered as incompressible for pressures below 100 MPa) is governed by the usual conservation equations of mass, momentum, and energy [1, 23], that is, ∇⋅𝐯=0,𝜌𝐯⋅∇𝐯=−∇𝑝+∇⋅𝝉,𝜌𝑐_𝑝𝐯⋅∇𝑇=𝑘∇²𝑇+𝝉∶∇𝐯, (44)(45)(46) where 𝐯 is the velocity vector, 𝑝 is the scalar pressure, 𝝉 is the extra stress tensor, 𝜌 is the density, 𝑐_𝑝 is the heat capacity, 𝑘 is the thermal conductivity, and 𝑇 is the temperature.

The above system of conservation equations is usually called the nonisothermal Cauchy or momentum equations in Fluid Mechanics.

4.2. Constitutive Equations

The above system of conservation equations is not closed for non-Newtonian fluids, such as polymeric liquids, due to the presence of the stress tensor 𝝉. The required relationship between the stress tensor 𝝉 and the kinematics (velocities and velocity gradients) must be given by appropriate rheological constitutive equations, and this is an eminent subject in theoretical rheology as explained in Section 1.1 [1, 23, 24]. To address this issue all the aforementioned efforts in developing the K-BKZ model were directed.

However, it is always instructive to consider the simplest possible models for purely viscous fluids, so that the differences due to viscoelasticity with the K-BKZ model become apparent. The simplest viscous rheological constitutive equation that relates the stresses 𝝉 to the velocity gradients is the generalized Newtonian model [1, 23] and is written as 𝝉=𝜂𝛾̇𝜸̇, (47) where 𝜸̇=∇𝐯+∇𝐯^𝑇 is the rate-of-strain tensor and 𝜂(𝛾̇) is the apparent viscosity given in its simplest form by the power-law model [1, 23] 𝜂𝛾̇=𝐾𝛾̇^𝑛−1, (48) where 𝐾 is the consistency index and 𝑛 is the power-law index (usually 0<𝑛<1, representing a degree of shear thinning). Other popular models for viscosity computations—among others—are the Carreau model [23] given by 𝜂𝛾̇=𝜂_∞+𝜂₀1+𝜆_𝐶𝛾̇^{2(𝑛−1)/2} (49) and the Cross model [53] given by 𝜂𝛾̇=𝜂_∞+𝜂₀1+𝜆_𝐶𝛾̇^1−𝑛. (50) In the above, 𝜂₀ is the zero-shear-rate viscosity, 𝜂_∞ is the infinite-shear-rate viscosity, 𝜆_𝐶 is a characteristic time, and 𝑛 is the power-law index. The magnitude 𝛾̇ of the rate-of-strain tensor is given by 𝛾̇=12II_̇𝛾=12𝛾̇∶𝛾̇, (51) where II_̇𝛾 is the second invariant of the rate-of-strain tensor. The Carreau model describes well the shear thinning behaviour of polymer solutions and melts for all shear rates and exhibits two plateaus for low and for high shear rates, while for intermediate to large shear rates it represents well the power law.

The effect of temperature on the viscosity is of primordial importance in polymer processing flows, where tight control of temperatures is required for a successful operation. The viscosity as a function of temperature is given by an exponential relationship, according to [1, 23, 53] 𝜂_𝑇=𝜂₀exp−𝛽_𝑇𝑇−𝑇₀, (52) where 𝛽_𝑇 is a temperature-shift factor in the expression that relates viscosity to temperature, and 𝜂₀ is the viscosity at a reference temperature, 𝑇₀. The values of 𝛽_𝑇 for polymers are usually in the range of 0.01–0.04/C but occasionally they may reach 0.1/C or more for some polymers.

Another expression for the temperature dependence of the viscosity is the Arrhenius law [1, 23, 53]: 𝜂_𝑇=𝜂₀𝐸exp^𝑎𝑅_𝑔1𝑇−1𝑇₀, (53) where 𝑅_𝑔 is the ideal-gas constant (=8.13 J/K⋅mol), 𝐸_𝑎 is the activation energy for flow (J/mol), 𝑇 is the absolute temperature (K), and 𝑇₀ is the absolute reference temperature (K).

Then combining (52) and (53) yields 𝛽_𝑇=𝐸_𝛼𝑅_𝑔𝑇𝑇₀. (54) It is convenient to define the Arrhenius temperature-shifting function 𝑎_𝑇 as 𝑎_𝑇(𝑇)=𝜂_𝑇𝜂₀𝐸=exp𝑅_𝑔1𝑇−1𝑇₀. (55)

For the integral K-BKZ model and its variants, the nonisothermal version comes about from applying the time-temperature superposition principle [1]. This simply consists of shifting the relaxation times 𝜆_𝑘 from the temperature history within the material’s internal time scale 𝑡^′ [54]. The equation used to shift the relaxation times in the material’s history is given by [54–56] 𝜆_𝑘𝑇^′𝑡^′=𝜆_𝑘𝑇₀𝑎_𝑇𝑇^′𝑡^′, (56) where 𝑇^′ is the temperature at time 𝑡^′. The viscoelastic stresses calculated by the nonisothermal version of the K-BKZ model enter in the energy equation (see (46)) as a contribution to the viscous dissipation term.

Similarly with the time-temperature superposition principle where the stresses are calculated at a different temperature using the shift factor 𝑎_𝑇, the time-pressure superposition principle can be used to account for the pressure effect on the viscosity and hence the stresses. In both cases of viscous or viscoelastic models, the new stresses are calculated using the pressure-shift factor 𝑎_𝑝 defined by the Barus equation [57, 58] 𝑎_𝑝≡𝜂_𝑝𝜂_𝑝0=exp𝛽_𝑝𝑝, (57) where 𝜂_𝑝 is the viscosity at absolute pressure 𝑝, 𝜂_𝑝0 is the viscosity at ambient pressure, and 𝛽_𝑝 is the pressure coefficient.

For viscous models, (57) is used to modify the viscosity. For viscoelastic models, such as the K-BKZ model, the pressure-shift factor modifies the relaxation moduli, 𝑎_𝑘, according to 𝑎_𝑘𝑝𝑡^′=𝑎_𝑘𝑝₀𝑎_𝑝𝑝𝑡^′. (58) This is equivalent to multiplying the viscoelastic stresses by 𝑎_𝑝.

4.3. Dimensionless Groups

Before proceeding with the boundary conditions, it is important to examine the relevant dimensionless numbers in polymer flows through dies used in rheometry and in processing equipment. The dimensionless groups are calculated at a reference temperature, here taken as the temperature of the process, 𝑇₀. As a characteristic length it is usually assumed the smallest dimension, for example, in a capillary tube its radius, 𝑅. As a characteristic speed it is usually assumed the average velocity defined by using the volumetric flow rate 𝑄 according to 𝑄𝑉=𝜋𝑅². (59)

A characteristic apparent shear rate is then defined according to 𝛾̇_𝐴=𝑉𝑅 (60) and a characteristic viscosity is given as a function of apparent shear rate and reference temperature, that is, 𝜂=𝜂𝛾̇_𝐴,𝑇₀. (61) For a power-law model, the characteristic viscosity can be found as 𝜂𝑇₀=𝜏𝛾̇_𝐴=𝐾𝛾̇_𝐴^𝑛−1, (62) where the material parameters 𝐾 and 𝑛 are calculated at 𝑇₀.

The relative importance of inertia forces in the equation of momentum is assessed by the Reynolds number, defined for Newtonian fluids by Re=𝜌𝑉𝐷𝜇, (63) where 𝐷 is the characteristic diameter (=2𝑅). It is noted that for most polymer melt flows the Re number is usually small, in the range of 0.0001 to 0.01. Therefore these flows are inertialess or creeping.

For viscoelastic fluids with a relaxation time 𝜆, several dimensionless groups can be defined, but these can be seen as being equivalent [59]. For example, the Deborah number (De) is defined as 𝜆De=𝑡_𝑝=𝜆𝛾̇, (64) where 𝜆 is a material relaxation time, 𝑡_𝑝 is a process relaxation time usually taken to be equal to 1/𝛾̇, and 𝛾̇ is a shear rate usually evaluated at the channel wall. The Weissenberg number (Wi, also written as We or Ws) is defined as 𝑉Wi=𝜆𝑅. (65) The recoverable shear or stress ratio (𝑆_𝑅) is defined as 𝑆_𝑅=𝑁_1,𝑤2𝜏_𝑤, (66) where 𝑁_1,𝑤=(𝜏₁₁−𝜏₂₂)_𝑤 is the first normal stress difference and 𝜏_𝑤 is the shear stress, both evaluated at the channel wall. The equivalence is evident when we take ̇𝑉𝛾=𝑅,𝑁_1,𝑤=Ψ₁𝛾̇²,𝜏_𝑤̇Ψ=𝜂𝛾,𝜆=¹,2𝜂 (67) where Ψ₁ is the first normal stress difference coefficient and 𝜂 is the shear viscosity. The case of  De=Wi=𝑆_𝑅=0 corresponds to inelastic fluids (𝜆=Ψ₁=0), while it is understood that  De=Wi=𝑆_𝑅=1 corresponds to the elastic effects being as important as the viscous effects, and for De=Wi=𝑆_𝑅>1 the elastic effects dominate the flow over the viscous effects.

The relative importance of surface tension effects (usually for polymer solutions) is assessed by the capillary number, defined by Ca=𝜇𝑉𝛾, (68) where 𝛾 is the surface tension. For very viscous fluids, such as polymer melts, the surface tension effects are negligible (Ca→∞), and the boundary terms, including a force balance with the capillary forces, can be set to 0.

The relative importance of each term in the energy equation is assessed through a variety of dimensionless groups [60].

The Peclet number is defined by Pe=𝜌𝐶_𝑝𝑉𝐷𝑘. (69)

The Peclet number is a measure of convective heat transfer with regard to conductive heat transfer. High Pe values indicate a flow dominated by convection. Another group related to Pe is the Graetz number defined by Gz=𝜌𝐶_𝑝𝑉𝐷²𝐷𝑘𝐿=Pe𝐿, (70) where 𝐿 is the axial length of the die. The Graetz number can be understood as the ratio of the time required for heat conduction from the center of the capillary to the wall and the average residence time in the capillary. As with Pe, a large value of Gz means that heat convection in the flow direction is more important than conduction toward the walls.

The Nahme number is defined by 𝛽Na=^𝑇𝜂𝑉²𝑘. (71)

The Nahme number is a measure of viscous dissipation effects compared with conduction; hence it is an indicator of coupling of the energy and momentum equations. For values of Na > 0.1–0.5 (depending on geometry and thermal boundary conditions), the viscous dissipation leads to considerable coupling of the conservation equations, and a nonisothermal analysis is necessary.

The relative importance of heat transfer mode at the boundaries is expressed in terms of the dimensionless Biot number defined by Biot=𝜕𝑇𝜕𝑟_𝑤𝑟_𝑖𝑇_𝑠−𝑇_𝑤, (72) where 𝑟_𝑖 is the local radius (gap), 𝑇_𝑠 is some temperature of the surroundings, and 𝑇_𝑤 is the local boundary wall temperature. A high value of Biot (Biot>100) approaches isothermal conditions (Biot=∞), while a low value of Biot (Biot<1) describes poor heat transfer to the surroundings (nearly adiabatic case, for which Biot=0). Usual Biot values in highly viscous flows inside dies or other processing equipment range between l0 and 100 [60].

For engineering calculations, the specific heat flux 𝐪 is often described by the Nusselt number ℎNu=^𝑇𝐷𝑘, (73) where ℎ_𝑇 is a heat transfer coefficient. However, as explained by Winter [60], the Nusselt number is not adequate for describing the wall heat flux in flows with significant viscous dissipation.

4.4. Boundary Conditions

The solution of the conservation equations (44)–(46) and constitutive equation (47) (or whatever viscoelastic K-BKZ model is used for the stresses) is only possible after a set of boundary conditions (BCs) has been imposed on the flow domain. Boundary conditions for flow analysis are highly dependent on the problem at hand, and as such they defy a complete description for all polymer processing applications. However, a rough guide that encompasses most of the types of boundary conditions used in the past follows.

For steady-state problems, the set of equations is elliptic for viscous flows and elliptic-hyperbolic for viscoelastic flows. Elliptic problems have boundary conditions everywhere in the perimeter of the domain, while hyperbolic problems are more difficult to determine their boundary conditions and may need some degree of trial and error. This is especially true at the outflow boundaries, which are more often than not “artificial” or “computational” boundaries, arbitrarily set to reduce the computational domain.

For unsteady-state problems, the set of equations is parabolic due to time, and it also requires initial conditions at time 𝑡=0. Then the solution proceeds in time until a steady state is reached (if such exists).

The boundary conditions can be either “fixed” or “natural”. These types refer, respectively, to the primary variables (velocities, pressures, and temperatures) or variables involving their derivatives (stresses, surface tractions, heat fluxes, etc.).

The usual flow boundary conditions in polymer processing are as follows:(a)along the domain entry, a fully developed velocity profile is imposed, according to the assumed constitutive model and corresponding to the apparent shear rate 𝛾̇_𝑎, which is in turn related to the volumetric flow rate 𝑄. If viscoelasticity is involved through some viscoelastic constitutive model, then the fully developed stress profiles have to be imposed as well;(b)along the centerline (if one exists, as is the case in axisymmetric flows), and because of symmetry, the radial velocity component is set to zero, as well as the shear stresses;(c)along solid walls, usually the no-slip velocity boundary condition is imposed, which states that the velocity of the fluid is the same as that of the boundary, that is, zero if the boundary is stationary, or non-zero if the boundary is moving. In cases where the fluid slips at the wall, as is the case with some polymers [61], then a slip law has to be assumed based on measurements, which relates the tangential velocity to the tangential components of the stress tensor, while the normal velocity is set to 0;(d)along the domain exit, it is not clear what are the correct or physical boundary conditions for all situations. The best candidate appears to be the “open” or “free” boundary condition (FBC) advocated by Papanastasiou et al. [62], which basically assumes an extrapolation of the governing equations to the artificial exit boundary. However, the vast majority of the computations assume a long enough domain, where they impose zero surface tractions and zero transverse velocity (assuming implicitly a fully developed profile at exit). For viscoelastic models, fully developed stresses based on the model at hand are used, unless there is a known force or velocity imposed at the exit boundary;(e)in problems with free surfaces, and especially for polymer melts due to their very viscous character (zero surface tension assumed), zero surface tractions are imposed along with a kinematic boundary condition of no flow normal to the surface, that is, 𝐧⋅𝐯=0, where 𝐧 is the unit outward normal vector to the surface.

For the thermal boundary conditions, the situation is even more difficult to describe due to the many types of thermal conditions used in polymer processing. However, in the vast majority of computations, the following set of boundary conditions is a good representation of what has been applied, based on the previous flow sets:(a)along the domain entry, a set (quite often constant) temperature profile is assumed, sometimes based on measurements or settings;(b)along the centerline (if one exists, as is the case in axisymmetric flows), and because of symmetry, the heat flux is set equal to zero;(c)along solid walls, the surfaces are assumed either as isothermal walls (a set temperature) or as adiabatic walls (heat flux set to 0) or a heat balance between the fluid and the solid boundary [60, 63]. In the latter case, a local Biot number can then be calculated, which is neither 0 (adiabatic walls) nor ∞ (isothermal walls), but may range between 10 and 100. Other types of wall thermal boundary conditions involve a known heat flux at the wall based on measurements of an effective heat transfer coefficient, but this is also tantamount to having a non-zero local Biot number;(d)along the domain exit, the same comments apply as above, regarding the flow BCs, with the best candidate being the “open” or “free” BC. However, the majority of computational problems have assumed a long enough domain so as to impose a zero heat flux at the exit;(e)in problems with free surfaces, a zero heat flux is usually imposed, that is, 𝐧⋅𝐪=0, or a heat balance according to 𝑞_𝑛=ℎ_𝑇𝑇_𝑓−𝑇_𝑎, (74) where ℎ_𝑇 is a convective heat transfer coefficient to the ambient cooling medium (e.g. air) of temperature 𝑇_𝑎, and 𝑇_𝑓 is the unknown free surface temperature.

For unsteady-state processes, initial conditions are also needed.

5. Method of Solution

5.1. A Brief Literature Survey

It is clear from the start that integral constitutive equations of the K-BKZ type are not trivial. On the contrary, they constitute a major numerical challenge as they have to be solved together with the conservation equations for real polymer flows. Obviously, such a daunting task has kept the computational rheology community on its feet for the better part of the last 30 or so years, starting in 1980 with the pioneering work of Viriyayuthakorn and Caswell [64], who devised the first finite element technique for integral viscoelastic models. Since then many groups around the world have tackled the issue and many dedicated people have contributed towards a much better state of affairs 30 years later. This was also made possible by the well-known advances in computational power offered by digital computers.

The review article by Keunings [65] captures in an eloquent and precise way the developments made with the Finite Element Method (FEM) for integral constitutive equations up to 2003. Here we follow Keunings [65] and pay attention to some interesting developments associated with real polymer flows.

After the work of Viriyayuthakorn and Caswell [64] who solved the integral Maxwell model, Papanastasiou et al. [66] devised a coupled streamline FEM (CS-FEM) to solve the K-BKZ/PSM model. In CS-FEM the governing equations are solved simultaneously for the whole set of primary variables with the Newton-Raphson iterative technique. The memory integral is evaluated along particle paths (pathlines, see Figure 3), that are unknown a priori. Problems with storage requirements were encountered and this elegant method was not pursued any further.

A decoupled Eulerian-Lagrangian approach was then pursued, wherein the computation of the viscoelastic stresses is performed separately from that of the flow kinematics. This requires a Picard (or direct substitution) iterative scheme, which is slow to converge as it has linear convergence characteristics. Dupont et al. [67] used this method to solve the integral Maxwell fluid and the Doi-Edwards model, while later Dupont and Crochet [68] solved with the same method the K-BKZ/PSM model for an LDPE melt. Meanwhile, Luo and Tanner [69] had developed a streamline element FEM (SE-FEM) to compute the stresses along unknown a priori streamlines, and with this method they solved the integral Maxwell model and the K-BKZ/PSM model for flows with open streamlines (no recirculation present). This work was considered a breakthrough because for the first time the authors were able to achieve convergence for high enough flow rates, as those encountered in real polymer melt flows. This was attributed to the realistic behavior of the K-BKZ model and the accuracy of their integration scheme. Luo and Tanner [31] went further to present and solve the K-BKZ/L-T model in extrusion flows of the IUPAC-LDPE melt.

The SE-FEM was modified by Luo and Mitsoulis [70] to account for both open and closed streamlines (recirculation present) by removing the requirement to construct streamlines in the flow field. Luo and Mitsoulis [71] used this method to also solve the now considered as a benchmark problem for integral constitutive equations, namely, the flow of the IUPAC-LDPE melt in extrusion dies with the integral K-BKZ/PSM model.

The same problem was also solved by Bernstein et al. [72] and Feigl and Ottinger [73] using similar integrating schemes for the viscoelastic stresses. At this point it is to be noted that at least 4 different groups around the world (Crochet’s in Belgium, Tanner’s in Australia, Mitsoulis’s in Canada, and Bernstein’s in the USA) employing different codes were able to solve the same benchmark problem of die flow of the IUPAC-LDPE melt with the K-BKZ/PSM or its variants, obtaining similar (if not identical) results. In this way a confidence was built in the different schemes for the integration of the viscoelastic stresses according to K-BKZ-type models.

At the same time, a Lagrangian Integral Method (LIM) was advocated by Hassager and Bisgaard [74]. As its name indicates, LIM is based on the Lagrangian formulation of the conservation laws. The primary kinematic variables are the time-dependent positions of the fluid elements, while the computational domain is a deforming material volume. The latter is discretized by means of a finite element mesh that deforms with the flow. A typical flow simulation with LIM amounts to solving a transient problem, given the initial particle positions and deformation pre history. LIM is thus a coupled method capable of solving transient flow problems. Steady-state flows are computed as the long-time limit of a transient problem. As all Lagrangian techniques, LIM readily applies to flows with free surfaces. Over the years, the method has been refined and extended to treat three-dimensional flows (Rasmussen and Hassager [75, 76], Rasmussen [77, 78]). To this date, it remains one of the few cases (if not the only one) where integral constitutive equations are solved for transient 3D problems.

In the review article by Keunings [65] there are also references to the method of deformation fields by Peters et al. [79] and to methods applied for integrodifferential models derived from molecular theory (namely, the Doi-Edwards models). However, these methods are beyond the scope of the present paper, and interested readers are kindly asked to consult the excellent review by Keunings [65].

5.2. The Decoupled Eulerian-Lagrangian Approach

In what follows we present in some detail the numerical technique implemented with Picard iterative schemes, namely the decoupled Eulerian-Lagrangian approach, which has been most widely used in numerical simulations of polymer flows with the variants of the K-BKZ model, namely, the PSM and the L-T models.

In the decoupled techniques the computation of the viscoelastic stresses is performed separately from that of the flow kinematics. These techniques are based on an Eulerian-Lagrangian approach, which combines the computation of Eulerian velocity and pressure fields with the Lagrangian evaluation of the strain history along particle trajectories. The decoupled Eulerian-Lagrangian techniques are aimed at steady-state flows only. They have been under development since the end of the 1970s, and most of the fundamental ideas have been reviewed in detail by Crochet et al. [80] and Keunings [81].

5.2.1. Iterative Strategy

The iterative strategy consists of the following steps:(1)using the velocity field computed at the previous iteration, integrate the constitutive equation to update the viscoelastic stresses;(2)using the viscoelastic stress computed in step 1, update the kinematics by solving the conservation laws;(3)check for convergence; if needed, return to step 1 for another iteration.

It was realized early on [82] that the viscoelastic stresses can be split into a viscous (Newtonian) component and a purely viscoelastic component by adding and subtracting the viscous component in the momentum equation. This scheme has been called the elastic-viscous split stress (EVSS) scheme and it is widely used in viscoelastic simulations [80, 81]. Effectively, the extra stress 𝝉 in the momentum equation is replaced by 𝝉̃(𝜔)=𝜔𝝉+(1−𝜔)𝜂_ref𝜸̇, (75) where 𝜂_ref is an arbitrary Newtonian viscosity, and 𝜔 is a control parameter between 0 and 1 (under-relaxation parameter). For 𝜔=0, the flow problem is purely Newtonian, while the original viscoelastic equation is recovered for 𝜔=1.

Let 𝝉^𝑛 denote the viscoelastic stress computed at the 𝑛th iteration. Then, for creeping flows, step2 amounts to solving ∇⋅−𝑝^𝑛+1𝜹+𝜂_𝑟𝑒𝑓𝜸̇^𝑛+1=−∇⋅𝜔𝝉^𝑛+1−𝜂_ref𝜸̇^𝑛,∇⋅𝐯^𝑛+1=0 (76) for the updated velocity and pressure fields 𝐯^𝑛+1 and 𝑝^𝑛+1, respectively. During the iterative process, the control parameter 𝜔 is progressively increased from 0 to 1, so that a solution to the original viscoelastic equations is finally obtained. This strategy is similar to the incremental loading procedure used in non-linear elasticity, where now the “load” is the viscoelastic stresses. In all cases, the solution starts from the Newtonian solution (𝜔=0) and proceeds from there.

The above strategy was able to reproduce the exact analytical solutions for the integral Maxwell fluid (IMX) in Poiseuille flow for values of Wi up to 2, but not higher. This is called the High-Weissenberg Number Problem (HWNP) [83] that had plagued viscoelastic simulations for almost all of the 1980s. A major breakthrough was achieved by the adaptive-viscous split stress (AVSS) scheme proposed from Tanner’s group [84], whereby the reference viscosity (𝜂_ref)_𝑒 at element level is given by (𝜂_ref)_𝑒=ℎ_𝑒𝝉^𝑒𝑙_max|𝐕|_max, (77) where ℎ_𝑒 is the element size, |𝝉^𝑒𝑙|_max is the maximum elastic stress value in the element, and |𝐕|_max is the maximum velocity value in the element. With this definition of changing 𝜂_ref per element, no upper limit of convergence was detected for the benchmark problem of IMX Poiseuille flow, and solutions were attainable at much higher Wi values for other problems as well [84].

5.2.2. Finite Element Equations

Step2 of the decoupled strategy is a Newtonian flow problem defined by (76), where −∇⋅(𝜔[𝝉^𝑛+1−𝜂_ref𝜸̇^𝑛]) is treated as a pseudo-body force that is known from the previous iteration. The corresponding Galerkin finite element equations read _Ω−𝑝^𝑎𝜹+𝜂_ref𝜸̇^𝑎⋅∇𝜙^𝑗𝑑Ω=−_Ω𝜔𝝉−𝜂_ref𝜸̇^𝑎⋅∇𝜙^𝑗𝑑Ω+_Γ𝐓𝜙^𝑗𝑑Γ,_Ω∇⋅𝐯^𝑎𝜓^𝑙𝑑Ω=0. (78) In the above, Ω is the flow domain, and Γ is part of the boundary with outward unit normal vector 𝐧 where the contact force (surface traction) 𝐓=𝝈⋅𝐧 is specified; the discrete velocity and pressure fields are approximated by 𝜐^𝑎=^𝑁𝜐_𝑖=1𝜐^𝑖𝜙^𝑖,𝑝^𝑎=^𝑁𝑝_𝑘=1𝑝^𝑘𝜓^𝑘, (79) where 𝜐^𝑖 and 𝑝^𝑘 are unknown nodal values, 𝜙^𝑖 and 𝜓^𝑘 are given shape functions, and 𝑁_𝜐 and 𝑁_𝑝 are the number of velocity and pressure nodes, respectively. In (78), 1≤𝑗≤𝑁_𝜐, 1≤𝑙≤𝑁_𝑝.

The pseudo-body force integral in the right-hand side of the Galerkin momentum equations is computed using the viscoelastic stress and velocity fields known from the previous iteration. To this end, we must somehow evaluate the viscoelastic stress 𝝉^𝑛 at all integration points of the finite element mesh; this defines Step 1 of the decoupled strategy, which we discuss next.

5.2.3. Integration of the Constitutive Model

Computation of the viscoelastic stress 𝝉^𝑛 at each integration point of the mesh is performed in three basic steps, using the Lagrangian formulation of the constitutive model.(1)Particle tracking: using the steady-state velocity field known from the previous iteration, compute the past (i.e. upstream) trajectory and the travel time of the integration point.(2)Strain evaluation: at selected past times, compute the deformation gradient tensor 𝐅 and from it the integrand of the constitutive model.(3)Stress evaluation: compute the memory integral numerically along the past trajectories, using the results of the first and second steps.

The last task is the simplest one. In a steady-state flow, the displacement function 𝐱^′ is independent of 𝑡 and may be written as 𝐱^′=𝐱^′(𝐱;𝑠), where 𝑠 is the time lapse 𝑡−𝑡^′(𝑠≥0). This implies that 𝐂⁻¹=𝐂(𝑠). With a single-exponential memory function (although multiple relaxation times are easily handled), the generic integral model 𝝉=^𝑡_−∞𝑚𝑡−𝑡^′𝑓₁𝐂⁻¹+𝑓₂𝐂𝑑𝑡^′ (80) becomes 𝑎𝝉=𝜆^∞₀exp−𝑠𝜆𝐒(𝑠)𝑑𝑠, (81) where 𝐒(𝑠)=𝑓₁(𝐂⁻¹)+𝑓₂(𝐂). Following the early work of Viriyayuthakorn and Caswell [64], a possible approach is to use a Gauss-Laguerre integration rule to approximate the memory integral: 𝑎𝝉≈𝜆^𝑁_𝑖=1𝑤^𝑁_𝑖𝐒𝜆𝑧^𝑁_𝑖, (82) where 𝑧^𝑁_𝑖 are the roots of the 𝑁th Laguerre polynomial, and 𝑤^𝑁_𝑖 are the weights of the quadrature rule. In practice, 𝑁 is of order 10. This approximation implies that the deformation 𝐒 need only be computed at the discrete time lapses 𝑠_𝑖=𝜆𝑧^𝑁_𝑖. Obviously, ever larger strains must be evaluated when the memory of the fluid increases. The Gauss-Laguerre approach has been extended to finite and infinite relaxation spectra by Malkus and Bernstein [85]. More flexible integration techniques, though also more expensive, have been introduced by Dupont et al. [67] and Luo and Tanner [69]. The time integral equation (81) is transformed into a line integral along the past trajectory; the latter is subdivided into segments through which the fluid particle travels in a time shorter than the relaxation time 𝜆, and a Gauss quadrature rule is used to compute the integral along each segment. The backward integration is stopped once the marginal contribution of a segment is less than some specified tolerance. Interestingly, small relaxation times require a finer segmentation in order to capture the recent deformation history.

In principle, tracking of past particle positions and computation of the strain history require the integration of the kinematic equations 𝐷𝐱𝐷𝑠^′(𝑠)=−𝜐𝐱^′(𝑠),𝐷𝐅(𝑠)𝐷𝑠=−∇𝜐^𝑇𝐱^′(𝑠)⋅𝐅(𝑠) (83)(84) backward along each particle trajectory, with the respective initial conditions 𝐱^′(0)=𝐱,𝐅(0)=𝜹. (85)

The particle tracking procedure necessary with the above equations presents a major difficulty in computations with the K-BKZ integral model. The early numerical schemes reviewed by Crochet et al. [80] had severe numerical problems with these two tasks. A first satisfactory solution was advanced by Bernstein et al. [86] for two-dimensional planar flows. Briefly, it consists in using a special low-order finite element made of four triangles to compute the velocity field. The velocity finite element approximation is such that both tracking and strain history calculations can be performed analytically within each triangle. The price to pay for this nice result is the relatively low accuracy of the computed velocity field, and also the fact that the special element has spurious pressure modes. Extension of this approach to axisymmetric flows is reported by Bernstein et al. [72].

Luo and Tanner [69] introduced another approach based on the concept of streamline finite elements. These are quadrilateral elements which have a pair of opposed sides that are aligned with streamlines. From the computation of the stream function 𝜓 using the current velocity field, which amounts to solving the Poisson equation, ∇²𝜓=(𝜕𝑢/𝜕𝑦−𝜕𝑣/𝜕𝑥), the mesh is updated at each iteration so that the nodes (and thus the integration points) always lie along selected streamlines. A parametric representation of past trajectories is thus readily available in each element. The authors then compute the deformation history by solving the kinematic equation (84) backward along each discrete streamline using a 4th-order Runge-Kutta method. This approach is both fast and accurate, but it cannot handle recirculation regions. This problem was circumvented by Luo and Mitsoulis [70], who did not use streamline elements but integrated the tracking problem (see (84)) by means of a second-order scheme, valid for both open and closed streamlines, which is independent of the type of element.

The tracking procedure developed by Goublomme et al. [87] is also independent of the type of element. Pathlines are computed on an element-by-element basis, by means of a 4th-order Runge-Kutta integration of (83) within the parent element. This eliminates most of the complicated geometrical problems of earlier tracking schemes. In order to compute the deformation history, the authors follow the work of Dupont et al. [67], who showed that, for a steady-state flow, the deformation gradient 𝐅 can be computed by means of the velocity vectors at Lagrangian times 𝑡 and 𝑡^′ and values of a scalar quantity 𝛽 which obeys a first-order differential equation involving the velocity field along the pathline. Similarly for the tracking itself, Goublomme et al. [87] integrate the equation for 𝛽 in the parent element. Very efficient implementations of this method on message-passing parallel computers have been proposed by Aggarwal et al. [88] and Henriksen and Keunings [89].

5.3. Some Issues Arising in Viscoelastic Computations

Implementation of a finite element formulation for the mass and momentum equations uses the primitive-variables approach, that is, velocities and pressure (in 2D flows this is called the 𝑢-𝑣-𝑝  formulation). For viscoelastic formulations, the stresses are also needed as discrete variables, and this leads to a mixed-variable approach (called 𝑢-𝑣-𝑝-𝜏 formulation or EVSS for elastic-viscous split stress). A better formulation requires also the rates of strain as extra variables, thus leading to an enhanced mixed-variable approach (called 𝑢-𝑣-𝑝-𝜏-𝑔 formulation or DEVSS-G for discontinuous elastic-viscous split stress-strain rate). Variations of these also exist. For the energy formulation, the temperature 𝑇 is the single primitive field variable. Special streamline-upwind/Petrov-Galerkin (SU/PG) techniques are used for highly convective flows (high Pe numbers; Mitsoulis et al. [63]). More details about the SU/PG scheme can be found in the landmark paper by Marchal and Crochet [90].

At this point it is perhaps instructive to give a feel of the complexity of solving flow problems in polymer processing with viscous or viscoelastic constitutive models, and in the latter case the differences between using differential or integral models to describe polymer melts. We take as an example a 2-D flow problem and the FEM employing a typical 9-node Lagrangian quadrilateral element, as shown in Figure 7. A viscous problem, based on the 𝑢-𝑣-𝑝 formulation, would require 22 dof (degrees of freedom) per element (9𝑢, 9𝑣, and 4𝑝 variables). If the flow problem is also coupled with the thermal problem, then we have 31 dof/elem (9𝑇 variables). A serendipity element is cheaper as it has only 8 nodes (it lacks the centroid node), thus giving 20(28) dof/elem (flow/+thermal problem).

Figure 7 

A typical Lagrangian 9-node quadrilateral element for the primitive-variables approach (𝑢-𝑣-𝑝-𝑇 formulation) [63].

Now for a viscoelastic problem, with one (1) relaxation mode, the stresses and the strain rates have to be added to the nodal unknowns. In 2-D flows, we distinguish between planar (𝜏_𝑥𝑥,𝜏_𝑦𝑦,𝜏_𝑥𝑦;𝑔_𝑥𝑥,𝑔_𝑦𝑦,𝑔_𝑥𝑦) and axisymmetric flows (𝜏_𝑟𝑟,𝜏_𝑧𝑧,𝜏_𝑟𝑧,𝜏_𝜃𝜃;𝑔_𝑟𝑟,𝑔_𝑧𝑧,𝑔_𝑟𝑧,𝑔_𝜃𝜃). These extra nodal unknowns are defined at the corner nodes (due to linear interpolation), thus adding 4×3=12 stresses and 4×3=12 strain rates, a total of 24 extra dof (planar); and 4×4=16 stresses and 4×4=16 strain rates, a total of 32 dof (axisymmetric). Thus, the total dof/element are 22+24=46(55) dof/elem (flow/+thermal planar problem) and 22+32=54(63) dof/elem (flow/+thermal axisymmetric problem). Note that these numbers are for a single relaxation model.

For the 8 relaxation modes describing the IUPAC-LDPE melt [9] or the polymer melts of Figures 4–6, the corresponding numbers become 214(223) dof/elem (flow/+thermal planar problem) and 278(287) dof/elem (flow/+thermal axisymmetric problem). These numbers are given collectively in Table 2 and are necessary for solving any problem with a differential viscoelastic model, for example, the popular “pom-pom” model [40].

Table 2 

Degrees of freedom (dof) analysis needed for differential versus integral constitutive models for the 2-D finite element of Figure 7.

It is then obvious that even for a sparse FEM mesh of a 1000 elements the total number of dof climbs to 𝑂(10⁺⁵). So it is not surprising that even in today’s computers there are only very few viscoelastic problems that have been solved with the full spectrum and differential viscoelastic models (like the “pom-pom”) even for simple 2-D flows.

The situation is different with integral constitutive models, such as the K-BKZ/PSM model or the K-BKZ/L-T model. In this case, the formulation remains 𝑢-𝑣-𝑝(-𝑇) with the same number of dof/element as in viscous flows. The stresses are then calculated a posteriori along streamlines according to the integral using a 15-point Gauss-Laguerre quadrature suitable for exponentially fading functions [69]. The summation used in the integration is very fast with digital computers, and it does not make much difference in computational time by using either 1 or 8 relaxation modes. This is a major advantage of integral constitutive equations! Most of the computational time is used for the solution of the 𝑢-𝑣-𝑝 problem, while the a posteriori computation of the stresses takes approximately an equal amount of CPU time (at least this is the case for the middle range of flow rates; at low flow rates the Newtonian-like behavior of polymer melts makes the computation of stresses extremely fast; at the other end of very high flow rates, the computation of stresses becomes most time consuming). Thus, it was possible to solve many 2-D flow problems by the year 2000 in the personal computers (PCs) of the day, as evidenced in a series of papers by Mitsoulis and his co-workers [91].

As an example we give here the case of flow through a contraction of the IUPAC-LDPE melt A, using the K-BKZ/L-T model with 8 relaxation modes and the data given by Tanner [9], and employing the computational method of streamline integration developed by Luo and Mitsoulis [70]. In the year 2000, working on a PC Pentium II at 300 MHz with 128 MB RAM, for a mesh with 640 elements it was taking 20 CPU seconds for the frontal solution of the 𝑢-𝑣-𝑝 system of equations and another 20 CPU seconds for the calculation of stresses for each iteration. To get a solution at an elevated flow rate (corresponding to an apparent shear rate 𝛾̇_𝐴=100 s⁻¹), about 1000 iterations were needed for a total of 11 CPU hours. In today’s computers (Intel Core2 Duo at 2×2.66 GHz with 2 GB RAM), the same problem can be solved in the same overall time but with 10 times the number of finite elements (6400) and the results are the same, just smoother. Thus, it has been possible to solve viscoelastic problems of multimode fluids, such as polymer melts, in more complex geometries, such as those encountered in polymer processing, as will be discussed below.

This K-BKZ/PSM or L-T integral model has been used in numerical flow simulations for a number of flow problems more or less successfully (see, e.g., Luo and Tanner [31], Dupont and Crochet [68], Bernstein et al. [72], Barakos and Mitsoulis [92], Sun et al. [84], and Ansari et al. [45, 46]). The review by Keunings [65] on the subject gives a list of problems solved with this model through numerical simulation, including many flows from polymer processing operations. Other flows solved with a number of different constitutive equations can be found in an excellent book on computational rheology [93].

6. Some Applications of the Model

6.1. Flow Patterns in Contraction Flows

The flow of polymer melts through abrupt contractions has been a benchmark problem in computational rheology ever since visualization data on polyethylenes appeared in a capillary rheometer more than 50 years ago from Bagley and Birks [94] in 1960. The differences in behavior between an LDPE and a HDPE melt regarding vortex development, as shown in Figure 8, were a big challenge for the theoreticians and the numerical analysts to reproduce by an appropriate rheological model this difference in behavior, which is obviously attributed to viscoelasticity. It suffices to say that any viscous model fails miserably to capture this behavior, in particular the big vortices of the LDPE melt.

(a)

(b)

Figure 8 

Differences in flow behavior of polyethylene melts in a 20 : 1 axisymmetric abrupt contraction at the same capillary wall shear stress of 𝜏𝑤=0.1 MPa (Bagley and Birks [94]): (a) branched LDPE, (b) linear HDPE. The flow is from top to bottom.

With the K-BKZ/PSM model and the method developed by Luo and Mitsoulis [70] this problem was successfully solved, see Mitsoulis [95]. For the same shear stress level at the wall (𝜏_𝑤=0.1 MPa), the simulations showed a remarkable resemblance to the experimental patterns of the two melts as shown in Figure 9. The data for the LDPE are those for the IUPAC-LDPE melt [9]. Those of HDPE are given in [95]. This is a very good indication of the ability of the K-BKZ model to predict unusual non-Newtonian flow behavior for two major industrial polymers. After all, this was the original idea behind developing such models.

(a)

(b)

Figure 9 

Differences in flow behavior of polyethylene melts in a 20 : 1 axisymmetric abrupt contraction at the same capillary wall shear stress of 𝜏𝑤=0.1 MPa (simulation results with the K-BKZ/PSM model): (a) branched LDPE at 𝛾̇𝐴=84 s−1, (b) linear HDPE at 𝛾̇𝐴=248 s−1. The flow is from top to bottom.

6.2. Pressure Drop in Contraction Flows

In capillary flow of molten polymers generated by means of a capillary rheometer, there is a large pressure drop associated with the flow in the entrance (major) and exit (minor) regions. If the pressure in the reservoir is the quantity measured to determine the wall shear stress in the part of the capillary where flow is fully developed, these two excess pressure losses, collectively known as end pressure, should be taken into account (Dealy and Wissbrun [57]). The end pressure is also used to estimate the extensional viscosity of polymers, a method well practised in industry (Cogswell [96, 97]). Therefore, it is significant to understand the origin of this excess pressure and consequently be able to predict it.

Figure 10 plots the axial pressure variation in a capillary die including both its entrance and exit regions. It can be seen that the total pressure drop, Δ𝑃, consists of two components and may be written as Δ𝑃=Δ𝑃_Cap+Δ𝑃_End=Δ𝑃_Cap+Δ𝑃_End+Δ𝑃_Exit, (86) where Δ𝑃 is the total pressure drop from the reservoir to the capillary exit, Δ𝑃_Cap is the pressure drop over the length of the capillary where the flow is fully developed, and Δ𝑃_End=Δ𝑃_Ent+Δ𝑃_Exit is the excess pressure drop due to the entrance and exit flows. Bagley [98] has suggested an empirical method of calculating this combined excess pressure Δ𝑃_End. The results from this procedure are customarily presented in terms of the Bagley correction or end correction, 𝑛_𝐵, defined as 𝑛_𝐵=Δ𝑃_End2𝜎_𝑤, (87) where 𝜎_𝑤 is the wall shear stress in the die for fully developed flow. In a capillary die of radius 𝑅 and length 𝐿 (see Figure 10), the shear stress, 𝜎_𝑤, can be related to the Bagley correction, 𝑛_𝐵, the overall pressure drop, Δ𝑃, and the derivative of the pressure, 𝑑𝑃/𝑑𝑧, along the die axis according to Bagley [98]: 𝜎_𝑤=𝑅2𝑑𝑃=1𝑑𝑧2Δ𝑃𝑛_𝐵.+𝐿/𝑅 (88) A graph of Δ𝑃 versus 𝐿/𝑅 for a constant apparent shear rate 𝛾̇_𝑎 is the well-known Bagley plot, where 𝑛_𝐵 can be found by extrapolation of straight lines to Δ𝑃=0.

Figure 10 

Schematic diagram of the excess pressure losses in flow through a capillary rheometer.

This important phenomenon exhibited by polymer melts had not been successfully predicted by any viscoelastic constitutive model. On the contrary, the early efforts by simple models gave less pressure losses than even the Newtonian fluids in contrast with experimental evidence.

The situation was recently overturned by the application of the K-BKZ/PSM or K-BKZ/L-T model for the two polyethylene melts of Section 3. Figure 11 shows the good agreement obtained for the two melts, namely, the LDPE melt [45] and the HDPE melt [46]. In doing the simulations, all effects have been accounted for, such as temperature (nonisothermal simulations), pressure dependence of viscosity (important for LDPE), and slip at the wall (important for HDPE). More information can be found in the two publications. Therefore, use of the K-BKZ model in its variant types (PSM or L-T) was successful in correctly simulating the pressure drops in rheometry, an important phenomenon that eluded proper prediction so far.

(a)

(b)

Figure 11 

Bagley plot for two polyethylene melts. The effect of 𝐿/𝐷 on the pressure at various values of apparent shear rate (all effects accounted for). Symbols are experimental data while lines are simulations with the K-BKZ/PSM model: (a) LDPE melt at 160°C [45], (b) HDPE melt at 190°C [46].

7. Conclusions and Future Prospects

The present paper reviewed the 50 years since the appearance of the K-BKZ model in the rheological community around 1962. This review followed and borrowed from Tanner’s article [9], which had reviewed the relevant publications for the 25th anniversary of the model in 1988. The numerical methods section followed and borrowed from Keunings’ article [65] in 2003. The updating of the works showed the following.(1)The K-BKZ model has had a strong influence in the rheological community in the last 50 years, with an increasing number of citations (reaching 2000) and a reputation for good results.(2)The interest by other groups peaked around the 25th anniversary of the article (1986–1990) with different groups around the world working towards a numerical implementation of the model and an emphasis on its predictions for real polymer flows.(3)There were important improvements of the model following the idea of irreversibility by Wagner et al. [28], a more complete description by Luo and Tanner [31], and a new damping function for planar flows by Olley [33].(4)The numerical implementation saw the emergence of new schemes that were more general (valid for recirculation patterns) [70] and were extended to transient flows (LIM) [74–76] and 3D flows [77]. To this effect the continuing increasing power of computers helped significantly.(5)Many works focused on simulating real data and phenomena associated with polymer melt flows in rheology and polymer processing, in most cases with considerable success [91].(6)Important new efforts focused on molecularly-based integral equations based on the Doi-Edwards theories [34–42].

The challenges ahead are likely to continue along the molecularly-based models, although efforts will also be directed towards numerical solutions of the K-BKZ/PSM or L-T variants for three-dimensional, time-dependent flows, perhaps with the Newton-Raphson method for better precision. The method developed by Dimakopoulos et al. [99] is extremely useful towards that end. The main advantage of integral models remains the inclusion of many relaxation times (or an infinite spectrum) without having extra degrees of freedom associated with the stresses.

As Tanner [9] wrote presciently in 1988, “I do not think the theory will fade away rapidly since it is the optimal type of single-integral equation if a sufficiently limited range of flows is to be modeled. Therefore I look forward to 2012 and the 50th anniversary of the KBKZ theory”. In the same spirit, we can also say that we are looking forward to 2022 and the 60th (Diamond Jubilee) anniversary of the K-BKZ theory.

Acknowledgment

Financial assistance from the Program for Basic Research “THALES” of the Ministry of Education of Greece is gratefully acknowledged.