Tangential boundary conditions for a potential function

Question

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

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]

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]

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.

Answer 1

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]

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

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]

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]

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

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