9
$\begingroup$

I am trying to describe the fluid flow within a radially symmetric tube and use the code as it can be found from Wolframs homepage as template (Stokes's equation).

I adapted the geometry and the inlet flow to my given system and get a reasonable solution which is nice (parabolic laminar velocity profile):

laminar flow

Now, with the goal to reduce computational demand and to use later on the common balance equations for radially symmetric systems I want to solve only the upper half of the system depicted above. By doing so the centre line of the tube, which is at r = 0.0024, is going to be the new x axis.

In principle this can be done by using the symmetry of the system and setting up Neumann boundary conditions at the bottom (r == 0):

$\frac{\partial u}{\partial r}\left ( r=0 \right )=\frac{\partial v}{\partial r}\left ( r=0 \right )=\frac{\partial p}{\partial r}\left ( r=0 \right )=0$

Unfortunately here I am at the point where I get stuck. By using the FE method NeumannValues need to be defined. According to the documentation NeumannValue specifies the value on the right side of this equation with u being the dependent variable:

$\vec{n}\left ( c\triangledown u+\alpha u-\gamma \right )=g-qu$

But how does it look like when I have PDEs with more than one dependent variable like in $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial r}=0$ ?

Although I am not sure about the NeumannValue definition for a first test I set it to zero using this code:

Needs["NDSolve`FEM`"]

Ω = Rectangle[{0, 0}, {2.0016, 0.0024}]; 
RegionPlot[Ω, AspectRatio -> 0.5]

op = {D[μ*D[u[x, r], x], x] + D[μ*D[u[x, r], r], r] - 
 D[p[x, r], x], 
D[μ*D[v[x, r], x], x] + D[μ*D[v[x, r], r], r] - 
 D[p[x, r], r], D[u[x, r], x] + D[v[x, r], r]} /. μ -> 0.0001; 

Subscript[Γ, N] = {NeumannValue[0, r == 0 && x > 0], 
NeumannValue[0, r == 0 && x > 0], 
NeumannValue[0, r == 0 && x > 0]}; 

Subscript[Γ, 
D] = {DirichletCondition[u[x, r] == 2*0.065*(1 - (r/0.0024)^2), 
x == 0.], 
DirichletCondition[{u[x, 0.0024] == 0., v[x, 0.0024] == 0.}, 
Inequality[0, Less, x, LessEqual, 2.0016]], 
   DirichletCondition[p[x, r] == 0., x == 2.0016]}; 

mesh = ToElementMesh[Ω, 
MaxCellMeasure -> {"Length" -> 0.0012}]; 

{xVel, yVel, pressure} = 
NDSolveValue[{op == Subscript[Γ, N], 
Subscript[Γ, D]}, {u, v, p}, Element[{x, r}, mesh], 
Method -> {"FiniteElement", 
 "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}]; 

ContourPlot[xVel[x, r], {x, 0, 2}, {r, 0, 0.0024}, 
AspectRatio -> 0.5, ColorFunction -> "TemperatureMap", 
PlotLegends -> Automatic]

But apparently this is wrong as there is a flux at the bottom border causing the axial velocity to be zero where it should be at its maximum:

NeumannValue zero

Any help is appreciated, thanks.

$\endgroup$
  • $\begingroup$ If you want to reduce the problem, can't you solve it over an axi-symetric region? wolfram.com/mathematica/new-in-10/pdes-and-finite-elements/… $\endgroup$ – Feyre Aug 3 '16 at 10:33
  • $\begingroup$ I think that is in principle what I did. In your example they use a Neumann boundary condition as well. Only in my case a NeumannValue of zero is not really preventing a flux of the single variables u, v and p which might be caused by the structure of the PDEs. Otherwise the result would look like the upper half of the first picture. For the first picture I used two Dirichlet conditions as suggested by the template from Wolfram. $\endgroup$ – Sebbo Aug 3 '16 at 11:29
  • $\begingroup$ Is your PDE condition $\frac{\partial u}{\partial x} + \frac{\partial v}{\partial r} = 0$ not equivalent to the condition $\mathbf{n} \cdot \nabla u = 0$.? Which is then NeumannValue[0, x < 0] ? $\endgroup$ – M Horsley Aug 3 '16 at 12:52
8
$\begingroup$

No need for Neumann BCs at all. The independent variables got messed up a bit. This works:

Needs["NDSolve`FEM`"]

xL = 2.0016; yL = 0.0024;
Ω = Rectangle[{0, 0}, {xL, yL}];
RegionPlot[Ω, AspectRatio -> 0.5]

op = {D[μ*D[u[x, r], x], x] + D[μ*D[u[x, r], r], r] - 
     D[p[x, r], x], 
    D[μ*D[v[x, r], x], x] + D[μ*D[v[x, r], r], r] - 
     D[p[x, r], r], D[u[x, r], x] + D[v[x, r], r]} /. μ -> 
    0.0001;
Subscript[Γ, D] = {
   DirichletCondition[u[x, r] == 2*0.065*(1 - (r/yL)^2), x == 0.], 
   DirichletCondition[{u[x, r] == 0., v[x, r] == 0.}, 
    0 < x <= xL && r == yL], 
   DirichletCondition[p[x, r] == 0., x == xL]};

mesh = ToElementMesh[Ω, 
   MaxCellMeasure -> {"Length" -> 0.0012}];

{xVel, yVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, Subscript[Γ, D]}, {u, 
    v, p}, Element[{x, r}, mesh], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

ContourPlot[xVel[x, r], {x, 0, xL/100}, {r, 0, yL}, 
 AspectRatio -> 0.5, ColorFunction -> "TemperatureMap", 
 PlotLegends -> Automatic]

enter image description here

$\endgroup$
  • 1
    $\begingroup$ This is great!! Thank you so much! :-) So, what I actually messed up was the Dirichlet boundary condition where I set the bottom to zero as well. So am I right, that in this case a no flux Neumann condition is applied as I originally wanted? $\endgroup$ – Sebbo Aug 3 '16 at 13:00
  • $\begingroup$ You are welcome. Do you have a specific application you work on? $\endgroup$ – user21 Aug 3 '16 at 13:02
  • $\begingroup$ I want to couple the flow field with a reaction and diffusion process. Through the nature of the reaction the viscosity in the system is not homogenous distributed which leads to interesting flow patterns that can cause problems like fouling. $\endgroup$ – Sebbo Aug 3 '16 at 13:09
  • $\begingroup$ Very interesting. Do you think it's possible to see the result once you are done with it? $\endgroup$ – user21 Aug 3 '16 at 13:24
  • $\begingroup$ That should be possible. Until then I may come up with more questions. :-) $\endgroup$ – Sebbo Aug 3 '16 at 13:59

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.