4
$\begingroup$

I wish to solve a potential flow problem where I impose a tangential boundary condition on part of the region. This is similar to this Wolfram blog post but can I do it using the Finite Element Method?

Here is a region

Needs["NDSolve`FEM`"];
L = 4;
r1 = ImplicitRegion[0 <= x <= L && 0 <= y <= L, {x, y}];
r2 = ImplicitRegion[x^2 + y^2 <= 1, {x, y}];
reg = RegionDifference[r1, r2];
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.1];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

Mathematica graphics

I am thinking of this as a fluid flow region and would like to put in a tangential velocity on the circular arc. The problem is that NeumannValue boundary conditions are normal to the edge so this is exactly in the wrong direction. Thus the wrong thing to do is

sol = NDSolveValue[{
   D[u[x, y], x, x] + D[u[x, y], y, y] ==
    NeumannValue[1, x^2 + y^2 == 1],
   DirichletCondition[u[x, y] == 0, x == L && y == 0]
   },
  u, {x, y} ∈ mesh
  ];

Which has a normal velocity on the circular arc. This can be seen if we work out the velocity and plot the stream function.

ClearAll[vel];
vel[x_, y_] := Evaluate[Grad[sol[x, y], {x, y}]]
StreamPlot[vel[x, y], {x, y} ∈ mesh]

Mathematica graphics

In fact this is an impossible solution because the flow is coming in on the circular arc and then leaving in the bottom right corner. How can it flow out of the corner when the normal velocity on all surfaces is zero everywhere? However, that is not the question. The problem is how to put a tangential velocity on the circular arc?

Thanks

Edit

Following useful comments from user21 the answer may be to do a viscous solution using a Stokes flow. Using a variant of his solution gives

ClearAll[u, v, p, x, y];
op = {
   Inactive[
      Div][{{-1, 0}, {0, -1}}.Inactive[Grad][u[x, y], {x, y}], {x, 
      y}] + Derivative[1, 0][p][x, y],
   Inactive[
      Div][{{-1, 0}, {0, -1}}.Inactive[Grad][v[x, y], {x, y}], {x, 
      y}] + Derivative[0, 1][p][x, y],
   Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]
   };
pde = op == {0, 0, 0};
bcs = {
   DirichletCondition[u[x, y] == -y/Sqrt[2], x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[v[x, y] == x/Sqrt[2], x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[u[x, y] == 0., x == 0 || x == L],
   DirichletCondition[v[x, y] == 0., y == 0 || y == L],
   DirichletCondition[p[x, y] == 0, x == 1 && y == 0]};
{xVel, yVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} ∈ 
    mesh,
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];
StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} ∈ mesh]

Mathematica graphics

What I have done is to put a tangential velocity along the circular arc and enforced zero normal velocities along the other surfaces but allowed tangential velocities. This is along the correct lines but is a viscous flow solution rather than a potential flow solution. Thanks again to user21.

$\endgroup$
0

1 Answer 1

2
$\begingroup$

I am not sure I completely understand what you would like to do but does this help:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
    DirichletCondition[u[x, y] == 1, x == 0], 
    DirichletCondition[u[x, y] == 0, x == L]}, 
   u, {x, y} \[Element] mesh];
ClearAll[vel];
vel[x_, y_] := Evaluate[-Grad[sol[x, y], {x, y}]]
StreamPlot[vel[x, y], {x, y} \[Element] mesh]

enter image description here

What I have done is to use an inflow and outflow Dirichlet boundary conditions (i.e. set a potential) and then I computed the negative gradient. This model implies a natural Neumann zero boundary condition along the quarter cylinder. Which in the gradient is of the potential shows tangential flow.

Here is another visualization:

Show[
 ContourPlot[vel[x, y], {x, y} \[Element] mesh, 
  ColorFunction -> "TemperatureMap"],
 StreamPlot[vel[x, y], {x, y} \[Element] mesh]
 ]

enter image description here

Update

Here is a version that goes round and round:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
    DirichletCondition[u[x, y] == 1, x == 0], 
    DirichletCondition[u[x, y] == 0, y == 0]}, 
   u, {x, y} \[Element] mesh];
ClearAll[vel];
vel[x_, y_] := Evaluate[-Grad[sol[x, y], {x, y}]]
StreamPlot[vel[x, y], {x, y} \[Element] mesh]

enter image description here

Update 2

Another idea is to use a Stokes flow:

op = {
   Inactive[
      Div][{{-1, 0}, {0, -1}}.Inactive[Grad][u[x, y], {x, y}], {x, 
      y}] + Derivative[1, 0][p][x, y], 
   Inactive[
      Div][{{-1, 0}, {0, -1}}.Inactive[Grad][v[x, y], {x, y}], {x, 
      y}] + Derivative[0, 1][p][x, y], 
   Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]};
pde = op == {0, 0, 0};
bcs = {
   DirichletCondition[u[x, y] == Cos[x*\[Pi]/2], 
    x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[v[x, y] == Sin[x*\[Pi]/2], 
    x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == L || y == L]
   , DirichletCondition[p[x, y] == 0, x == 1 && y == 0]
   };
{xVel, yVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} \[Element] 
    mesh, Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];
StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] mesh]

enter image description here

Update 3

The above does not have inflow or outflow conditions. At x==0 and y==0 there are free flow conditions. The rotating cylinder enforces that fluid motion.

You can also see this if the full geometry is used like so:

L = 4;
r1 = Rectangle[{-L, -L}, {L, L}];
r2 = Disk[];
reg = RegionDifference[r1, r2];
bmesh = ToBoundaryMesh[reg, "MaxBoundaryCellMeasure" -> 0.1];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

enter image description here

And the solve and visualize the system with the following boundary conditions. At the outer walls we have zero velocity in u and v direction. At inner cylinder we have a driving force and there is a reference pressure at x==1 && y==0.

bcs = {
   DirichletCondition[u[x, y] == If[y >= 0, 1, -1]*Cos[x*\[Pi]/2], 
    x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[v[x, y] == If[x >= 0, -1, 1]*Cos[y*\[Pi]/2], 
    x^2 + y^2 - 1 <= 10^-3],
   DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
    x == L || y == L || x == -L || y == -L]
   , DirichletCondition[p[x, y] == 0, x == 1 && y == 0]
   };
{xVel, yVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} \[Element] 
    mesh, Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
     "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];

Show[
 ToBoundaryMesh[mesh]["Wireframe"],
 StreamPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] mesh]]

enter image description here

$\endgroup$
10
  • $\begingroup$ Thank you for looking at this. You have managed to get a tangential flow. However, can we do this without having an inflow and outflow? I would like the flow just to go round-and-round within the region. I have been looking at this too and will post a bit more later today. $\endgroup$
    – Hugh
    Mar 12, 2018 at 8:48
  • $\begingroup$ @Hugh, you'd need something to drive the flow. $\endgroup$
    – user21
    Mar 12, 2018 at 9:15
  • $\begingroup$ A circulation drives the flow. The circulation can be driven by sliding surface. However, you are correct to notice that the circulation for Laplace's equation is always zero. This is what I have been thinking on. $\endgroup$
    – Hugh
    Mar 12, 2018 at 9:44
  • $\begingroup$ @Hugh, you could also use a Stokes flow. See update. $\endgroup$
    – user21
    Mar 12, 2018 at 10:12
  • 1
    $\begingroup$ This Stokes flow is wonderful. Being ambitions I would like to do a transient problem. I can easily add the du/dt terms. Do we have an example somewhere? It seems different to the swinging beam problem. Perhaps start with a plate oscillating in its own plane (i.e. Stokes second problem). $\endgroup$
    – Hugh
    Mar 12, 2018 at 20:28

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.