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.
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 tryYourPFunction[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$