Extensional Instability in Complex Fluids: A Computational Study

(1)

Extensional Instability in Complex Fluids : A

Computational Study

Author: M. A. Abdulrazaq Supervisor: Luca Brandt Examiner: Luca Brandt

Master’s degree project

Stockholm, Sweden, October 2020 Department of Mechanical Engineering KTH Royal Institute of Technology

(2)

i

Abstract

In this study, instability and failure in complex fluids (Elastoviscoplastic fluids1) is explored using the classic Considère ( ˙F < 0) and stress curvature (¨σ < 0) criteria. Employing the Saramito model, numerical simulation of the extensional protocol on non-Newtonian fluids is carried out. Validation is firstly performed (with a purely viscoelastic model) and in gen-eral found to be in agreement with previous works.

Parameter variation of the Bingham number (Bi), capillary number (Ca) and extension rate ( ˙) is then undertaken. It is found out that for Oldroyd-B based fluids, the stress curva-ture condition almost always occurs from inception of the flow for all cases. Additionally, increasing surface tension has a stabilizing effect on the extending fluid when it is below a critical value, above which it aids breakup. Increasing the yield stress, though, delays the onset of instability, but reduces the final length of the extending filament. At mild to high extension rates, the Considère criterion and the extension at the maximum stress are suit-able indicators of the final extension at strain-to-break(ST B). Furthermore, the rate of the of

necking instability till final breakup varies with the ST B at moderate to high ˙.

Keywords: Complex Fluids, Computational Rheology, Elastoviscoplastic Fluids, Insta-bility.

(3)

ii

Abstrakt

I denna studie undersöks instabilitet och misslyckande i komplexa vätskor (Elastoviskoplas-tiska vätskor2) med den klassiska Considère ( ˙F < 0) och stresskurvatur (¨σ < 0) kriterier. Genom att använda Saramito-modellen utförs numerisk simulering av det utökade protokol-let på icke-newtonska vätskor. Valideringen utförs först (med en rent viskoelastisk modell) och i allmänhet visar sig överensstämma med tidigare verk.

Parametervariation av Bingham-numret (Bi), kapillärnummer (Ca) och förlängningshastighet ( ˙) genomförs sedan. Det har visat sig att för Oldroyd-B-baserade vätskor, uppträder stresskrökn-ingstillståndet nästan alltid från början av flödet i alla fall. Dessutom har ökande ytspänning en stabiliserande effekt på den utsträckande vätskan när den är under ett kritiskt värde, över vilket den underlättar uppbrytning. En ökning av sträckgränsen fördröjer dock instabiliteten men minskar den slutliga längden på det utsträckta filamentet. Vid milda till höga utvidgn-ingshastigheter är Considère-kriteriet och förlängningen vid maximal spänning lämpliga indikatorer för den slutliga förlängningen vid spänning till brott (ST B). Vidare varierar

frekvensen av instabilitet i halsen till slutlig upplösning med ST B vid måttlig till hög ˙.

Nyckelord: Komplexa vätskor, Beräkningsreologi, Elastoviskoplastiska vätskor, Insta-bilitet.

(4)

Introduction

1.1 Background

Fluid is all around us. The air we breathe, blood flowing through our veins etc. Hence, it is not surprising that scientists have devoted effort to unravel the physics, behaviour, prop-erties etc. of fluids. Puzzling however is that most of these have concentrated on "well-behaved"1 fluids. Non-Newtonian fluids had been sidelined until recently, justifiably, due to their complex nature. More often than not, equations and concepts in rheology are "im-ported" from ideas of solid mechanics and fluid dynamics.

Paradoxically, most fluids ranging from biological fluids2, industrial3, geological 4 are complex (Saramito , 2016). Rather, almost all fluids (except of course air and water) show non-Newtonian properties. Unfortunately, most of the research over the years have been concentrated on deriving the material properties of Newtonian fluids.

From a more practical perspective, the material response and dynamics of complex fluids are necessary in several commercial processes including (but not limited to) coatings, oil re-covery, drag reduction and lubrication, film blowing, extrusion and fiber spinning (Mckinley & Hassager, 1999), ink jet printing (Fielding , 2011), atomization of fertilizers and pesticides and electro-spinning processes (Mckinley , 2005).

Obtaining a suitable means to extract the mechanical material properties of complex flu-ids (most commonly the extensional viscosity µ+_E) has been a long challenging issue among rheologists with several disappointing results obtained (Walters , 1992). Various efforts rang-ing from the melt extensiometer (Meissner & Hostettler , 1994), filament stretchrang-ing rheome-ter(Matta & Tytus , 1990; Tirtaatmadja & Sridhar , 1993), etc. were devised previously with the sole aim of accurately measuring the extensional viscosity of non-Newtonian fluids.

Indeed, over the years, the filament stretching procedure(Nieuwkoop & Czernicki, 1996; Yao & McKingley, 1998; McKinley & Sridhar , 2002; Sridhar et al., 1991; Tirtaatmadja & Srid-har , 1993; Spiegelberg et al., 1991) has turned out to be quite effective in obtaining the mate-rial property which was previously challenging for experimentalists.

1_{This refers to Newtonian fluids where the relationship between stress and strain rate is linear i.e. shear stress}

∝ shear rate

2_{Blood and biological tissues}

3_{Liquid foams, polymer solutions, hot metal forming (Saramito , 2016), film blowing, fiber spinning (Mckinley}

& Hassager, 1999), etc.

4_{Meteoric and tectonic dynamics of the earth}

(7)

2 CHAPTER 1. INTRODUCTION

1.2 Literature Review

A thorough review of relevant research works dealing with the stability of complex fluids in extension is somewhat ambiguous. Rather, one does find research that contribute to various parts of the "puzzle" - for example - stability in Newtonian fluids (without focus on non-Newtonian fluids), computational study of extensional complex fluids (with no relevance to its stability), etc. All of these works are crucial, as will be seen later in this study. Instead, a cursory survey is presented.

In Newtonian fluids, the pioneering research of Plateau (Plateau , 1873) paved way for research in stability in fluids. Rayleigh later showed (by a linear stability) analysis that a liquid column loses stability when the wave length of disturbances exceeds the length of the circumference of the filament (Strutt & Rayleigh , 1878). The analytic works of Weber (Weber , 1931) and Tomotika (Tomotika , 1935); and experimental research of Taylor (Taylor , 1934) also formed pioneering works in this respect.

Stability and breakup in non-Newtonian fluids (notably viscoelastic jets) has also been investigated. Early references are (Bousfield et al., 1986; Goldin et al., 1969) who investigated the effect of surface tension on the breakup of viscoelastic fluids. In extensional or stretching filaments, however, early contributions were made by (Chang & Lodge , 1971) evincing that elongational flow of viscoelastic fluids are more stable than that of Newtonian fluids. (Ide & White , 1976) also studied the linear stability of uniaxial extension in Maxwell fluids. An important inference from their work is that Newtonian fluid and a Maxwell fluid at Deborah numbers De > 1/2 will show ductile failure at some finite Hencky strain. Stability of linear 2-D flows of Olyroyd-B based viscoelastics was also investigated by (Lagnado et al., 1969). The results from this work showed that uniaxial extension in an infinite domain (based on a 2-D setup) is unstable. (Pearson & Connely , 1982) also analyzed failure in viscoelastic fluids in terms of the Considère criterion.

In more recent times, considering the exponential rise in computational resources, re-searchers have gradually turned to numerical computation of extension in non-Newtonian fluids. Due to the transient evolving nature of the fluid interface, early research efforts (Yao & McKingley, 1998; Kolte et al., 1997; Gaudet & Mckinley, 1998; Rasmussen & Hassager, 1996) adopted a Lagrangian approach where the mesh evolves with the fluid interface. In addition, most of the works used the finite element approach. A major drawback of this ap-proach is the need to continuously re-mesh (Rasmussen & Hassager, 1996; Yao & McKingley, 1998) which increases the computational time. Another disadvantage is that finite element simulations can be problematic in predicting the failure of filaments at high Hencky strains (Hassager et al., 1998). Even more recent are the research papers from (Hoyle & Fielding , 2016a; Fielding , 2011) on stability of non-Newtonian viscoelastic fluids which is referred to heavily in this work.

1.3 Objective

The primary objective of this study is to understand how input parameters - namely the sur-face tension, yield stress and extension rate affect the development of instability in complex fluids (using the Considère and stress curvature criterion) and their final failure strain - strain to break. The focus here would be the Oldroyd-B based fluids (which includes the Saramito model).

(8)

CHAPTER 1. INTRODUCTION 3

These goals outlined are accomplished via direct simulation of an exponentially extend-ing filament incorporatextend-ing novel computational techniques - MTHINC VOF and a semi-implicit low Re solver.

(9)

Chapter 2

Theoretical Framework

2.1 Extensional Protocol in Rheology

Broadly speaking, there are two basic flows that are traditional to non-Newtonian fluid me-chanics - shear flows and extensional flows. In classical rheological texts, the discussion on shear flows is prevalent (Bird et al., 1987; Tanner , 2000). Generally, shear flows are "weak" due to the fact that they are dictated by a more linear input (i.e. the fluid particles are linearly separated as time evolves) (Mckinley & Sridhar, 2002). Extensional flows on the other hand are dictated by a non-linear exponential input. Thus, impose more stringent kinematic con-ditions on the flow (especially during constant exponential separation), hence, several non linear features only become apparent in extension1_{. This partly justifies the superiority of}

an extensional protocol (Hoyle & Fielding , 2016a). Another interesting observation that de-scribes the diffrence between these two protocols is the impracticability of having an infinite extensional process. The idea of an homogeneous, uniaxial, infinite extension is more of an "academic fairytale" than reality. Shear flows (in principle), however, can be undertaken for as long as experimentalist decides (Malkin & Petrie, 1997)

The concept behind extensional rheometry has been explored by many researchers and is the idea behind the so called "Filament stretching rheometers".

2.2 The Oldroyd-B Model

The Oldroyd-B model (also called the convected Jeffreys model) is perhaps the simplest model for describing viscoelasticity in fluids (Yao & McKingley, 1998; Bird et al., 1987; Mckin-ley & Sridhar, 2002). It belongs to the class of phenomenological models2 which assumes the the fluid continuum to be composed of massless infinitely extensible dumbbells which describe the dynamics and interaction between the non-Newtonian component (or polymer) and the base solvent.

Due to the assumption of infinite extensibility of the constituent dumbbells, when the De is high (specifically when De > 0.5), the transient viscosity grows without bound. This un-physical phenomenon has been termed "extensional catastrophe" (Hoyle & Fielding , 2016a). This "coil-stretch transition" first described by (deGennes , 1974) has also been confirmed by various experimental studies (Perkins et al., 1997; Smith & Chu , 1998). Other more realistic

1_{e.g. strain hardening that occurs in the middle column of the polymeric fluid} 2_{Models derived empirically and not from "scientific first principles"}

(10)

CHAPTER 2. THEORETICAL FRAMEWORK 5

models like the Giesekus model, finitely extensible non-linear elastic model (FENE-CR), etc. have been proposed. These alternatives are neither discussed nor considered in this study.

2.3 Elastoviscoplasticity and the Saramito Model

It is likely to assume that yield stress fluids lack extensibility, and merely possesses just vis-cous and plastic properties. Perhaps, this justifies why many of the simulations for yield stress fluids, for example - (Balmforth et al. , 2010), ignore its existence and consequent ef-fects. However, recent observations and studies have proved the existence of significant extensibility experienced by some viscoplastic materials. Phenomenon such as the fore-aft asymmetry and rapid velocity increase observed experimentally by (Holenberg et al., 2012) for the flow of yield stress fluids around circular cross-sections are only explainable when the extensional property of yield stress fluids are considered. In yet another study, (Balmforth et al. , 2014) ascribed the discrepancies in the comparison between experimental and com-putational results in Kaolin (a yield stress fluid) to the possibility of extensional properties in the fluid. Even more recently, (Nelson et al., 2018) documented high extensibility in yield stress fluids. These studies all point to the fact that even yield stress fluids have intrinsic elastic properties which should be considered.

With these realizations, several attempts have been put forward to mathematically de-scribe the behavior of EVP fluids (while accounting for its elastic property). (Benitoet al., 2008) presented a constitutive model for EVP fluids which predicts solid viscoelasticity (e.g. like the Kelvin-Voigt model) below a critical stress and fluid-like viscoelasticity (e.g. like the Maxwell model) above this critical stress. (Saramito , 2007) proposed a tensorial constitutive model wherein the material assumes a Kelvin-Voigt solid before a critical stress is reached and an Oldroyd-B viscoelastic fluid above this critical stress. An important distinguishing feature between the two previously listed models (after yielding) is that the former does not assume any background solvent viscosity µs, however, the latter assumes a polymeric

substance consisting of a background Newtonian fluid µsand viscoelastic part µp.

To account for the shear thinning 3 behaviour (or the less commonly observed shear thickening 4) observed in some materials after yield, (Saramito , 2009) proposed another model which predicts a Kelvin-Voigt solid before yielding and a non-linear, viscoelastic, Her-schel–Bulkley fluid after yielding. The yielding criteria in this model (and also in (Saramito , 2007)) is determined by the von Mises criterion. Unlike some models for EVP developed previously, the Saramito Oldroyd-B and Saramito Herschel-Bulkley fulfill the thermody-namic framework for generalized standard materials defined in (Halphen & Nguyen , 1975) and based on the second law of thermodynamics5_{. Additionally, these models predict a}

smoother and continuous behaviour as the material transitions through the yield/critical stress i.e. τapplied < τcritical → τapplied> τcritical.

In this study, the Saramito Oldroyd-B model (which will henceforth be referred to as the Saramito model) is employed. This model has been extensively validated, and recently; has been adopted widely in the simulation of EVP fluids. (Rosti et al., 2018) employed the Saramito model in comparing the effect of the Bingham number on pressure driven flow in transitioning to turbulence. (Rosti et al., 2019) also utilized the Saramito model in the

simula-3_{like Carbopol (Piau, 2007) and liquid foams (Katgert et al., 2008)} 4_{for example (Iordan et al., 2008)}

5_{This property is based on the convexity of the free energy function and dissipation potential which}

(11)

6 CHAPTER 2. THEORETICAL FRAMEWORK

tion of multiphase polymers in Couette flows to examine the effect of varying capillary num-ber on the effective viscosity of the solution. (Izbassarov et al., 2018) developed a numerical code for simulating the flow of deformable and rigid particles in viscoelastic and EVP fluids using the Saramito model. (Cheddadi et al., 2011) also simulated the flow of an EVP fluid around a circular object using the Saramito model. Their results correlated quite well with experimental findings reported in (Dollet & Graner , 2007). Interestingly, they were also able to capture the fore-aft asymmetry and rapid velocity increase after the circular barrier, which was not predicted by purely viscoplastic models (like Bingham, Herschel-Bulkley etc.). Fur-thermore, (Fraggedakis et al. , 2016) used the Saramito model to simulate the physics of a sedimenting particle in Carbopol. Their results agreed quite well with experimental find-ings reported in (Holenberg et al., 2012).

2.4 Tensile Stress Measures

In research works on extensional rheology, it is not uncommon to observe the use of terms like apparent Trouton ratio, engineering stress etc. Though, the difference may seem subtle, the actual measurement by these stresses (or forces) could differ substantially. These different stress measures (Hoyle & Fielding , 2016a; Fielding , 2011) are discussed below

-Engineering Stress - The ratio between the axial extending force Fz(t) and the initial

cross sectional area A0 of the extending material. It is also common to divide this by the

extension rate ˙ to obtain the engineering stress growth coefficient µ+ eng -τ_eng+ (t) = Fz(t) A0 µ+_eng(t) = Fz(t) A0˙ (2.1) This stress measure does not adequately depict the material response to the imposed constant extension rate since the cross sectional area throughout the experimental procedure is fixed to A0, the material’s initial area.

Apparent Stress- In this case, the elongating force Fz(t)is normalized by the a uniform,

exponentially reducing area A(t) = A0L0/L(t) = A0e−. The assumption is based on a

uni-form extending material that does not experience necking. The force and stress are expressed respectively as -τ_app+ (t) = Fz(t) A(t) µ + app(t) = Fz(t) A(t) ˙ (2.2)

This is generally a better quantifier of the stress than the previous measure, though, it does not take into account the area variation once necking or instability sets in.

True Stress- This is obtained by dividing the force Fz(t)by the actual area A(z, t) at that

filament cross section. Since it does not assume uniform extension, the area varies along the filament axis’ length z i.e. A(z, t). The force and stress are then expressed as

-τ_E+(z, t) = Fz(t) A(z, t) µ + E(z, t) = Fz(t) A(z, t) ˙ (2.3)

This true stress represents the most complete and accurate definition of stress, and there-fore will be used throughout this study. The axial position for the measure will be at the mid

(12)

point of the filament i.e.

τ_E+(zmid, t) = Fz(t) A(z, t) µ + E(zmid, t) = Fz(t) A(z, t) ˙ (2.4) Where zmid= L(t)/2

2.5 Elongation Strategies in Extensional Rheology

Generally, the elongation of filaments can be accomplished in several ways. Following the discourse in (Kolte et al., 1997), they are enumerated below:

Type 1 experiment- In this experimental procedure, the end plates are separated exponentially with a constant extension rate, ˙ according to the equation

-L(t) = L0e− ˙t (2.5)

The axial force needed to achieve this constant extension rate ( ˙) is measured with respect to time. The transient Trouton ratio is calculated thus:

T r ( ˙, ) = Fz A()µ0˙

(2.6) Where Fz is the axial extending force, A() = A0e− is taken to be the cross-sectional

area at the mid-filament section (assuming an homogeneous extension). The Hencky strain is = ˙t. This approach is usually adopted by most experimentalist due to its ease.

Type 2 experiment- This approach is similar to the experimental procedure highlighted in Type 1, however, the elongation rate at the mid section of the extending filament is also noted i.e. ˙s= − 2 Rs dRs dt (2.7)

In this case, the extension rate ˙sis a variable and changes to reflect the non-homogeneity in

the extending material. In this case, the Trouton ratio is then -T r ( ˙s, s) =

Fz

Asµ0˙s

(2.8) Fzis the axial force driving the extension process, A(s)is taken to be the cross-sectional

area at the mid-filamnet section (without the assumption of an homogeneous extension). The Hencky strain in this case is (which is not equivalent to the Type 1) is given as

-s= Z t t0 ˙sdt0 = −2 ln D(t) D0 6= ˙st (2.9)

Type 3 experiment- In this technique, the focus is on ensuring the mid-filament reduces exponentially according to the homogeneous relationship

-Rs(t) = R0e−

1

(13)

The extension rate ˙shere is kept constant while the end plates are varied accordingly. The

Trouton ratio in this case based on the ˙sat the midfilament is

-T r ( ˙s, s) =

Fz

Asµ0˙s

(2.11) The Hencky strain is then

-s= ˙st (2.12)

The approach is advantageous as the Trouton ratio is will be quite close to the "True" Trouton ratio, however, it is arduous to achieve experimentally.

In line with experimental procedure, Type 1 is adopted in this study.

2.6 Instability Signatures in Non-Newtonian Fluids

The concept of instability is a precursor to subsequent break-up and imminent failure in the fluid. In the literature, there are generally three prototypical "extension-type" flows that are used to study the breakup in fluids (Mckinley , 2005; Basaran , 2002). Reference is made to figure 2.1.

Figure 2.1: Prototypical flows for studying breakup fluids. From left to right - (a) continuous jetting instability (b) drops from a nozzle (c) necking and breakup of a liquid bridge. Figures reproduced from (Mckinley , 2005)

For Newtonian fluids, these flows have been reviewed and studied in details. Relevant reviews are presented in (Basaran , 2002; Eggers & Villermaux , 2008). However, for complex fluids, there has been relatively (and understandably) less research devoted. Worse still is even the study of instability in EVP’s. In reality, similarities occur in EVP’s and purely vis-coelastic fluids. This will motivate a strong assumption which will form the foundation of this study. The prospect of extending the present theories on instability (or necking, thinning, break-up etc.) in viscoelastic models to EVP’s. In the literature, there are three commonly discussed modes of failure or instability in viscoelastic materials. As viscoelastic materials tend to show properties which "fluctuate" between viscous, liquid-like materials and elastic, solid-like materials, it will not be unusual to notice an extension of ideas on instability in

(14)

solid or liquids to viscolastic materials6_{. We now briefly summarize three basic benchmark}

ways of determining failure in non-Newtonian fluids.

2.6.1 Considère Criterion

Historically, the Considère criterion (Considere , 1885) was originally developed as a con-cept in solid mechanics to predict the failure of materials. The pioneering work of (Vin-cent , 1960) on polymeric tensile test paved the way for its application to non-Newtonian fluid mechanics. Since then, several other works have applied this failure criterion to non-Newtonian fluids, for example - (Pearson & Connely , 1982) in estimating the extensional viscosity, (Cogswell & Moore , 1974) in understanding the behaviour of polymer melts in the "solid-like" rubbery state etc.

Formally, the Considère failure criterion states that a material undergoing uniaxial ho-mogeneous elongation is stable provided the present strain is lower than the strain when the maximum axial force (∂F/∂ = 0) occurs. In other words, stability is guaranteed when < F, where F is the extension when ∂F/∂ = 0. To further formalize this mathematical

expression,

dF d =

d

d(σzz− σrr· A()) (2.13)

For stability, dF/d > 0. Besides, an homogeneous extension presumes that A() = A0e−. Consequently, dF d = d d (σzz− σrr) · A0e −_{≥ 0} (2.14) Non dimensionalizing the expression above by µ0˙0and further simplifying,

dT r+ d − T r

+ _{≥ 0 ⇒} d ln T r+

d ≥ 1 (2.15)

Where the transient Trouton ratio is T r+_{= (σ}

zz− σrr) / (µ0˙0)

The implication of the above formulation is that for stability to be maintained, the tran-sient Trouton ratio must increase at least exponentially with the true strain (to balance or surpass the effect of exponential decrease in area). Otherwise, the material would undergo dynamic failure (Mckinley & Hassager, 1999). One issue with the expression above is the as-sumption of homogeneous extension which is strictly not valid for a realistic filament stretch-ing experiment due to the imposition of the no slip condition at the end plates.

While the stress curvature mode (discussed in the next section 2.6.2), is associated with a slower, "fluid-like" transition to instability; the Considère criterion is defined by a rapid occurrence of instability to final failure. Due to its rapid development to failure, it is rational to associate this failure criterion with higher elongation rates ˙ (or larger relaxation time λ or larger W e). At higher ˙, the flow dynamics evolves rapidly. Here, the work done to extend the filament is stored majorly elastically and the time scales (quantified by the inverse of extension rates 1/ ˙) are so small such that molecular relaxation does not occur. A visual in-terpretation of this evolving fast dynamics is the relative discontinuity and sudden decline in the stress plot (as would be seen in chapter 5 on results) when this failure mode is dominant.

(15)

From the major criticisms of this criterion for viscoelastic materials is its independence from the strain rate ˙ since the criterion expresses a relationship between force, F and exten-sion, which is an "elastic" relationship (Hutchinson & Obrecht , 1978). Moreso, the rate of necking - whether rapid or gradual - cannot be deduced from this criterion (Fielding , 2011).

2.6.2 Stress Curvature Criterion

This mode is elucidated in (Hoyle & Fielding , 2016a; Fielding , 2011). This mode is defined by the downward curvature of the constitutive stress plot7_{. It is associated with a slower,}

ductile failure of the material where the transition of the EVP to final failure is gradual rather than sudden and rapid (as in the Considère criterion).

This criterion can be defined in terms of a constant extension rate protocol or a constant stress protocol. In this study, we focus on a constant extension rate protocol. Consequently, further discussions will be tailored towards that. This criterion proposes that (initial) insta-bility will occur when the second differential of the stress function of the respective constitu-tive stress relationship as a function of time (or strain)8is negative i.e (¨σ ≤ 0). The derivation of this criterion is based on linear stability analysis and it is also independent of the con-stitutive model (e.g Oldroyd-B, FENE-CR etc.) adopted. Conversely, the EVP is stable if ¨

σ > 0.

To further explain, if stress σ is expressed in terms of time t, then, for stability of the non-Newtonian fluid in extension, the convex relationship equation (2.16) should hold for the constitutive stress equation

-σ

(αtn+ (1 − α)tn+1) <

σ

(tn) + (1 − α)

σ

(tn+1) (2.16)

where α ∈ [0, 1] and tn≤ tn+1

For non-Newtonian materials exposed to small or medium values of constant strain rate ˙(or small relaxation time λ or small W e = ˙λ), this form of instability is expected to dictate and dominate the dynamics of breakup of the non-Newtonian fluid. More so, at this condi-tion (of small or medium strain rate), it is expected that the break-up mechanism is slow and steady as opposed to being rapid.

One important qualitative point to derive from the stress condition discussed is the speed of necking with regards to the slope of the constitutive stress curve. From (Hoyle & Fielding , 2016a), the flatter the constitutive curve the faster the rate of necking at that section and conversely a higher stress slope curve implies slower rate of necking. This concept will be used in the discussion (in section 5) when comparing the rates at which necking occurs with regards to the parameter variation. It is also important to mention that in all cases, even when the force criteria dictates the actual failure, the stress curvature will have occurred earlier. In this case, the force criterion can then be perceived as an amplifier to an already existing instability process (Hoyle & Fielding , 2016a).

2.6.3 Elastic Considère Criteria

The elastic Considère mode is developed and modified from the conventional Considère mode to apply to viscoelastic fluids. Therefore, the derivative of the force is applied to the

7_{This definition is valid for a constant extension rate protocol which is the focal point of this study}

8_{Converting between strain and time derivative (at constant strain rate) can be achieved thus: = ˙t ⇒ d =}

(16)

elastic strain rather than the general strain. Mathematically, stability is ensured when -∂F ∂ elastic > 0 (2.17)

Stability analysis from (Hoyle & Fielding , 2016a) indicates the eigenvalue of this failure criterion has an amplitute of order O(G/µ) where G is the viscoelastic modulus. Expectedly, this mode is quite similar to the classic Considère criterion. In short - ∂F_∂

elastic → ∂F

∂ as

λ ˙ → ∞.

The authors of the original research on the elastic Considère criteria (Hoyle & Fielding , 2016a) do attest some ambiguity in determining the elastic Considère criteria (especially from an experimental perspective). Moreover, the elastic Considère criteria is similar to the classic Considère criteria. Due to these reasons, only the latter is opted for in this study.

It is emphasized that there is an important difference between the strain-to-break exten-sion (which is designated as ST B) and the strain deduced from the stress curvature failure

mode (which is assigned as σ). The σdoes characterizes the ST B(Hoyle & Fielding , 2016a)

and is thus a pointer to it, but it is definitely not equivalent (in quantitative terms) to the ST B.

To drive home the point, the σ according to the Oldroyd-B based models occurs from the

inception of flow at very small values of ˙, however, it is inconceivable to think that the ST B

also occurs at the inception of flow. Rather, the ST B is delayed at low strain rates, though,

instability would have occurred far earlier.

2.7 Modes of Occurrence of Instability and Failure

The signatures of instability in filament stretching indicate that failure is imminent. How-ever, this does not address the question of where it occurs. That is the direction of discussion in this section. Some of what is mentioned here is derived from (Spiegelberg et al., 1991; Yao & McKingley, 1998). It may be commonly assumed that failure occurs at the mid-plane of the elongating sample, but this is not always true. Failure/necking/instabilities may be - (a) Necking or capillary instabilities or (b) Endplate instability arising from the deformation of the fluid at a solid boundary interface. These common modes are explained below:

Capillary failure of a viscous filament

Often, this is the mode of failure for viscous, Newtonian, capillary driven flows. Due to the increasingly large capillary pressure at the mid plane, the filament section "necks" and separates into two distinct parts. In the absence of capillarity in Newtonian fluids, however, it has been shown (Renardy , 2004; Hassager et al., 1998) that theoretically the sample will undergo infinite extension.

Elastic breakup and rupture

This mode of failure is experienced by non-Newtonian fluids which are strongly "strain-hardening". At long times from the extension process, strongly elastic fluids experience "strain-hardening" where the central part of the fluid undergoing extension becomes increas-ingly resistant to extension. Various viscoelastic models (e.g. Oldroyd-B, Geisekus etc.) are capable of capturing this "strain hardening" phenomenon.

For the Oldroyd-B model, the strain hardened section is more axially uniform than its Newtonian counterpart. This axial uniformity and reduced curvature tends to obliterate the

(17)

effect of capillary at the midsection. It is hence not surprising that via linear stability analysis, the filament is stable to failure or perturbations when the elastic extension (De) is ≥ 0.5 (Olagunju , 1997; Hoyle & Fielding , 2016a). This result is independent of the surface tension at the interface. At this high De, the mode of purely elastic failure will be the dominant mode. In truth, at high De, the instability would likely be dictated by the force criterion. Some authors (Mckinley & Sridhar, 2002; Yao et al., 1998) have opined that a higher De actually speeds up the process of instability and failure through this mode due to elastic recoil near the end regions.

Endplate instability and decohesion

In the previous section, it was explained that the Oldroyd-B model at high elastic strain is theoretically stable to perturbation. Another mode of failure which is also preponderant in elastic fluids is via decohesion or adhesion instability at the end plates. Due to rising pressure gradients and increased elastic tensile stresses at the end plates, a progressive inward flow directed at the extending filament is induced. This causes (and expedites) local buckling at that section.

Further more, it is important to mention that a corollary effect to this mode of instability is the observation of "fluid reservoirs" at the endplates, where the effect of exponential increase in length and progressive strain hardening of the mid filament section leads to drainage of fluid at the end plates birthing the formation of quasi-static fluid reservoirs at the end plates (reference to Figure 2.2) which finally ends in the formation of a "dimple" at the free surface of the end plates. This effects are well supported by experimental observations (Spiegelberg & Mckinley, 1996) and numerical computations (Yao & McKingley, 1998; Kolte et al., 1997). It has also been hypothesized by (Gupta et al., 2000) that this failure mode can be controlled by changing the shape and geometry (to reduce the curvature) of the end plates.

2.8 Problem Description and Idealisation

In this study, the extensional procedure adopted is described in (Yao & McKingley, 1998; Mckinley & Sridhar, 2002). The schematic of this type of extension is presented in Figure 2.3. An initial fluid sample of length L0present between two parallel and coaxial plates are

sepa-rated instantaneously by imposing a unidirectional axial velocity V0. The extension rate ˙0in

the axial direction is constant. This condition translates to an exponentially increasing axial velocity. Mathematically, this imposed velocity field and tensorial velocity-gradient can be expressed respectively as:

u = V0   0 0 1   ∇.u = 1 2˙0   −1 0 0 0 −1 0 0 0 2   (2.18)

Due to the imposed kinematics described above, the length of the sample increases ex-ponentially with time and the cross sectional width (or diameter) decreases similarly in an exponential manner with time i.e L(t) = L0 exp( ˙0t) and D(t) = D0 exp(− ˙0t). This is a

consequence of the conservation of mass for the extending filament. The actual extensional deformation of the filament at time t is quantified by the Hencky strain (or the logarith-mic strain). The generally used engineering strain would be unsuitable because of the large deformation which is typical in extensional rheometry. The Hencky strain takes into cogni-sance the previous deformation of the extended fluid. Thus, it is the summation of all the

(18)

Figure 2.2: Mode of failure due to endplate instability. Pictorial illustrations from the simu-lations of (Kolte et al., 1997), (Yao & McKingley, 1998) and this present study respectively.

(19)

Figure 2.3: Procedure for extensional chosen in this study. Reproduced from (Yao & McK-ingley, 1998)

deformation experienced by the fluid. Mathematically, the Hencky strain is expressed as - = Z t t0 ˙ 0dt0= ˙0t = ln L(t) L0 = −2 ln D(t) D0 (2.19) As has been mentioned earlier, one of the primary motives behind extensional rheometry is to obtain the material function of the non-Newtonian fluid. Chiefly, the transient exten-sional viscosity is desired. In the literature, there occurs a myriad of technical terms to qualify this viscosity - transient extensional viscosity (Mckinley & Sridhar, 2002; Kolte et al., 1997), tensile-stress growth coefficient (Mckinley & Sridhar, 2002), the elongational stress growth viscosity (Kolte et al., 1997). Generally, they imply the same thing and are quantitatively equivalent. Often, the transient extensional viscosity is further divided by µ0, and called the

transient Trouton ratio9where µ0is the zero-shear viscosity (Kolte et al., 1997). Formally, the

transient extensional viscosity10µ+_Eand the transient Trouton ratio, T r+_are

-µ+_E( ˙0, t) = τ11− τ22 ˙ 0 T r+( ˙0, t) = µ+_E µ0 (2.20) At extended times, this coefficient may reach a steady-state. In this case:

µE( ˙0) = limt→∞µ+_E( ˙0, t)

(2.21) It is necessary to mention that the fluid considered in this study are the so called Boger fluids where the zero-shear viscosity µ0is assumed constant.

9_{This ratio is equal to three in Newtonian fluids and in the limit of very small extension rates} 10_{The "+" subscript is put to characterize the time varying nature of this viscosity measure}

(20)

Chapter 3

Governing Equations and Models

In modelling the extensional dynamics of the EVP fluid, we utilize three mathematical for-mulations, which are enumerated below. This is while keeping in mind the recommendation in (Yao & McKingley, 1998) that an accurate simulation of extensional free surface flows in non-Newtonian fluids should include:

• A suitable model to capture the free surface of the fluid

• A flexible algorithm to capture the transient extensional effects, and • An appropriate constitutive model

Thus, our mathematical formulation is based on these models which we explain after-wards

• Mass conservation equation • Cauchy momentum equation • non-Newtonian constitutive model

• Multidimensional tangent of hyperbola for interface capturing (MTHINC VOF method)

3.1 Mass Conservation Equation

We consider the mass conservation equation as an incompressibility constraint on the momentum equation and write it in a non conservative form as

-∇.u = 0 (3.1)

3.2 Cauchy Momentum Equation

The force balance in the non-Newtonian fluid is achieved by resorting to the Cauchy mo-mentum equation. Here, the tensorial stress field, velocity field and the scalar pressure are coupled by balancing the contributions from the transient advective/inertia term, the pres-sure gradient, the diffusion term, the non-Newtonian stress (sometimes referred to as the extra stress), surface tension force and an external force field (most commonly gravity is con-sidered).

(21)

16 CHAPTER 3. GOVERNING EQUATIONS AND MODELS

Its form here is in a non-conservative form similar to the format written in (Dodd & Ferrante, 2014) and (Li et al., 2012). This non-conservative format is also recommended by (Tryggvason, Scardovelli & Zaleski , 2011) in the treatment of varied density multiphase flows. ∂u ∂t + ∇. (u u) = 1 ρi −∇p + ∇. µi ∇u + ∇uT + ∇.τ + fσ+ g (3.2) In the equations above, u represents the velocity field, p is the isotropic scalar pressure, ρi and µi are the density and dynamic viscosity of the fluid respectively. The EVP medium

is taken to be i = 1, while the external fluid medium (air) is i = 2. In this work, a density ratio between the EVP and air is set to 500. Furthermore, τ is the extra stress tensor obtained considering only the non-Newtonian part of the constitutive model, fσdenotes the force per

unit mass due to the surface tension and g is an imposed forcing term.

3.3 Non-Newtonian Constitutive Model

To describe the nonlinearity between the deformation rate and the applied stress, the Saramito model (Saramito , 2007) capable of describing the extensible property of yield stress fluids is used and presented below

-λ ˙τ + max 0,|τd| − τ0 |τd| τ − 2µmD(u) = 0 (3.3)

In the equations above, an Eulerian framework is employed since it is expected that dis-placements will be significant. This necessitates an invariant derivative of the stress tensor i.e. the conventional Lagragian derivative is not appropriate. The Gordon Schowalter material derivative (Gordon & Schowalter , 1972) is used, and is defined as

-˙ τ = ∂τ

∂t + (u.∇τ ) + τ W (u) − W (u)τ − a (τ D(u) + D(u)τ )

where ˙ = D(u) = ∇u + ∇uT_{/2 is the rate of deformation tensor and W (u) =}

∇u − ∇uT_{/2 is the vorticity tensor. For computational studies, the upper convected}

deriva-tive where a = 1 is often adopted (Saramito , 2016).

Returning to equation (3.3), λ is the relaxation time, τddenotes the yielding criteria which

we (may) obtain from the von Mises criterion (Shaukat et al., 2012), τ0 is the yield stress1of

the viscoplastic fluid, µmrepresents the dynamic viscosity of the polymeric part of the fluid

and G = µm

λ is the relaxation modulus of the polymeric fluid.

3.4 MTHINC VOF Scheme

To capture the dynamics of the interface of the EVP and the external fluid, we use the recently developed multidimensional tangent of hyperbola for interface capturing method, which we refer to as MTHINC VOF scheme. This was originally proposed by (Xiao et al., 2005) and extended to multidimensions (or 3D in this case) by (Li et al., 2012).

We define a hyberbolic tangent function H (x, t) (where x = [X, Y, Z] ) which approx-imates the indicator function to describe the discontinuous jump in properties across the

1_{We define the yield stress in fluids even while acknowledging the debate about the existence of yield stress}

(22)

CHAPTER 3. GOVERNING EQUATIONS AND MODELS 17

interface of the varying fluids. We further define the primary fluid domain of interest as Ωf l.

Therefore, x ∈ Ωf lwithin the fluid domain and x ∈ Ωexat the external phase. Thus;

H (x, t) =

1 x ∈ Ωf l

0 x ∈ Ωex

The advection of H (x, t) is governed by the equation for the Lagrangian derivative

dH(x,t)

dt in an Eulerian framework. This is written as:

dH (x, t)

dt =

∂H (x, t)

∂t + u.∇H (x, t) = 0 (3.4)

As before, we define u as the velocity vector. φ (x, t) will denote the averaged value of the indicator function within a defined control volume δΩ. This is easily obtained by applying an integration on the indicator function as shown below

-φ (x, t) = 1 |δΩ (x, t) |

Z

δΩ(x,t)

H x0, t dV0 (3.5)

Combining (3.4) and (3.5), we could rewrite the equation (in divergence form) as: ∂φ

∂t + ∇. (uH) = φ∇.u (3.6)

We note that for the grids present at the boundary interface ∂Ωf l, the values of φ (x, t)

range between 0 ≤ φ (x, t) ≤ 1. Otherwise, typically2, φ should be 0 or 1.

The actual indicator function H◦(x, t)defined over the whole domain (including Ωf land

Ωex) is approximated by a hyperbolic tangent function:

H (x, t) ≈ H◦(x, t) = 1

2[1 + tanh (β (g (x) + d))] (3.7) g (x)is a function defining the surface in multidimensional space. This defining function could be linear or quadratic.

A rigorous analysis of the details of the linear and quadratic functions can be found in (Xiao et al., 2005; Li et al., 2012). β is a parameter to determine the sharpness of the transition between Ωf l and Ωexdomains. Ideally, we expect a discontinuous change in the fluid

prop-erties at ∂Ωf l. This would require an infinitely small mesh size which is not computationally

feasible. The recommendation of (Xia et al., 2011) where β = 2.5 is acknowledged. However, for practical purposes, the value of β = 1.0 is applied in this study. This variation is not ex-pected to change the results significantly. More so, some previous simulations (for example (Rosti et al., 2019)) have also adopted the value of β to be unity.

2_{In reality, φ would be 0 + ε or 1 − ε where varies (but nevertheless should ideally be a small number). This}

(23)

Chapter 4

Computational Methodology

We now shift our discussion to the computational schemes used in discretizing the equations we presented earlier - equations 3.1 to 3.7.

Combining all the equations, there are 10 unknown quantities which are to be obtained - the tensorial non-Newtonian stress field, which is a symmetric second order tensor, thus only 6 values are unknown - τxx, τxy, τxz, τxy, τyy, τzz. Additionally, there are 3 components

for the velocity vector u − [u, v, w], the scalar pressure p and the volume of fluid parameter φ.

A staggered uniform grid system (Ferziger & Peric, 2002; Hirsch , 2007) is applied where the flow velocity - u, v, w are represented at their respective cell faces and the other flow parameters, namely - φ, τxx, τxy, petc. are represented at the cell centre.

The numerical solution for the governing equations is firstly based on the formulations in (Chorin, 1968; Temam, 1969) for single phase flows, and (Harlow & Welch , 1965) for interfacial free surface flows, sometimes all jointly referred to as the projection method for the Navier Stokes equation. All these combinations are then applied to the specific case studied here.

The time advancement of the equations is based on a hybrid-explicit scheme via a frac-tional step of the Adams Bashforth scheme (Kim & Moin, 1985; Dodd & Ferrante, 2014; Izbas-sarov et al., 2018). D(un) = 1 2 h ∇un+ (∇un)Ti (4.1) dun dt = −∇. (u n_un₎ | {z } Fn adv −∇p n ρn + 1 ρn∇. (2µ n_D(un₎₎ | {z } F_{dif f}n + σ ρnκ n_∇Cn | {z } F_capin + 1 ρn∇.τ n | {z } Fstrn (4.2)

The equation above is in semi discrete form and time advancement is from n → n + 1. The equation handles below - Fn

adv, Fdif fn , Fcapin and Fstrn represent the force per unit mass

due to the inertial-advective term, diffusve term, surface tension and non Newtonian stress (obtained from the constitutive equations in (3.3)) terms respectively.

The surface tension force is computed as fσ = σκ∇φ based on the continuum surface

approach (Brackbill et al., 1992). σ is designated as the surface tension constant between the two fluids, κ is the interface curvature and φ is the value for the volume fraction of fluid.

In advancing from time n → n + 1, firstly compute an intermediate velocity at n +1₂ with the pressure term excluded.

(24)

CHAPTER 4. COMPUTATIONAL METHODOLOGY 19 un+12 − un ∆t = 3 2 h

F_advn + F_{dif f}n + F_capin + F_strn+1 i

−1 2

h

F_advn−1+ F_{dif f}n−1+ F_capin−1+ F_strn i

(4.3) Thereafter, correct the velocity field via the imposition of the solenoidal velocity field invoking the pressure correction procedure in (Dodd & Ferrante, 2014) which is a two step procedure

Firstly, solve a Poisson equation for the pressure -∇2pn+1= ∇. 1 − ρ0 ρn+1 ∇p◦ + ρ0 ∆t∇.u n+1 2 (4.4)

Then correct the velocity with the newly obtained pressure -un+1 = un+12 + ∆t 1 ρ0 ∇pn+1+ 1 ρn+1 − 1 ρ0 ∇p◦ (4.5) Note that p◦ is a linear extrapolation of the pressure at time step n + 1 based on the previous two time steps n and n − 1, therefore p◦= 2pn− pn−1_{. Also, ρ}

0 = min(ρ1, ρ2)

For the stress equations, use an Adams Bashforth scheme i.e. the time advancement of the stress equation is similar to the momentum equation. It is also necessary to mention that the stress τ used in the semi discrete equation (4.2) is obtained via the procedure outlined below. The equations (3.3) is rewritten in a more convenient and workable form as:

Dτ Dt = 1 λ max 0,|τd| − τ0 |τ_d| τ − 2µmD(u) | {z } τn sar (4.6) Where D

Dtis the upper convected derivative operator obtained from the Gordon Schowal-ter derivative when a = 1 as mentioned in Section 3.3. It is also necessary to mention that the above equations should be taken as discrete equations applying to the current time step n. The superscript n was removed from the terms of the equation to ensure clarity. The upper convected derivative for the stress term is now explicitly reproduced below as

-Dτ Dt = ∂τ ∂t + [u.∇τ ] − [τ .∇u] − h (∇u)T .τi (4.7)

With the more convenient form (4.6), in addition to the expression above, the semi discrete equation for stress in a form more suited for an explicit time integration is

-dτn dt = − [u n_.∇τn_{] + [τ}n_.∇un_{] +}h_(∇un₎T_.τni | {z } τn adv +τ_sarn (4.8)

Recall that the last term in equation (4.8) - τn

sar is the additional stress term defined in

(4.6) which is specific to the constitutive model used. Finally, apply the Adams Bashforth scheme as below to integrate equation (4.8)

-τn+1− τn ∆t = 3 2 τ n adv+ τ n sar − 1 2 τ_advn−1+ τ_sarn−1 (4.9) All spatial derivatives are discretized with second order central difference, except the ad-vection term - u.∇τ in the Gordon Schowalter derivative of the non-Newtonian constitutive

(25)

20 CHAPTER 4. COMPUTATIONAL METHODOLOGY

equation 3.3. The absence of a diffusion term in this equation poses a problem of uncon-ditional instability. Different solutions have been proffered. One common approach is to use an upwinding scheme for the velocity (Dupret & Marchal , 1986). This however intro-duces artificial dissipation to the system (Dupret & Marchal , 1986). A third order upwind scheme turns out to be a very viable option (Min et al., 2001). Perhaps, a more preferred op-tion is to use an explicit fifth order weighted essential non oscillatory (WENO) scheme (Shu , 1999) which is adopted in this work (The WENO scheme is explained in further details in a subsequent section). It is relatively less computationally expensive and it also has similar accuracy and performance. It has also been applied successfully by (Izbassarov et al., 2018; Rosti & Brandt , 2017)

The volume-of-fluid method which is adopted is based on the work of (Xiao et al., 2005), with improvements by (Xia et al., 2011) and (Li et al., 2012). The advection equation (3.6) is integrated for the MTHINC VOF via the procedure outlined by (Li et al., 2012), also used successfully by (Rosti et al., 2019). The directional splitting approach, used here, involves a sequential discretization of equation (3.6) applied in three stages based on the 3D Cartesian coordinates -                 φ?_ijk= φn_ijk− _∆x1 f_i+n 1 2,j,k − fn i−1₂,j,k +_∆x∆tφ?_ijk u_i+1 2,j,k − u_i−1 2,j,k

φ??_ijk= φ?_ijk− _∆y1 g_i,j+? 1 2,k − g? i,j−1₂,k +_∆y∆tφ??_ijk v_i,j+1 2,k − v_i,j−1 2,k

φ???_ijk = φ??_ijk− _∆z1 h??_i,j,k+1 2 − h?? i,j,k−1₂ + ∆t_∆zφ???_ijk w_i,j,k+1 2 − wi,j,k− 1 2 φn+1_ijk = φ???_ijk − ∆t φ?_ijk _u i+ 1₂,j,k−ui− 1₂,j,k ∆x + φ??_ijk _v i,j+ 1₂,k−vi,j− 1₂,k ∆y + φ???_ijk _w i,j,k+ 1₂−wi,j,k− 1₂ ∆z (4.10) The φ field should be corrected to ensure it is solenoidal. This is achieved by applying the last equation in (4.10)

The fluxes fn, g?and h??are calculated through an integration (with respect to space and time) of the of the product of the indicator function and the appropriate velocity

-               f_i±n 1 2,j,k = 1 ∆y∆z Z δtn Z ∆y Z ∆z (uH◦)_i±1 2,j,kdy dz dt g_i,j±? 1 2,k = 1 ∆x∆z Z δt? Z ∆x Z ∆z (uH◦)_i,j±1 2,kdx dz dt h??_i,j,k±1 2 = 1 ∆x∆y Z δt?? Z ∆x Z ∆y (uH◦)_i,j,k±1 2 dx dy dt (4.11)

The integration of the fluxes (with respect to time) can be converted to one solely based on spatial variables. More so, this integration should be along the upwind propagation path of the velocity field. For instance, the pathway along the x direction will be defined as:

   ∆x+= h x_i+1 2 − ∆tu_i+1 2,j,k, xi+ 1 2 i u_i+1 2,j,k ≥ 0 ∆x−= h x_i+1 2, xi+ 1 2 − ∆tui+ 1 2,j,k i u_i+1 2,j,k < 0 (4.12) The numerical fluxes can finally be computed by making reference to the series of formulae defined in equation (4.13) below

(26)

-CHAPTER 4. COMPUTATIONAL METHODOLOGY 21                                              f_i+n 1 2,j,k = 1 ∆y∆z Z ∆x+ Z ∆y Z ∆z (H◦)n_i,j,kdx dy dz u_i+1 2,j,k ≥ 0 f_i+n 1 2,j,k = − 1 ∆y∆z Z ∆x− Z ∆y Z ∆z (H◦)n_i+1,j,kdx dy dz u_i+1 2,j,k < 0 g?_i,j+1 2,k = 1 ∆x∆z Z ∆x Z ∆y+ Z ∆z (H◦)?_i,j,kdx dy dz v_i,j+1 2,k ≥ 0 g?_i,j+1 2,k = − 1 ∆x∆z Z ∆x Z ∆y− Z ∆z (H◦)?_i,j+1,kdx dy dz v_i,j+1 2,k < 0 h??_i,j,k+1 2 = 1 ∆x∆y Z ∆x Z ∆y Z ∆z+ (H◦)??_i,j,kdx dy dz w_i,j,k+1 2 ≥ 0 h??_i,j,k+1 2 = − 1 ∆x∆y Z ∆x Z ∆y Z ∆z− (H◦)??_i,j,k+1dx dy dz w_i,j,k+1 2 < 0 (4.13)

It is mentioned that in equation (4.10) to (4.13), the symbols ? and ?? represent the inter-mediate time steps between n and n + 1, such that |tn+1− t??_{| = |t}??_{− t}?_{| = |t}?_{− t}n_|.

The relevant fluid properties (density and viscosity scalar fields) should then be updated after obtaining the volume-of-fluid parameter φ at time n + 1 i.e. φn+1_{. This is actualized by}

the formulae:    µn+1_i,j,k = 1 − φn+1_i,j,k µ2+ φn+1_i,j,k µ1

ρn+1_i,j,k =1 − φn+1_i,j,kρ2+ φn+1_i,j,k ρ1

(4.14)

4.1 Implicit Formulation for Low Reynolds Number Flow

The extensional protocol adopted in this study assumes a low Reynolds number flow to describe the stretching process (Yao & McKingley, 1998; McKinley & Sridhar , 2002; Sridhar et al., 1991; Tirtaatmadja & Sridhar , 1993; Spiegelberg et al., 1991). It is not uncommon to even assume a Stokes flow (Hoyle & Fielding , 2016a; Fielding , 2011; Hoyle & Fielding , 2016b)1_{. This imposes some constraints on the computation as the differential equation}

becomes stiff due to the viscous term F_{dif f}n in (4.2). For a stable solution for an explicit time stepping scheme (Tryggvason, Scardovelli & Zaleski , 2011), the time step is governed by -∆t = Re∆x₆ 2. This usually leads to prohibitively small time steps to achieve a stable solution. One proposed solution (Dodd & Ferrante, 2014) to this (which is adopted in this study) is to resort to a Crank-Nicolson second order semi-implicit method described here.

Firstly, the viscous term Fn

dif fin (4.2) is split into a variable part and a constant part

-F_{dif f}n = 1 ρn+1∇. 2µ n+1_D(u?_{) −→ ν} 0∇2u?+ 1 ρn+1∇. 2µ n+1_D(un_{) − ν} 0∇2un (4.15)

All the terms retain their previous meanings, however, ν0 = 1₂

µ1 ρ1 + µ2 ρ2 . It is also men-tioned that u?and un+12 are equivalent, and both symbols are used interchangeably.

1_{Though this assumption is questionable either at high extension rates or large strains when the convective}

inertia force u.∇u become significant (Yao & McKingley, 1998; Mckinley & Sridhar, 2002) and thus, a correction needs to be applied in these cases (Szabo & Mckinley, 2003; Yao & McKingley, 1998)

(27)

Thereafter, a Helmoltz equation for u?_{is formulated while considering the contributions}

to the original momentum equation (4.2) from all the various terms. ∇2u?− 2u ? ν0∆t = −2u n ν0∆t − 1 ν0 3Hn_{− H}n−1₋ 2 ν0 1 ρn+1∇. 2µ n+1_D(un_{) + ∇}2_un (4.16) Where Hn = F_advn + F_capin + F_strn .

This system is solved effectively with an FFT-based elliptic solver (Costa , 2018) which is also used here.

4.2 Higher Order Spatial Discretization for the Stress Term

The weighted essentially non-oscillatory (WENO) schemes for spatial discretization are dis-cussed here. To study these schemes, first consider a simple scalar hyperbolic equation in conservative form (Xing & Shu, 2005).

ut+ f (u)x= 0 (4.17)

One can use a conservative flux to discretize f (u)xas follow on a staggered grid:

f (u)x = ˆ f_i+1 2 − ˆf_i−1 2 ∆x (4.18) Where ˆf_i+1 2 and ˆfi− 1

2 are numerical fluxes and are computed through the neighbouring point

values. Considering fj = f (uj), the detailed steps one needs to follow to calculate (2k − 1)th

order WENO scheme is presented in this section.

At the first step, k numerical fluxes (of the order k) should be calculated using k different stencils (Sr(i), where r = 0, 1, 2, . . . , k − 1).

ˆ f(r) i+1₂ = k−1 X j=0 crjfi−r+j (4.19)

Considering the fifth order WENO scheme, which is used in the present study, means that k = 3. Therefore, one can write three different third order fluxes as below:

ˆ f(0) i+1₂ = 1 3fi+ 5 6fi+1− 1 6fi+2 ˆ f(1) i+1₂ = − 1 6fi−1+ 5 6fi+ 1 3fi+1 ˆ f_i+(2)1 2 = 1 3fi−2− 7 6fi−1+ 11 6 fi (4.20)

Finally, the (2k − 1)th order flux of WENO is the combination of all the above mentioned fluxes. ˆ f(r) i+1₂ = k−1 X j=0 wrf_i+1 2 (4.21)

Where wr represents the non-linear weights and can be calculated through the following

(28)

CHAPTER 4. COMPUTATIONAL METHODOLOGY 23 wr= αr Pk−1 s=0αs αr= dr ( + βr)2 (4.22)

Before introducing the parameters in the relations, it is important to notice that wr ≥ 0 and

Pk−1

j=0wr = 1. In 4.22, r is the linear weight, βr is the smoothness indicator of the stencil

Sr(i), and is a small constant (typically 10−6). Simply, dr which makes the scheme to be

of the (2k − 1)th _{order of accuracy, β}

r considers the smoothness of the flux function, and

guarantees the non-zero denominator.

For the fifth order WENO scheme, the linear weights and the smoothness indicators can be written as (Xing & Shu, 2005):

d0 = 3 10, d1= 3 5, d2 = 1 10 β0 = 13 12(fi− 2fi+1+ fi+2) 2₊ 1 4(3fi− 4fi+1+ fi+2) 2 β1 = 13 12(fi−1− 2fi+ fi+1) 2₊1 4(fi−1− fi+1) 2 β2 = 13 12(fi−2− 2fi−1+ fi) 2₊ 1 4(fi−2− 4fi−1+ 3fi) 2 (4.23)

4.3 Algorithmic Procedure

In the previous sections, we furnished the reader with notes on the schemes used and other necessary numerical details (in a somewhat isolated fashion). We now explain how the whole code is combined together and implemented within the framework of a computer code. The programming language used is FORTRAN 90 though the algorithmic implementation is not program dependent.

Step

1

Impose the condition at the initial time t = t0

                   un=0= uinit= 0 x ∈ Ω τn=0= τinit = 0 x ∈ Ω φn=0= 1 x ∈ Ωf l φn=0= 0 x ∈ Ωex pn=0= 0 x ∈ Ω (4.24)

For the boundary condition, we use a periodic boundary condition (Hoyle & Fielding , 2016a) except for the top and bottom planes of the rectangular computational domain in which the Neumann boundary condition is applied for the pressure and volume-of-fluid parameter i.e. ∂p_∂~_n = 0 and _∂~∂φ_n = 0; while the no slip, no-penetration is employed for the velocity i.e. un_{= 0.}

Step

2 φn→ φn+1

The MTHINC VOF relationship is solved to obtain φn+1_{with the appropriate inputs from}

(29)

function (H◦)n at the current time step. This is done while making reference to equations (4.10), (4.11),(4.12) and (4.13).

Step

3 µn_i,j,k → µn+1_i,j,k and ρn i,j,k→ ρ

n+1 i,j,k

Update the viscosity µn+1_i,j,k and density ρn+1_i,j,kfields via equation (4.14).

Step

4 τn → τn+1

Integrate the non-Newtonian constitutive stress relationships - equations (4.6) - (4.9) are integrated with the second order Adams-Bashforth scheme.

Step

5 un→ un+12

The momentum equation is advanced by an Adams-Bashforth scheme. To achieve this, initially consider all the terms in equations (4.1) - (4.3) excluding the diffusion term Fn

dif f.

These terms are summed up for the present and previous time steps -Hn = F_advn + F_capin + F_strn+1

Hn−1 = F_advn−1+ F_capin−1+ F_strn (4.25) Thereafter, through an implicit decomposition of the diffusion term, the Helmholtz equa-tion (4.16) is computed with an elliptic solver to obtain un+12.

Step

6 un+12 → un+1and pn→ pn+1

The solenoidal condition on the velocity field is enforced. This is actualized by solving a Poisson equation for pressure (equation (4.5)) with an FFT-based elliptic solver (Costa , 2018). The velocity field is corrected with this newly obtained pressure (equation (4.4)). Impose the boundary conditions as in 1 and repeat steps 2 - 6 until final necking occurs.

These steps are repeated until final necking occurs in the non-Newtonian filament. Differ-ent criteria have been presDiffer-ented to determine when the EVP has "necked". For a computational-based study as this (as opposed to an experimental study), the EVP is assumed to have necked when the minimum necking width along the EVP filament hEV P,min is less than a

specified criteria i.e. hEV P,min < ψ, where the necking criteria is ψ. Definitely, ψ would be

set by the (minimum) mesh size of the simulation. An estimation of ψ is a non dimensional mesh size i.e. ˜∆x = ∆x_L where L is a characteristic length scale. For most simulations run, the mesh size used2is ∆x = 1/32. Hence, ψ ≈ ˜∆x = ∆x_L = 1/32₂ = 0.02. This is about two times finer than the criteria set in (Balmforth et al. , 2010) which is hmin < ψ = 0.05for halting the

numerical simulation in their simulation of a viscoplastic fluid.

(30)

Chapter 5

Results

In this chapter, the results from the simulation is compared with previous related studies on instability in extensional filament stretching of non-Newtonian fluids. Specifically, the analytic results from (Saramito , 2007) for homogeneous, uniaxial extension; and the numer-ical simulation of instability in viscoelastic fluids from (Hoyle & Fielding , 2016a) are quite insightful for comparison.

5.1 Dimensionless Scaling

The initial length, L0of the EVP and the inverse of the imposed extension rate, 1/ ˙0are

cho-sen as suitable scales to reprecho-sent the characteristic length (Hoyle & Fielding , 2016a; Szabo , 1997; Hoyle & Fielding , 2016b) and time scale (Szabo , 1997; Yao & McKingley, 1998) respec-tively. It is not uncommon to observe that the length and time are scaled with the radius (Yao & McKingley, 1998; Mckinley & Sridhar, 2002) and relaxation time (λ) respectively (Hoyle & Fielding , 2016a; Mckinley & Hassager, 1999; Hoyle & Fielding , 2016b). Since our model is based on a quasi-3D rectangular idealization of the EVP, the initial length (L0) seems more

appropriate to represent the length scale. The extension rate ( ˙0) of the EVP primary dictates

the dynamics of the flow. This justifies why the extension rate is used to scale the time. By inference, the characteristic velocity is L0˙0 and stress is scaled with µ0˙0 where µ0 1 is the

dynamic viscosity of the EVP.

Following the definition of the characteristic scales, some important dimensional parame-ters which will motivate the result presentation are also clarified below (Mckinley & Sridhar, 2002; Szabo , 1997; Saramito , 2007). The Weissenberg number2(W e), Reynolds number (Re), Bingham number (Bi) and Capillary number (Ca).

W e = λ ˙0 Re = ρ0˙0L20 µ0 Bi = τ0 µ0˙0 Ca = µ0˙0L0 σ (5.1)

In all our simulations, the Re was kept at ≈ 0.01. This is consistent with the experimental procedure for extensional rheometry (Nieuwkoop & Czernicki, 1996; Yao & McKingley, 1998;

1_{Note that µ}

0= µs+ µpwhere µsand µpare the solvent and polymer viscosities of the EVP respectively 2_{There exist a subtle difference in terminological usage between the Deborah number, De and the Weissenberg}

number, W e. The use intended here is with regards to the ratio between the elastic force to the viscous force (Poole , 2012). The papers from (McKinley & Sridhar , 2002; Yao & McKingley, 1998) use De to represent this ratio. It is however not uncommon to find studies (Saramito , 2007,0) that represent this ratio with W e

(31)

26 CHAPTER 5. RESULTS

McKinley & Sridhar , 2002; Sridhar et al., 1991; Tirtaatmadja & Sridhar , 1993; Spiegelberg et al., 1991; Rasmussen & Hassager, 1996).

We report the stress values at the midsection of the filament3. This is generally acceptable in the literature (Hoyle & Fielding , 2016a; Yao & McKingley, 1998; Szabo , 1997; Mckinley & Sridhar, 2002) as representative for reporting the stress since the stress field is heterogeneous along the length of the filament. However, the strain rate reported is the globally imposed strain rate (sometimes referred to as ¯˙) and not the local strain rate obtained at the midsection. In the results that follow, four failure strains are important and reported, following earlier discussions. The strain at final failure (extension at strain-to-break), the strain from the Con-sidère force criterion, the strain from the stress curvature criterion and the strain at maximum stress designated as ST B, F, σ and σ,max respectively. When these indices are denoted

with an overline i.e. ST B, F, σand σ,max- it means they are obtained from the simulations

(and not the uniaxial, homogeneous assumption). Practically, extension at strain-to-break (ST B) is reported only from the simulations, as the homogeneous solution never predicts

final failure.

One final note to mention is that because the primary focus of this work is in non-Newtonian fluids, the retardation parameter (Saramito , 2007) - defined as α = µp/µ0where

µ0 = µs+ µp; and µsand µpare the solvent and polymer viscosities of the EVP respectively

-is kept high in the simulations (0.9). Th-is implies that bulk elastic stresses dominate and are more significant than viscous stresses.

5.2 Instability in Viscoelastic Fluids

Figure 5.1: Plot of the temporal evolution of the transient extensional viscosity µ+_E and non dimensional axial force F with extension. The blue broken line is the result from the present simulation, the magenta broken line is the result from (Hoyle & Fielding , 2016a) and the red broken line is the analytic results assuming homogeneous extension. The input parameters are ˙ = 0.1, λ = 1, Re = 0.01, Bi = 0, Ca = ∞.

(32)

CHAPTER 5. RESULTS 27

Figure 5.2: Plot of the temporal evolution of the transient extensional viscosity µ+_E and non dimensional axial force F with extension. The blue broken line is the result from the present simulation, the magenta broken line is the result from (Hoyle & Fielding , 2016a) and the red broken line is the analytic results assuming homogeneous extension. The input parameters are ˙ = 1, λ = 1, Re = 0.01, Bi = 0, Ca = ∞.

Figure 5.3: Plot of the temporal evolution of the transient extensional viscosity µ+_E and non dimensional axial force F with extension. The blue broken line is the result from the present simulation, the magenta broken line is the result from (Hoyle & Fielding , 2016a) and the red broken line is the analytic results assuming homogeneous extension. The input parameters are ˙ = 10, λ = 1, Re = 0.01, Bi = 0, Ca = ∞.

Extensional Instability in Complex Fluids: A Computational Study