I'm attempting to take the Navier Stokes Equation and coming up with an expression that will allow me to numerically determine the velocity of non-Newtonian fluid flow. The text I'm using is Cengel and Cimbala (Chpt. 9) which is available for download HERE.
I would really appreciate it if others here can help out by reviewing my maths and assisting with sanity checking.
The textbook states that the viscous stress tensor for an incompressible Newtonian fluid with constant properties is as follows:
$$\tau_{ij}=2\mu\varepsilon_{ij}$$
...where $\mu$ is the viscosity coefficient and $\varepsilon_{ij}$ is the shear strain rate tensor...
$$ \varepsilon_{ij}=\begin{bmatrix} \varepsilon_{xx}&\varepsilon_{xy}&\varepsilon_{xz}\\\varepsilon_{yx}&\varepsilon_{yy}&\varepsilon_{yz}\\\varepsilon_{zx}&\varepsilon_{zy}&\varepsilon_{zz} \end{bmatrix} =\begin{bmatrix} \frac{\partial u}{\partial x}&\frac{1}{2}(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})&\frac{1}{2}(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x})\\\frac{1}{2}(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y})&\frac{\partial v}{\partial y}&\frac{1}{2}(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})\\\frac{1}{2}(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z})&\frac{1}{2}(\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z})&\frac{\partial w}{\partial z} \end{bmatrix} $$
I'm modelling for an industrial oil-based lubricant, used in power transformers, which is a shear thinning fluid. I'm arguing that, for this fluid, the viscous stress tensor becomes:
$$\tau_{ij}=2\mu\sqrt{\varepsilon_{ij}}$$
When the fluid is in motion, the stress tensor is as follows:
$$\sigma_{ij}=\tau_{ij}+\lambda_{ij}= \begin{bmatrix} 2\mu\sqrt{\frac{\partial u}{\partial x}}&\sqrt{2}\mu\sqrt{\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}}&\sqrt{2}\mu\sqrt{\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}}\\\sqrt{2}\mu\sqrt{\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}}&2\mu\sqrt{\frac{\partial v}{\partial y}}&\sqrt{2}\mu\sqrt{\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}}\\\sqrt{2}\mu\sqrt{\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}}&\sqrt{2}\mu\sqrt{\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}}&2\mu\sqrt{\frac{\partial v}{\partial y}} \end{bmatrix} + \begin{bmatrix} -P&0&0\\0&-P&0\\0&0&-P \end{bmatrix}$$
where $\lambda$ is the hydrostatic stress term and P is the hydrostatic/thermodynamic pressure actually inwardly and orthogonally on the unit cell (control volume).
Next, I need to selectively plug these shear stresses into Cauchy's equation. The three Cartesian components of Cauchy equation, as stipulated in the text, is as follows:
x-component:
$$\rho\frac{Du}{Dt}=\rho g_{x} + \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{yx}}{\partial y} + \frac{\partial \sigma_{zx}}{\partial z}$$
y-component:
$$\rho\frac{Dv}{Dt}=\rho g_{y} + \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y} + \frac{\partial \sigma_{zy}}{\partial z}$$
z-component:
$$\rho\frac{Dw}{Dt}=\rho g_{z} + \frac{\partial \sigma_{xz}}{\partial x} + \frac{\partial \sigma_{yz}}{\partial y} + \frac{\partial \sigma_{zz}}{\partial z}$$
Pit-stop question: The term $\frac{Du}{Dt}$ is acceleration in the x direction and since $\rho$ is mass per volm, then $\rho\frac{Dw}{Dt}$ should be force in the x direction per unit volume?
Anyhow, after making the relevant substitutions...
x-component:
$$\rho\frac{Du}{Dt}=\rho g_{x} + \frac{\partial}{\partial x}(2\mu\sqrt{\frac{\partial u}{\partial x}}-P) + \frac{\partial}{\partial y}(\sqrt{2}\mu\sqrt{\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}}) + \frac{\partial}{\partial z}(\sqrt{2}\mu\sqrt{\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}})$$
y-component:
$$\rho\frac{Dv}{Dt}=\rho g_{y} + \frac{\partial}{\partial x}(\sqrt{2}\mu\sqrt{\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}}) + \frac{\partial}{\partial y}(2\mu\sqrt{\frac{\partial v}{\partial y}}-P) + \frac{\partial}{\partial z}(\sqrt{2}\mu\sqrt{\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}})$$
z-component:
$$\rho\frac{Dw}{Dt}=\rho g_{z} + \frac{\partial}{\partial x}(\sqrt{2}\mu\sqrt{\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}}) + \frac{\partial}{\partial y}(\sqrt{2}\mu\sqrt{\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}}) + \frac{\partial}{\partial z}(2\mu\sqrt{\frac{\partial v}{\partial y}}-P)$$
And, of course, this is where things take an ugly turn. That pesky root prevents the three Cauchy equations from converging neatly into the familiar NS we all know. Any suggestions as to where to go where here?