# Non-Newtonian Momentum Eqn in a rectangular pipe

I am trying to (numerically) solve the momentum eqn for a non-Newtonian fluid in a pipe with rectangular cross section using Mathematica.

Here are the assumptions:

So, the flow is only in z direction (axial direction) and varies only in x and y directions. Sides of the rectangle are 2a and 2b in x and y directions, respectively.

The z-momentum equation is:

The boundary conditions are:

The non-Newtonian model I am using is called Carreau model, in which the fluid viscosity is a function of scalar shear rate as:

All the parameters in above model are constants, except for gamma dot, which is the scalar shear rate, evaluated as:

Update: This scalar strain rate is off, please see the answer by @Tim Laska for the corrected version.

So, I am using NDSolve and NDSolveValue to solve it, but had no luck (error: The spatial derivative order of the PDE may not exceed two).

I am attaching the notebook I put together so far in this link (https://community.wolfram.com/groups/-/m/t/1832082). Can you please take a look at it and let me know how I can work around it.

Important: Note that in the notebook, I am actually trying to solve the eqn only for 1/4 of the geometry due to the symmetry (0<=x<=a and 0<=y<=b).

Edit: The dUz/dx and dUz/dy should actually be partial derivatives in my description above.

Thanks

• It would be preferable if you actually copy/pasted the relevant code with all definitions in your question instead (see: How to copy code from Mathematica so it looks good on this site, and also what we mean by a minimal working example). I am wary of external files in general, and particularly of ones that contain executable code of unknown origin, so I prefer not to download and open your file. – MarcoB Nov 26 '19 at 21:28
• Crossposted here. @MarcoB The crosspost has an attached notebook. – Rohit Namjoshi Nov 26 '19 at 21:51
• The link is just a crosspost from wolfram forum with an attached mathematica notebook. I am assuming it is pretty safe for you to open it. See the comment by Rohit Namjoshi above. – Muhamad Mohaqeq Nov 26 '19 at 22:36
• You may want to check your definition of the scalar shear rate. I added an answer that I think corrects it, but the other answers did the heavy lifting. – Tim Laska Nov 29 '19 at 3:39
• @TimLaska Yes, my definition of scalar strain rate was off. Thanks for pointing that out. I Updated my questions and mentioned it. – Muhamad Mohaqeq Nov 29 '19 at 16:59

We can use the FEM solver from here Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM or from here https://community.wolfram.com/groups/-/m/t/1433064

Let's make a small code modification

a = 0.005  (*m ; range in x direction*);
b = 0.005 (*m ; range in y direction*);
\[Rho] = 1060;(*kg/m^3*);
\[Mu]N = 0.0035;(*Pa.s ; Newtonian viscosity*)
Q = 0.00001; (*m^3/s ; flow rate*)
dPdzN = -10 ; (*Pa/m ; Newtonian pressure gradient*)
dPdzNN = -25; (*Pa/m ; non-Newtonian pressure gradient*)
n = 0.3568; (*exponent in Carreau viscosity model*)
\[Mu]0 = 0.056 ; (*Pa.s ; 0-shear viscosity*)
\[Mu]\[Infinity] = 0.0035; (*Pa.s ; infinite-shear viscosity*)
\[Lambda] = 3.313;(*s ; coefficient in Carreau model*)

\[Mu]eff = \[Mu]\[Infinity] +(\[Mu]0 - \[Mu]\[Infinity])*(1 + (\
\[Lambda]*\[Gamma]dot)^2)^((n - 1)/2); (*Carreau model*)
\[Gamma]dot =
0.5*Sqrt[0.5*((D[Uz[t - t0][x, y], x])^2 + (D[Uz[t - t0][x, y],
y])^2)]; (*Scalar shear rate*)

eqn =
D[\[Mu]eff*(D[u[x, y], x]), x] + D[\[Mu]eff*(D[u[x, y], y]), y] -
dPdzNN ;
bc = DirichletCondition[u[x, y] == 0, True];
Uz[0][x_, y_] := 0

t0 = 1/10; nn = 10; reg =
ImplicitRegion[-a <= x <= a && -b <= y <= b, {x, y}]; Do[
Uz[t] = NDSolveValue[{eqn == (u[x, y] - Uz[t - t0][x, y])/t0, bc},
u, {x, y} \[Element] reg,
Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}];, {t, t0,
nn t0, t0}]


The numerical solution converges quickly, so we can use nn=10

Table[ContourPlot[Uz[t][x, y], {x, y} \[Element] reg, Contours -> 20,
PlotLegends -> Automatic, PlotLabel -> Row[{"t = ", t}],
ColorFunction -> "Rainbow"], {t, t0, nn t0, 3 t0}]

ListPlot[Table[Uz[t][0, 0], {t, 0, nn t0, t0}]]


• Thanks Alex, I can try it when I get back to my laptop. However, can you provide more explanation/details on the modifications you did. – Muhamad Mohaqeq Nov 27 '19 at 12:24
• The method of the false transient with step t0 is used there. We calculate \[Mu]eff in the previous step t-t0. With a small step t0<1/50, this model describes the formation of a flow in a channel. Then t is the time. – Alex Trounev Nov 27 '19 at 12:51
• Can you tell me why the NDSolve could not do it without the false method? – Muhamad Mohaqeq Nov 27 '19 at 17:22
• This equation can be solved without using the method of the false transient with the low-level FEM tools. See mathematica.stackexchange.com/questions/202446/… – Alex Trounev Nov 27 '19 at 18:27

Here is a much simpler way to solve this:

a = 0.005;
b = 0.005;
reg = Rectangle[{0,0}, {a, b}];
dPdzN = -10; dPdzNN = -25; n = 0.3568; mu0 = 0.056; muinf = 0.0035; lam = 3.313;


(* gdot = 1/2* Sqrt[1/2*((Derivative[1, 0][u][x, y])^2 + (Derivative[0, 1][u][x,y])^2)]; *)
(* corrected according to Tim Laska *)
gdot = Sqrt[(Derivative[0, 1][u][x, y]^2 + Derivative[1, 0][u][x, y]^2)];
mueff = muinf + (mu0 - muinf)*(1 + (lam*gdot)^2)^((n - 1)/2);
eqn = Inactive[Div][mueff Inactive[Grad][u[x, y], {x, y}], {x, y}];


Look at this equation: It reads like what you have written.

sol = NDSolveValue[{eqn == dPdzNN,
DirichletCondition[u[x, y] == 0, x == a || y == b]}, u, {x, y}\[Element] reg,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}];


This will solve the PDE much, much faster that what Alex has shown. Visualize the result.

ContourPlot[sol[x, y], {x, y} \[Element] reg, Contours -> 20,
PlotLegends -> Automatic, ColorFunction -> "Rainbow"]


sol[0, 0]
0.011507510331532756

• I think the shear rate definition is off in the OP. I added an answer that I think corrects the definition, but I mostly used your code. I'm kind of excited to have a MMA implementation of a non-Newtonian fluid. – Tim Laska Nov 29 '19 at 3:36
• I am glad my post attracted some interests :-) – Muhamad Mohaqeq Nov 29 '19 at 13:53

It appears that the shear rate in the OP did not use a symmetrized strain rate tensor and that lowered the actual shear rate used by the Carreau model, thereby increasing the apparent viscosity. The symmetrized strain rate tensor is given by:

$${\Delta _{ij}} = \left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}} \right)$$

The shear rate is expressed in terms of the scalar tensor product of the symmetrized strain rate tensor.

$$\dot \gamma = \sqrt {\frac{1}{2}{\mathbf{\Delta :\Delta }}} = \sqrt {\frac{1}{2}\sum\limits_i {\sum\limits_j {{\Delta _{ij}}{\Delta _{ji}}} } }$$

We can use Mathematica to estimate the term under the radical

crds = Array[x, 3];
vels = Array[u @@ crds, 3] /. {u[x[1], x[2], x[3]][1] -> 0,
u[x[1], x[2], x[3]][2] -> 0};
Del = {D[#, x[1]], D[#, x[2]], D[#, x[3]]} &;
g = (Del[#] & /@ vels) /. D[u[x[1], x[2], x[3]][3], x[3]] -> 0;
gsym = g + Transpose@g;
1/2 Sum[gsym[[i, j]]*gsym[[j, i]], {i, 1, 3}, {j, 1, 3}] // Expand


$$\left( {{{\left( {\frac{{\partial u}}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial u}}{{\partial y}}} \right)}^2}} \right)$$

So, the corrected shear rate should look like:

$$\dot \gamma = \sqrt {\left( {{{\left( {\frac{{\partial u}}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial u}}{{\partial y}}} \right)}^2}} \right)}$$

Now, we can plug the new shear rate definition into @user21's answer and obtain:

a = 0.005;
b = 0.005;
reg = Rectangle[{-a, -b}, {a, b}];
dPdzN = -10; dPdzNN = -25; n = 0.3568; mu0 = 0.056; muinf = 0.0035; \
lam = 3.313;
gdot = Sqrt[(Derivative[0, 1][u][x, y]^2 +
Derivative[1, 0][u][x, y]^2)];
mueff = muinf + (mu0 - muinf)*(1 + (lam*gdot)^2)^((n - 1)/2);
eqn = Inactive[Div][mueff Inactive[Grad][u[x, y], {x, y}], {x, y}];
sol = NDSolveValue[{eqn == dPdzNN,
DirichletCondition[u[x, y] == 0, True]}, u, {x, y} \[Element] reg,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}];
ContourPlot[sol[x, y], {x, y} \[Element] reg, Contours -> 20,
PlotLegends -> Automatic, ColorFunction -> "Rainbow"]
Plot[{sol[x, 0], sol[0, x]}, {x, -a, a}]
Plot[{sol[x, x]}, {x, -a, a}]
sol[0, 0] (* 0.0115075 *)


The maximum velocity of the solution with the corrected shear rate compares favorably to the solutions of the commercial CFD solvers COMSOL (see below) and AcuSolve (0.01145).

# Update: Quarter Symmetry Case for a Cylindrical Pipe

As stated in the comments, since the solve times are short, it maybe easier to apply quarter symmetry to a disk versus trying to develop an axisymmetric formulation. The following is quarter symmetry case of a disk applied to a thicker and more shear thinning fluid.

a = 0.005;
reg = Disk[{0, 0}, a, {0, Pi/2}];
dPdzN = -10; dPdzNN = -1000; n = 0.03568; mu0 = 5.6; muinf = 0.0035; \
lam = 3.313;
gdot = Sqrt[(Derivative[0, 1][u][x, y]^2 +
Derivative[1, 0][u][x, y]^2)];
mueff = muinf + (mu0 - muinf)*(1 + (lam*gdot)^2)^((n - 1)/2);
eqn = Inactive[Div][mueff Inactive[Grad][u[x, y], {x, y}], {x, y}];
sol = NDSolveValue[{eqn == dPdzNN,
DirichletCondition[u[x, y] == 0, x^2 + y^2 == a^2]},
u, {x, y} \[Element] reg,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}];
velavg = NIntegrate[2 Pi x sol[x, 0], {x, 0, a}]/(\[Pi] a^2);
ContourPlot[sol[x, y], {x, y} \[Element] reg, Contours -> 20,
PlotLegends -> Automatic, ColorFunction -> "Rainbow"]
Plot[{sol[x, 0], sol[0, x]}, {x, 0, a}, PlotRange -> All]
Plot[{sol[x, 0], 2 velavg (1 - (x/a)^2)}, {x, 0, a},
PlotLegends -> {"NN", "Newt"}]
sol[0, 0]


The solution has close agreement with COMSOL's axisymmetric formulation, but only required that we change the region from rectangle to disk and adjusted the Dirichlet condition to follow the arc.

# Update 2: Axisymmetric Approximation Using a $$5^\circ$$ Wedge

Some CFD solvers (e.g., openFOAM and AcuSolve) only have a 3D Cartesian formulation. To model an "axisymmetric" case, one usually just performs the simulation on a $$5^\circ$$ wedge and applies the appropriate symmetry boundary conditions. I tried that approach with Mathematica and it turned out to be quite fast compared to the quarter symmetry case.

I like to view the computational mesh, so I imported the FEM package.

Needs["NDSolveFEM"]


Here is the code to build the mesh and do some post processing on the solution:

a = 0.005;
reg = Disk[{0, 0}, a, {0, Pi/72}];
(mesh = ToElementMesh[reg])["Wireframe"]
dPdzN = -10; dPdzNN = -1000; n = 0.03568; mu0 = 5.6; muinf = 0.0035; \
lam = 3.313;
gdot = Sqrt[(Derivative[0, 1][u][x, y]^2 +
Derivative[1, 0][u][x, y]^2)];
mueff = muinf + (mu0 - muinf)*(1 + (lam*gdot)^2)^((n - 1)/2);
eqn = Inactive[Div][mueff Inactive[Grad][u[x, y], {x, y}], {x, y}];
sol = NDSolveValue[{eqn == dPdzNN,
DirichletCondition[u[x, y] == 0, x^2 + y^2 == a^2]},
u, {x, y} \[Element] mesh];
velavg = NIntegrate[2 Pi x sol[x, 0], {x, 0, a}]/(\[Pi] a^2);
Plot[{sol[x, 0], 2 velavg (1 - (x/a)^2)}, {x, 0, a},
PlotLegends -> {"NN", "Newt"}]
ParametricPlot[a t { Cos[th], Sin[th]}, {t, 0, 1}, {th, 0, 2 Pi},
Mesh -> 10, MeshFunctions -> {sol[a #3, 0] &},
ColorFunction -> (ColorData["Rainbow"][sol[a #3, 0]/sol[0, 0]] &),
Axes -> {False, False, True}]
RevolutionPlot3D[sol[a t, 0], {t, 0, 1},
ColorFunction -> (ColorData["Rainbow"][sol[a #4, 0]/sol[0, 0]] &),
PlotRange -> All]
sol[0, 0]
`

• Great. Now we just need gdot for the Navier-Stokes equation ;-) and a whole new world opens – user21 Nov 29 '19 at 6:24
• @TimLaska Thanks a lot. It really helps. – Muhamad Mohaqeq Nov 29 '19 at 13:51
• @TimLaska Thanks for correcting the scalar strain rate. We are solving the same non-Newtonian flow in a circular cross section, which is a simpler version of it, because use of the cylindrical coordinate leaves us only with one variable (r). So we have only delUzdelr as non-zero terms (z being the axial direction). However, I am using delUzdelr (dUz/dr) as the scalar strain rate. I ma not sure if that is correct? Cab you verify that?! (Since it is in cylindrical coordinated, I am still not sure). – Muhamad Mohaqeq Nov 29 '19 at 13:58
• @MuhamadMohaqeq In my opinion, unless it is already formulated, trying to come up with your own axisymmetric implementation often is more trouble than it is worth (1/r blows up at r=0). In your current quarter sym case, your solve times are quite short. Why not change the domain to a quarter disk and then you do not need to mess with the shear rate? Also, you could use the quarter disk solution to verify your axisymmetric formulation. – Tim Laska Nov 29 '19 at 21:58
• @TimLaska Thanks. It really helps. – Muhamad Mohaqeq Dec 2 '19 at 18:40