4
$\begingroup$

I know that there are lots of questions on this forum regarding mathematica FEM, but I haven't quite found what I am looking for. I will be working with data obtained through PIV (particle image velocimetry) which will give me the velocity field for a given 2D (and eventually 3D) region everywhere at t=0. Assuming for now this region is rectangular, I also know the velocity boundary conditions for all time at the edges of the region. I am looking to make a 2D flow solver that can propagate this initial velocity condition for a certain time T. The governing equations are the incompressible Navier/Stokes. Most of the solutions on this forum start by setting the whole field (pressure/velocity) to zero and allowing the boundary conditions to induce a flow, whereas I have non-zero initial conditions.

Effectively the question is as follows: given the initial velocity distribution in a rectangular region, and knowing the boundary velocity conditions at all time, how could one propagate this velocity field in time?

In 1D this is very simple using NDSolve as follows (using the viscous Burger's equation, not asking NDSolve to use FEM yet). This starts with an initial sinusoidal velocity field which decays due to viscosity over time as expected. I know that I am giving NDSolve a continuous initial condition, but this could be achieved from discrete data using interpolation. 1D basic example:

Needs["NDSolve`FEM`"]
samples = 60;
nu = 0.1;
tmax = 0.2;
xmax = 1;
(*define and visualise the initial velocity condition*)
u0[x_] := Sin[3Pi*x/xmax];
Plot[u0[x],{x,0,xmax}]


(*1D propagation, solve the governing equation for the known starting conditions*)
(*for now set the boundary conditions at x=0 and x=xmax to 0; to set non zero boundary conditions you can remove the "0*" before the Sin[100t]*)
bcdirichlet = {DirichletCondition[u[x, t] == u0[x],t==0],
               DirichletCondition[u[x,t]==0,x==0],
                DirichletCondition[u[x,t]==0*Sin[100t],x==xmax]};

s = NDSolve[{D[u[x,t],t]+u[x,t]*D[u[x,t],x]- nu * D[u[x,t],{x,2}]==0,bcdirichlet}, u[x,t],{x,0,xmax},{t,0,tmax}];

(*plot the solution of the propagation *)
actualsol = Plot3D[Evaluate[u[x,t]/.%],{x,0,xmax},{t,0,tmax},PlotRange->All,MaxRecursion->5,PlotTheme->"Scientific",AxesLabel->{Framed[Style["x (m)",20,Italic],FrameStyle->None],Framed[Style["t (s)",20,Italic],FrameStyle->None],Framed[Style["u (m/s)",20,Italic],FrameStyle->None]},TicksStyle->Directive[15]]

In 2D we have to worry about pressure as well. I know there is the SIMPLE algorithm (https://en.wikipedia.org/wiki/SIMPLE_algorithm), and I am looking into that but I wonder if mathematica has ways of approaching this.

Thanks!

** EDIT **

So here is an example of what I would like to achieve, adapted from "https://www.wolfram.com/language/12/nonlinear-finite-elements/transient-navier-stokes.html?product=language":

(*define a square working region*)
rules = {length -> 100/100, height -> 100/100};
\[CapitalOmega] = 
  Rectangle[{-length, -height}, {length, height}] /. rules;

(*visualise the region*)
Show[
 BoundaryDiscretizeRegion[\[CapitalOmega]],
 bmesh = HighlightMesh[
   BoundaryDiscretizeRegion[\[CapitalOmega], 
    AccuracyGoal -> 5], {Style[1, Black], Style[2, None]}],
 ImageSize -> Medium]

(*define the initial velocities, ensuring the the divergence is zero*)

velocityU[x_,y_]:=Cos[x*Pi/2]Cos[y*Pi/2]
velocityV[x_,y_]:=Sin[x*Pi/2]Sin[y*Pi/2]
(*this P function is zero at the moment as it is just the divergence*)
Pfunction[x_,y_]:=\!\(\*SuperscriptBox[\(velocityU\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x,y]+\!\(\*SuperscriptBox[\(velocityV\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x,y];

n=1
(*visualise the flow field*)
VectorPlot[{velocityU[x,y],velocityV[x,y]},{x,-1n,1n},{y,-1n,1n},VectorPoints->"Regular"]
Plot3D[Pfunction[x,y],{x,-1n,1n},{y,-1n,1n}]


(*navier stokes*)
op = {
    \[Rho] 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, 
       y] + \[Rho] {u[t, x, y], v[t, x, y]}.Inactive[Grad][
        u[t, x, y], {x, y}] + 
     Inactive[Div][(-\[Mu] Inactive[Grad][u[t, x, y], {x, y}]), {x, 
       y}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y], \[Rho] 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, 
       y] + \[Rho] {u[t, x, y], v[t, x, y]}.Inactive[Grad][
        v[t, x, y], {x, y}] + 
     Inactive[Div][(-\[Mu] Inactive[Grad][v[t, x, y], {x, y}]), {x, 
       y}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y], 
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y]} /. {\[Mu] -> 10^-3, \[Rho] -> 1};

(*define some boundary conditions, is this correct?, are there enough?*)
topBC = 
DirichletCondition[{u[t,x,y] == velocityU[x,y],v[t,x,y] == velocityV[x,y]},-length<x<length];

(*I thought to set the pressure at the sides (x=1,-1) to zero as the velocity has no curvature there*)

sideBC = 
DirichletCondition[p[t,x,y] ==0 ,{x==length, x==-length}]

bcs = {topBC,sideBC} /. rules;
ic = {u[0, x, y] == velocityU[x,y], v[0, x, y] == velocityV[x,y], p[0, x, y] == Pfunction[x,y]};


Monitor[AbsoluteTiming[{xVel, yVel, pressure} = 
    NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v, 
      p}, {x, y} \[Element] \[CapitalOmega], {t, 0, 6},
     Method -> {"PDEDiscretization" -> {"MethodOfLines",
         "SpatialDiscretization" -> {"FiniteElement", 
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
           "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}])];], currentTime]

In this current format, I know I haven't matched the pressure field to the initial conditions properly, but I just want to see if this format is correct. With this setup, we immediately get errors (I'm not sure if this is over or underdefined as well). Any thoughts? I'm getting zero pivot errors and Set::shape: Lists {xVel,yVel,pressure} and NDSolveValue[<<1>>] are not the same shape. The latter error is a little unexpected.

Last point: I know my temporal boundary conditions are currently fixed in time, so I'm not expecting anything to move, but I'd like to be able to see it calculate the stationary flow anyway.

Improve this question
$\endgroup$
12
  • $\begingroup$ How important are turbulence and boundary layers in your set-up? What Reynolds number are you simulating? $\endgroup$
    – Hugh
    Feb 4, 2021 at 11:30
  • 1
    $\begingroup$ Is this question about SIMPLE algorithm implementation in 2D or about initial condition for NSE? $\endgroup$
    – Alex Trounev
    Feb 4, 2021 at 11:48
  • $\begingroup$ @AlexTrounev someone who is working on the project with me has an implementation of the SIMPLE algorithm in 2D in matlab, but I am wondering if mathematica can emulate this algorithm for my specific problem or use its own methods. The main question is to propagate turbulent (incompressible) flow given an intial velocity field. $\endgroup$
    – Adrian del Ser
    Feb 4, 2021 at 12:30
  • 1
    $\begingroup$ I do not think there is such an example in the documentation. When you use the code use something like this as an initial condition ics={u[0,x,y]==YourUFunction[x,y], v[0,x,y]==YourVFunction[x,y], p[0,x,y]==YourPFunction[x,y]} If you do not have a YourPFunction then you could try YourPFunction[x,y]:=D[YourUFunction[0,x,y],x]+D[YourVFunction[0,x,y],y] and see where that leads. But note there must be some driving force to drive the fluid. So if you set your BCs also to with YourU/V/P/Function then not much can happen as the system then is in balance. Post your code here with all data. $\endgroup$
    – user21
    Feb 5, 2021 at 9:33
  • 1
    $\begingroup$ @user21 I've added some (currently not working code with zero divergence initial velocity field. I'm not entirely sure how to match the boundaries yet or what would make it overdefined. Surely the equation for YourPFunction you suggest is just 0, since it represents the divergence of the field? I'll try use the pressure poisson equation to get a better pressure field but it would be good to see the code at least work (even if the pressure field is bad) $\endgroup$
    – Adrian del Ser
    Feb 7, 2021 at 9:12

2 Answers 2

Reset to default
5
$\begingroup$

For 2D time dependent incompressible flow we can use linear FEM, as it discussed on Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM. This solver has been tested and compared with other different approaches. The linear FEM solver shows also nice solution for the problem for turbulent like initial data. We use as initial data the random, divergence free vector field given by
(define a square working region)

 rules = {length -> 100/100, height -> 100/100};
    \[CapitalOmega] = 
        Rectangle[{-length, -height}, {length, height}] /. rules;
    
    
    (*define the initial velocities, ensuring the the divergence is zero*)
    
    n=1;
    SeedRandom[1234]; x0 = RandomReal[{-1, 1}, 10]; y0 = 
     RandomReal[{-1, 1}, 10]; amp = RandomReal[{-1/2, 1/2}, 10]; k = 
     RandomReal[{0, 2}, 10];
    velocityU[x_, y_] := 
      Sum[amp[[i]] Cos[k[[i]] x*Pi/2 + x0[[i]]] Cos[
         k[[i]] y*Pi/2 + y0[[i]]], {i, Length[x0]}];
    velocityV[x_, y_] := 
      Sum[amp[[i]] Sin[k[[i]] x*Pi/2 + x0[[i]]] Sin[
         k[[i]] y*Pi/2 + y0[[i]]], {i, Length[x0]}];
    StreamDensityPlot[{velocityU[x, y], velocityV[x, y]}, {x, -1 n, 
      1 n}, {y, -1 n, 1 n}, ColorFunction -> "Rainbow", 
     PlotLegends -> Automatic]

Figure 1

We run solver on 10 steps with time step of t0=1/100 on v.12.2

UX[0][x_, y_] := velocityU[x, y]; 
VY[0][x_, y_] := velocityV[x, y]; 
P0[0][x_, y_] := 0; nn = 10; t0 = 1/100; 
Do[{UX[i], VY[i], P0[i]} = NDSolveValue[{{Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
          Derivative[1, 0][p][x, y] + (u[x, y] - UX[i - 1][x, y])/t0 + UX[i - 1][x, y]*D[u[x, y], x] + VY[i - 1][x, y]*D[u[x, y], y], 
         Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y] + 
          (v[x, y] - VY[i - 1][x, y])/t0 + UX[i - 1][x, y]*D[v[x, y], x] + VY[i - 1][x, y]*D[v[x, y], y], 
         Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]} == {0, 0, 0} /. \[Mu] -> 10^(-3), 
      {DirichletCondition[{u[x, y] == velocityU[x, y], v[x, y] == velocityV[x, y]}, y == length || y == -length], 
        DirichletCondition[p[x, y] == P0[i - 1][x, y], x == length || x == -length]} /. rules}, {u, v, p}, Element[{x, y}, \[CapitalOmega]], 
     Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], 
   {i, 1, nn}];

Finally we can define mesh and plot solution on every step

mesh = UX[1]["ElementMesh"]

(*Out[]= NDSolve`FEM`ElementMesh[{{-1., 1.}, {-1., 
   1.}}, {NDSolve`FEM`QuadElement["<" 4096 ">"]}]*)
Table[{Show[
   DensityPlot[
    Norm[{UX[i][x, y], VY[i][x, y]}], {x, -1 n, 1 n}, {y, -1 n, 1 n}, 
    ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
    PlotLabel -> i], 
   VectorPlot[{UX[i][x, y], VY[i][x, y]}, Element[{x, y}, mesh]]],
  Plot3D[P0[i][x, y], {x, -1 n, 1 n}, {y, -1 n, 1 n}, Mesh -> None, 
   ColorFunction -> "Rainbow"]}, {i, 1, nn}]

Picture looks the same one on every step since boundary conditions not depend on time, but we restore pressure from zero at t=0 to some function dependent on velocity distribution Figure 2

Improve this answer
$\endgroup$
7
  • $\begingroup$ I am not sure this answers the question. How do you propose to use your result to advance the initialized solver in time? $\endgroup$
    – user21
    Feb 8, 2021 at 6:16
  • $\begingroup$ I'll be looking at this a little bit later today or tomorrow, but I assume that if you force a certain number of the boundary conditions to vary in time in a way that doesn't violate the conservation equations, the flow should propagate correspondingly. The question is, how many of the boundary conditions are indepenent of each other? I assume we just have to be careful to conserve mass, ie, the net flow into the control volume must equal the flow out when looking at incompressible flow. I'll investigate. $\endgroup$
    – Adrian del Ser
    Feb 8, 2021 at 10:07
  • $\begingroup$ @user21 This answers the question in this particular case with time independent boundary conditions. My example even more complex then that in the question states. But I don't understand why nonlinear FEM can't run on this example. $\endgroup$
    – Alex Trounev
    Feb 8, 2021 at 10:35
  • $\begingroup$ @AlexTrounev indeed. I wanted to see if the flow solver would work with a correct initial setup. Of course, the boundary conditions will be time dependent as the system is open but I believe that the system should be able to handle that anyway (I assume we would need to define how 3 of the 4 boundaries change in time, otherwise it is overdefined?). When referring to turbulent-like intial conditions, are you referring to the initial velocity being expressed as an average field, with zero mean disturbances added on top? $\endgroup$
    – Adrian del Ser
    Feb 8, 2021 at 10:44
  • $\begingroup$ @AdriandelSer Boundary conditions for small pies of turbulent flow are unknown. In some applications like CUDA, DNS, LES, mathematical theorems and so on there are periodic bc have used. $\endgroup$
    – Alex Trounev
    Feb 8, 2021 at 10:50
3
$\begingroup$

I'm posting this as an answer only because it is good enough for my requirements at the moment. I have taken the code provided by Alex Trounev in the first answer and modified it to decay the two dirichlet conditions in time. As the top and bottom velocities decay, the flow starts resembling 2D channel flow (as the top and bottom effectively approach the no slip condition, and represent walls). I have not imposed any additional driving forces on the flow.

As before, we define the working region:

    \[CapitalOmega] = 
        Rectangle[{-length, -height}, {length, height}] /. rules;
    
(*define the initial velocities, ensuring the the divergence is zero*)
    
n=1;
SeedRandom[1234];
x0 = RandomReal[{-1, 1}, 10];
y0 = RandomReal[{-1, 1}, 10];
amp = RandomReal[{-1/2, 1/2}, 10];
k = RandomReal[{0, 2}, 10];
velocityU[x_, y_] := 
      Sum[amp[[i]] Cos[k[[i]] x*Pi/2 + x0[[i]]] Cos[
         k[[i]] y*Pi/2 + y0[[i]]], {i, Length[x0]}];
         
velocityV[x_, y_] := 
      Sum[amp[[i]] Sin[k[[i]] x*Pi/2 + x0[[i]]] Sin[
         k[[i]] y*Pi/2 + y0[[i]]], {i, Length[x0]}];
         
StreamDensityPlot[{velocityU[x, y], velocityV[x, y]}, {x, -1 n, 1 n}, {y, -1 n, 1 n}, ColorFunction -> "Rainbow", PlotLegends -> Automatic]

stream plot of initial velocity field

Here I have modified nn, t0, mu and the boundary conditions to decay in time:

UX[0][x_, y_] := velocityU[x, y]; 
VY[0][x_, y_] := velocityV[x, y]; 
P0[0][x_, y_] := 0; nn = 90; t0 = 5/100; 
Do[{UX[i], VY[i], P0[i]} = NDSolveValue[{{Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
          Derivative[1, 0][p][x, y] + (u[x, y] - UX[i - 1][x, y])/t0 + UX[i - 1][x, y]*D[u[x, y], x] + VY[i - 1][x, y]*D[u[x, y], y], 
         Inactive[Div][{{-\[Mu], 0}, {0, -\[Mu]}} . Inactive[Grad][v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y] + 
          (v[x, y] - VY[i - 1][x, y])/t0 + UX[i - 1][x, y]*D[v[x, y], x] + VY[i - 1][x, y]*D[v[x, y], y], 
         Derivative[1, 0][u][x, y] + Derivative[0, 1][v][x, y]} == {0, 0, 0} /. \[Mu] -> 10^(-1), 
      {DirichletCondition[{u[x, y] == (1/(1+10(i-1)*t0))*velocityU[x, y], v[x, y] == (1/(1+10(i-1)*t0))*velocityV[x, y]}, y == length || y == -length], 
        DirichletCondition[p[x, y] == P0[i - 1][x, y], x == length || x == -length]} /. rules}, {u, v, p}, Element[{x, y}, \[CapitalOmega]], 
     Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}], 
   {i, 1, nn}];

Finally plot every 4 frames (to make sure the exported gif is below 2mb):

mesh = UX[1]["ElementMesh"]

twodflow = Table[{Show[
   DensityPlot[
    Norm[{UX[i][x, y], VY[i][x, y]}], {x, -1 n, 1 n}, {y, -1 n, 1 n}, 
    ColorFunction -> "Rainbow", (*ColorFunctionScaling\[Rule]False,*) PlotLegends -> Automatic, 
    PlotLabel -> i], 
   VectorPlot[{UX[i][x, y], VY[i][x, y]}, Element[{x, y}, mesh]]],
  Plot3D[P0[i][x, y], {x, -1 n, 1 n}, {y, -1 n, 1 n}, Mesh -> None, 
   ColorFunction -> "Rainbow",PlotRange->All]}, {i, 0, nn,4}]
   
Export["2dflow.gif",twodflow,"DisplayDurations"->0.2]

EDIT I'm not sure if the gif below is looping appropriately (you should be able to click on it to at least see it run once through), if anyone knows why it may not be looping automatically let me know.

gif of transient 2D flow

Improve this answer
$\endgroup$
1
  • $\begingroup$ It is very nice example (+1). We can also include force in this model. It woks up to $Re=10^4$ for some problems. $\endgroup$
    – Alex Trounev
    Feb 8, 2021 at 21: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.