3
$\begingroup$

Background: I'm looking to have a 2D re-entrant channel of a geophysical flow that is forced by a wind stress at the surface ("taux"; below) and experiences planetary rotation ("cor"; below). I'm using linear, 2D Navier Stokes with bottom stress (basically Stommel 1948; 1) and whilst I've got the closed box working (not shown; [2]), I'm having difficulty with having a re-entrant boundary condition.

Problem. Running the below code I get a few errors, but the one that concerns me is regarding the PeriodicBoundaryCondition, i.e., when I call NDSolveValue I get: DynamicNDSolveValue::fembpf: The boundary condition PeriodicBoundaryCondition[v[0,y]==v[1,y]] is not valid. And again the same error for BC1.

Here is the code, the errors arise after the 'NDSolveValue' line:

 SetOptions[EvaluationNotebook[],CellContext->Notebook]
Ω=Region[Rectangle[{0,0},{1,1}]];
BC0=PeriodicBoundaryCondition[u[0,y]==u[1,y]];
 BC1=PeriodicBoundaryCondition[v[0,y]==v[1,y]];
 BC3 = DirichletCondition[{u[x, y] == 0, v[x, y] == 0},  y==0];
BC4 = DirichletCondition[{u[x, y] == 0, v[x, y] == 0},  y==1];
 bcs = {BC0,BC1,BC3,BC4};
 cor[y_]:=-Sin[((y *Pi/180)*90)+Pi/2]
 taux[y_]:=-Cos[(Pi/180*y)*270]
(*x-momentum
y-momentum
incompressible*)

navierstokes = {Derivative[1, 0][p][x, y] - cor[y]*v[x, y] - taux[y] + r*u[x, y], 
    Derivative[0, 1][p][x, y] + cor[y]*u[x, y] + r*v[x, y], Derivative[1, 0][u][x, y] + 
     Derivative[0, 1][v][x, y]}; 

 op = navierstokes /. { r -> 5} ;
 {uVel, vVel,pressure} =   NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v,p}, {x, y} ∈ Ω,Method -> {"FiniteElement",      "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},"MeshOptions" -> {"MaxCellMeasure" -> 0.00025}}];
 picture=StreamPlot[{uVel[x,y],vVel[x,y]},{x,y}∈Ω,Sequence[PlotRange->All,StreamScale->{Medium,n,Automatic},Frame->None,Axes->None,AspectRatio->Automatic,PlotLegends->Automatic]];
Show[picture,Graphics[RegionBoundary[Ω]],ImageSize->Medium]

I suspect this problem is trivially solved but as I'm very new to Mathematica, I couldn't quite find an example that fixes this problem. Thus, please explain the syntax if you're able to point me in the right direction. Any general tips for fixing FEM issues much appreciated (I've read a lot of SO posts!)

1 Stommel, H. (1948). The westward intensification of wind-driven ocean currents. Transactions, American Geophysical Union, 29(2), 202. https://doi.org/10.1029/TR029i002p00202

[2] the closed box version just swaps BC0 and BC1 with DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0] and DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 1], respectively.


Update Following @AlexTrounev excellent feedback, I learnt the syntax that works for the numerical approach to my problem.

coff=0;
ruu=5;
\[CapitalOmega] = Region[Rectangle[{0, 0}, {1, 1}]];
  bcs = {PeriodicBoundaryCondition[u[x,y],x==0,Function[x,x+{1,0}]],
PeriodicBoundaryCondition[v[x,y],x==0,Function[x,x+{1,0}]],
PeriodicBoundaryCondition[p[x,y],x==0,Function[x,x+{1,0}]],
 DirichletCondition[v[x, y] == 0,  y==1&&0<x<1],
 DirichletCondition[ v[x, y] == 0, y==0&&0<x<1]};
  cor[y_] = -Sin[((y *Pi/180)*70) + Pi/2];
  taux[y_] = -Cos[(Pi/180*y)*270] + coff;
  navierstokes = {(p^(1,0))[x, y] - cor[y]*v[x, y] - taux[y] + ruu*u[x, y],
       (p^(0,1))[x, y] + cor[y]*u[x, y] + ruu*v[x, y],
       (u^(1,0))[x, y] + (v^(0,1))[x, y]};
  {uVe, vVe, pressure} =    NDSolveValue[{navierstokes == {0, 0, 0}, bcs}, {u, v, p}, {x, y} \[Element] \[CapitalOmega], 
    Method -> {"FiniteElement",      
      "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.00025}}];
  puv = StreamPlot[{uVe[x, y], 
     vVe[x, y]}, {x, y} \[Element] \[CapitalOmega], 
    Sequence[PlotRange -> All, StreamScale -> {Medium, n, Automatic}, 
     Frame -> None, Axes -> None, 
     PlotLegends -> Placed[Automatic, Below]], PlotLabel -> "UV", 
    ImagePadding -> {{Automatic, Automatic}, {None, None}}]

Picture of the u and v velocities

$\endgroup$
7
  • $\begingroup$ Is p also periodic? $\endgroup$ Mar 28, 2023 at 16:41
  • $\begingroup$ Yes. Although, I haven't been setting a pressure boundary condition as I'm only interested in gradients of pressure. Is that a problem here? (I think I got unexpected solutions in the non-reentrant case when setting a p boundary condition.) $\endgroup$
    – Chris
    Mar 28, 2023 at 18:47
  • $\begingroup$ Update: this feels relevant not sure how to apply it here though. $\endgroup$
    – Chris
    Mar 28, 2023 at 18:56
  • $\begingroup$ Have you any picture related to this flow? $\endgroup$ Apr 5, 2023 at 16:23
  • $\begingroup$ @AlexTrounev, I just uploaded one above. Thanks for the extremely helpful feedback! I did try your exciting algebraic strategy but couldn't find a way to enforce no normal flow (bc1/bc2) on the northern and southern boundaries to get a solution more like pictured above. Very interested to hear if I missed something. I'm having similar issues extending the problem to 3d in a cube link. $\endgroup$
    – Chris
    Apr 6, 2023 at 14:53

1 Answer 1

3
$\begingroup$

There are several scenario to solve this problem. First, we solve

navierstokes={Derivative[1,0][p][x,y]-cor[y]*v[x,y]-taux[y]+r*u[x,y],Derivative[0,1][p][x,y]+cor[y]*u[x,y]+r*v[x,y],Derivative[1,0][u][x,y]+Derivative[0,1][v][x,y]};
sol0=Solve[Drop[navierstokes,-1]=={0,0},{u[x,y],v[x,y]}]
(*Out[]= {{u[x,y]->-((-r taux[y]+cor[y] (p^(0,1))[x,y]+r (p^(1,0))[x,y])/(r^2+cor[y]^2)),v[x,y]->-((cor[y] taux[y]+r (p^(0,1))[x,y]-cor[y] (p^(1,0))[x,y])/(r^2+cor[y]^2))}}*) 

Second, we use last equation to derive elliptic PDE for pressure as follows

cor[y_] := -Sin[((y*Pi/180)*90) + Pi/2];
taux[y_] := -Cos[(Pi/180*y)*270];
 eqn=navierstokes[[3]]/.{(v^(0,1))[x,y]->D[v[x,y]/.sol0[[1]],y], (u^(1,0))[x,y]->D[u[x,y]/.sol0[[1]],x]}//Simplify
(*Out[]= -(1/(2 (1+2 r^2+Cos[\[Pi] y])^2))(-5 \[Pi] Sin[\[Pi] y]-4 \[Pi] r^2 Sin[\[Pi] y]-4 \[Pi] Sin[2 \[Pi] y]-8 \[Pi] r^2 Sin[2 \[Pi] y]-\[Pi] Sin[3 \[Pi] y]+4 \[Pi] r Sin[\[Pi] y] (p^(0,1))[x,y]+4 r (1+2 r^2+Cos[\[Pi] y]) (p^(0,2))[x,y]+\[Pi] Sin[(\[Pi] y)/2] (p^(1,0))[x,y]-4 \[Pi] r^2 Sin[(\[Pi] y)/2] (p^(1,0))[x,y]+\[Pi] Sin[(3 \[Pi] y)/2] (p^(1,0))[x,y]+4 r (p^(2,0))[x,y]+8 r^3 (p^(2,0))[x,y]+4 r Cos[\[Pi] y] (p^(2,0))[x,y])*) 

Third, we derive boundary conditions

bc1=v[x,y]/.sol0[[1]]/.y->0
(*Out[]= -((1+r (p^(0,1))[x,0]+(p^(1,0))[x,0])/(1+r^2))*)
bc2=v[x,y]/.sol0[[1]]/.y->1
(*Out[]= -((p^(0,1))[x,1]/r)*)
bc3=u[x,y]/.sol0[[1]]/.y->0
(*Out[]= -((r-(p^(0,1))[x,0]+r (p^(1,0))[x,0])/(1+r^2))*)
bc4=u[x,y]/.sol0[[1]]/.y->1
(*Out[]= -((p^(1,0))[x,1]/r)*)

Since we have one second order PDE with periodic boundary condition on x, we need to take 2 boundary conditions from bc1,bc2,bc3,bc4 only. Also we can try to combine all of them. Note, that in a case of bc1, bc3 we have

sol1=Solve[{bc1==0,bc3==0},{(p^(1,0))[x,0], (p^(0,1))[x,0]}]
(*Out[11]= {{(p^(1,0))[x,0]->-1,(p^(0,1))[x,0]->0}}*)

Effectively we can combine it as

DirichletCondition[p[x, y] == -x, y == 0 && 0 < x < 1]

And bc2, bc4 we can combine as

DirichletCondition[p[x, y] == (1 - y)^2, y == 1 && 0 < x < 1]

Since function x is not periodic, we test last combination

sol2 = NDSolve[{eqn == 0 /. r -> 5, 
   PeriodicBoundaryCondition[p[x, y], x == 0, 
    Function[x, x + {1, 0}]], 
   DirichletCondition[p[x, y] == 0, y == 1 && 0 < x < 1]}, p, 
  Element[{x, y}, Rectangle[{0, 0}, {1, 1}]]]

This solution looks as follows

frame = Show[
  DensityPlot[
   Evaluate[
    Norm[{u[x, y], v[x, y]}] /. sol0[[1]] /. sol2[[1]] /. r -> 5], 
   Element[{x, y}, Rectangle[{0, 0}, {1, 1}]], 
   ColorFunction -> "Rainbow", PlotPoints -> 100], 
  StreamPlot[
   Evaluate[{u[x, y], v[x, y]}] /. sol0[[1]] /. sol2[[1]] /. r -> 5, 
   Element[{x, y}, Rectangle[{0, 0}, {1, 1}]], 
   StreamColorFunction -> None, StreamStyle -> LightGray, 
   StreamPoints -> 20]]; Show[
 MapAt[Translate[#, {1, 0}] &, frame, 1], frame, 
 MapAt[Translate[#, {-1, 0}] &, frame, 1], PlotRange -> All]

Figure 1

Second method is follows. We add artificial viscosity of about 10^-3, then we have

Needs["NDSolve`FEM`"];

mesh = ToElementMesh[Rectangle[{0, 0}, {1, 1}], AccuracyGoal -> 5, 
  PrecisionGoal -> 5, "MaxCellMeasure" -> 0.0005]

\[CapitalOmega] = Region[Rectangle[{0, 0}, {1, 1}]];
BC2 = PeriodicBoundaryCondition[u[x, y], x == 0, 
   Function[x, x + {1, 0}]];
BC1 = PeriodicBoundaryCondition[v[x, y], x == 0, 
   Function[x, x + {1, 0}]];
BC4 = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, (y == 0) && 
   0 < x < 1]; BC3 = 
 PeriodicBoundaryCondition[p[x, y], x == 0, 
  Function[x, x + {1, 0}]]; BC5 = 
 DirichletCondition[p[x, y] == 0, (y == 1) && 0 < x < 1];
bcs = {BC2, BC1, BC4, BC5};
cor[y_] := -Sin[((y*Pi/180)*90) + Pi/2]
taux[y_] := -Cos[(Pi/180*y)*270]
(*x-momentum y-momentum incompressible*)

navierstokes = {Derivative[1, 0][p][x, y] - cor[y]*v[x, y] - taux[y] +
     r*u[x, y] - 1/Re0 Laplacian[u[x, y], {x, y}], 
   Derivative[0, 1][p][x, y] + cor[y]*u[x, y] + r*v[x, y] - 
    1/Re0 Laplacian[u[x, y], {x, y}], 
   Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y] - 
    1/Re0 Laplacian[p[x, y], {x, y}]};

op = navierstokes /. {r -> 5, Re0 -> 10^3};
{uVel, vVel, pressure} = 
  NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, 
   Element[{x, y}, mesh], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];

Visualization

picture = 
 StreamDensityPlot[{uVel[x, y], 
   vVel[x, y]}, {x, y} \[Element] \[CapitalOmega], 
  ColorFunction -> "Rainbow", StreamColorFunction -> None, 
  StreamStyle -> LightGray, PlotLegends -> Automatic]


Show[MapAt[Translate[#, {1, 0}] &, picture[[1]], 1], picture[[1]], 
 MapAt[Translate[#, {-1, 0}] &, picture[[1]], 1], PlotRange -> All]

Figure 2

$\endgroup$
1
  • $\begingroup$ Thanks @AlexTrounev. We get the same (see my update) answer now, if I set your bcs to {PeriodicBoundaryCondition[v[x, y], x == 0, Function[x, x + {1, 0}]], PeriodicBoundaryCondition[u[x, y], x == 0, Function[x, x + {1, 0}]], PeriodicBoundaryCondition[p[x, y], x == 0, Function[x, x + {1, 0}]], DirichletCondition[v[x, y] == 0, (y == 0) && 0 < x < 1], DirichletCondition[v[x, y] == 0, (y == 1) && 0 < x < 1]}. Nice to see how the artificial viscosity is done too (note the typo). $\endgroup$
    – Chris
    Apr 17, 2023 at 11:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.