Turbulent flow non-Newtonian pipe flow

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.optimize as spo
from scipy.optimize import fsolve

from rheoflow import friction_factor_property as f

import warnings
warnings.filterwarnings("ignore")

Technical background - friction factor calculations for non-Newtonian fluids

Von Karmen equation for f

\(\frac{1}{\sqrt{f}} = 4.0 log_{10} \left(Re \sqrt{f} \right) - 0.4\)

The function f_vonkarmen returns the value of the equation given \(f_{guess}\) and Re

\(\frac{1}{\sqrt{f_{guess}}} - 4.0 log_{10} \left(Re \sqrt{f_{guess}} \right) + 0.4\)

The scipy function brentq is used here to solve the equation above for f by setting the equation to 0. brentq requires three arguments: 1. A function with one argument that returns value of equation to be solved 2. Minimum value of variable to be solved for 3. Maximum value of variable top be solved for

The function f_vonkarmen has two arguments, \(f_{guess}\) and Re. A lambda function is used here to create a function with a single argument \(f_{guess}\) for a specified value of Re.

lambda \(f_{guess}\): f_vonkarmen(\(f_{guess}\),Re)

Non-Newtonian Von Karmen equation (Dodge and Metzner)

The modified Von Karmen equation from Dodge and Metzner is

\(\frac{1}{\sqrt{f}} = \frac{4}{n'^{0.75}} Log \left( f^{1-\frac{n'}{2}} Re \right) - \frac{0.4}{n'^{1.2}}\)

This equation is solved in the same manner as the Newtonian version with scipy.brentq.

Calculation of n’

In the series of Metzner papers, n’ is defined as

\(n' = \frac{d Log \left( \tau_w \right)}{d Log \left( \dot{\gamma}_a \right) }\)

where \(\tau_w\) is the wall shear stress and \(\dot{\gamma}_a\) is the apparent wall shear rate.

\(\tau_w = \frac{D}{4} \frac{\Delta P}{L}\)

\(\dot{\gamma}_a = \frac{8U}{D}\)

To make this calculation model independent, it is computed numerically now. In the future automatic differentiation would be nice. For now:

\(n' = \frac{Log (\eta (x+\Delta x)(x+\Delta x)) - Log (\eta(x-\Delta x )(x-\Delta x))}{Log(x+\Delta x) - Log(x-\Delta x)}\)

where x = \(\dot{\gamma}_w\)

Equations to solve

\(\frac{1}{\sqrt{f}} = \frac{4}{n'^{0.75}} Log \left( f^{1-\frac{n'}{2}} Re \right) - \frac{0.4}{n'^{1.2}}\)

\(Re = \frac{\rho D U}{\eta(\dot{\gamma}_w)}\)

\(\tau_w = \eta(\dot{\gamma}_w) \dot{\gamma}_w\)

\(\tau_w = \frac{D}{4} \frac{\Delta P}{L}\)

\(f = \frac{\Delta P D}{2 \rho^2 L}\)

with unknowns f, Re, \(\tau_w\), \(\dot{\gamma}_w\), \(\Delta P\) assuming \(\eta(\dot{\gamma})\), \(\rho\), D, and L are specified

Note on friction factor definitions

There are at least three friction factor definitions: Darcy, Fanning, and Newton. The Fanning friction factor (used here) is

\(f = \frac{D \Delta P}{2 \rho U^2 L}\)

The resulting laminar flow equation is

\(f = \frac{16}{Re}\)

The Darcy friction factor is

\(f_D = \frac{ 2 D \Delta P}{\rho U^2 L}\)

and the resulting laminar flow equation is

\(f_D = \frac{64}{Re}\)

Finally there is the Newton number. We don’t use that here.

Example power-law viscosity model calculations

Define a viscosity model

def power_law(k,n,gammadot):
    return k*(gammadot+1.e-9)**(n-1.)

k = .1
n = .5
viscosity = lambda x: power_law(k,n,x)

Create (instantiate) an object of class friction_factor

For this case the pipe diamter is 0.09387m, the pipe length is 100m, and the fluid density is 1000 kg/m^3.

a = f.friction_factor(name='test',rho=1000.,d=.09738,l=100.,viscosity=viscosity)

Case 1 - mass flow rate = 200 kg/min

The mass flow rate is 200 kg/min. The average velocity, u, is calculated anf the attribute a.u is set to this.

mdot = 200.0 # kg/min
density = 1000. # kg/m^3
q = mdot/60.0/density
u = q / (3.14159*(a.d/2.0)**2)

a.u=u
print(a)

Name= test
Diameter = 0.09738
Length = 100.0
Density = 1000.0
U = 0.44755837705026
Pressure drop = 2958.7020853480603
Friction factor = 0.007191866741121529
Reynolds number = 3139.283065162426
Wall shear rate = 51.882636029577945
Wall shear stress = 0.7202960226779853

Case 2 - Pressure drop is 2958.7 Pa

a.pressure_drop = 2958.7
print(a)

Name= test
Diameter = 0.09738
Length = 100.0
Density = 1000.0
U = 0.4475581150143437
Pressure drop = 2958.7
Friction factor = 0.007191870093534452
Reynolds number = 3139.279014555334
Wall shear rate = 51.88256289391151
Wall shear stress = 0.720295515

Accessing each attribute

a.u

0.4475581150143437

a.pressure_drop

2958.7

a.f

0.007191870093534452

a.d

0.09738

a.l

100.0