# 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. Nov 26, 2019 at 21:28
• Crossposted here. @MarcoB The crosspost has an attached notebook. Nov 26, 2019 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. Nov 26, 2019 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. Nov 29, 2019 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. Nov 29, 2019 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. Nov 27, 2019 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. Nov 27, 2019 at 12:51
• Can you tell me why the NDSolve could not do it without the false method? Nov 27, 2019 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/… Nov 27, 2019 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. Nov 29, 2019 at 3:36
• I am glad my post attracted some interests :-) Nov 29, 2019 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 Nov 29, 2019 at 6:24
• @TimLaska Thanks a lot. It really helps. Nov 29, 2019 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). Nov 29, 2019 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. Nov 29, 2019 at 21:58
• @TimLaska Thanks. It really helps. Dec 2, 2019 at 18:40