2
$\begingroup$

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:

enter image description here

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:

enter image description here

The boundary conditions are:

enter image description here

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

enter image description here

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

enter image description here

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

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – MarcoB Nov 26 '19 at 21:28
  • $\begingroup$ Crossposted here. @MarcoB The crosspost has an attached notebook. $\endgroup$ – Rohit Namjoshi Nov 26 '19 at 21:51
  • $\begingroup$ 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. $\endgroup$ – Muhamad Mohaqeq Nov 26 '19 at 22:36
  • 1
    $\begingroup$ 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. $\endgroup$ – Tim Laska Nov 29 '19 at 3:39
  • $\begingroup$ @TimLaska Yes, my definition of scalar strain rate was off. Thanks for pointing that out. I Updated my questions and mentioned it. $\endgroup$ – Muhamad Mohaqeq Nov 29 '19 at 16:59
6
$\begingroup$

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}]]

Figure 1

|improve this answer|||||
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Muhamad Mohaqeq Nov 27 '19 at 12:24
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev Nov 27 '19 at 12:51
  • $\begingroup$ Can you tell me why the NDSolve could not do it without the false method? $\endgroup$ – Muhamad Mohaqeq Nov 27 '19 at 17:22
  • $\begingroup$ 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/… $\endgroup$ – Alex Trounev Nov 27 '19 at 18:27
4
$\begingroup$

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;

Specify your equations:

(* 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"]

enter image description here

sol[0, 0]
0.011507510331532756`
|improve this answer|||||
$\endgroup$
  • 2
    $\begingroup$ 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. $\endgroup$ – Tim Laska Nov 29 '19 at 3:36
  • $\begingroup$ I am glad my post attracted some interests :-) $\endgroup$ – Muhamad Mohaqeq Nov 29 '19 at 13:53
4
$\begingroup$

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).

COMSOL solution

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]

Quarter Symmetry Disk

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. COMSOL Axisymmetric solution

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["NDSolve`FEM`"]

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]

5 degree wedge

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Great. Now we just need gdot for the Navier-Stokes equation ;-) and a whole new world opens $\endgroup$ – user21 Nov 29 '19 at 6:24
  • $\begingroup$ @TimLaska Thanks a lot. It really helps. $\endgroup$ – Muhamad Mohaqeq Nov 29 '19 at 13:51
  • $\begingroup$ @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). $\endgroup$ – Muhamad Mohaqeq Nov 29 '19 at 13:58
  • $\begingroup$ @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. $\endgroup$ – Tim Laska Nov 29 '19 at 21:58
  • $\begingroup$ @TimLaska Thanks. It really helps. $\endgroup$ – Muhamad Mohaqeq Dec 2 '19 at 18:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.