(*Basic code from https : // www.wolfram.com/language/12/nonlinear-finite-elements/transient-navier-stokes.html*)
w = 2.2; h = 0.41; (*Sizes*)
geometry1 = RegionDifference[Rectangle[{0, 0}, {w, h}], Disk[{2/5, 1/5}, 1/20]];
BoundaryDiscretizeRegion[geometry1]

eq = {
        \[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};

tmax = 12; (*maximum time for the simulation*)
flow[t_] := 1/(1 + Exp[-1.6 (t - 5.5)]); (*how fast the input velocity rises*)

(*boundary conditions*)
inflowBC = DirichletCondition[{u[t, x, y] == flow[t]*4*1.5*y*(h - y)/h^2, v[t, x, y] == 0}, x == 0];
outflowBC = DirichletCondition[p[t, x, y] == 0., x == w];
wallBC = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, 0 < x < w];
bcs = {inflowBC, outflowBC, wallBC};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
(*Solve*)
Monitor[AbsoluteTiming[{xVel1, yVel1, pressure1} = NDSolveValue[{eq == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} \[Element] geometry1, {t, 0, tmax}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];], currentTime]

centre = 1/5; l = 1/20;
geometry2 = RegionDifference @@ (BoundaryDiscretizeRegion /@ {Rectangle[{0, 0}, {w, h}], Rectangle[{2 centre - l, centre - l}, {2 centre + l, centre + l}] })
Monitor[AbsoluteTiming[{xVel2, yVel2, pressure2} = NDSolveValue[{eq == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} \[Element] geometry2, {t, 0, tmax}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];], currentTime]

geometry3 = RegionDifference @@ (BoundaryDiscretizeRegion /@ {Rectangle[{0, 0}, {w, h}], ParametricRegion[0.065 {r Cos[t], r Sin[t] Sin[t/2]^1} + {0.415, 1/5}, {{t, 0, 2 \[Pi]}, {r, 0, 1}}]})
Monitor[AbsoluteTiming[{xVel3, yVel3, pressure3} = NDSolveValue[{eq == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} \[Element] geometry3, {t, 0, tmax}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];], currentTime]

p0 = Table[
   GraphicsColumn[{
     Show[
      DensityPlot[Norm[{xVel2[t, x, y], yVel2[t, x, y]}]/2, {x, 0, 2.2}, {y, 0, 0.41}, PlotPoints -> 50, PlotRange -> {0, 2.1}, AspectRatio -> Automatic, Frame -> None, ColorFunction -> "TemperatureMap", ColorFunctionScaling -> False]
      ,
      VectorPlot[{xVel2[t, x, y], yVel2[t, x, y]}, {x, 0.05, 2.15}, {y, 0.02, 0.4}, AspectRatio -> Automatic, Frame -> None, VectorStyle -> Black]
      ,
      Graphics[{White, Rectangle[{2 centre - l, centre - l}, {2 centre + l, centre + l}] }]
      ]
     ,
     Show[
      DensityPlot[Norm[{xVel1[t, x, y], yVel1[t, x, y]}]/2, {x, 0, 2.2}, {y, 0, 0.41}, PlotPoints -> 50, PlotRange -> {0, 2.1}, AspectRatio -> Automatic, Frame -> None, ColorFunction -> "TemperatureMap", 
ColorFunctionScaling -> False]
      ,
      VectorPlot[{xVel1[t, x, y], yVel1[t, x, y]}, {x, 0.05, 2.15}, {y, 0.02, 0.4}, AspectRatio -> Automatic, Frame -> None, VectorStyle -> Black]
      ,
      Graphics[{White, Disk[{2/5, 1/5}, 1/20], Black, Circle[{2/5, 1/5}, 1/20]}]
      ]
     ,
     Show[
      DensityPlot[Norm[{xVel3[t, x, y], yVel3[t, x, y]}]/2, {x, 0, 2.2}, {y, 0, 0.41}, PlotPoints -> 50, PlotRange -> {0, 2.1}, AspectRatio -> Automatic, Frame -> None, ColorFunction -> "TemperatureMap", 
ColorFunctionScaling -> False]
      ,
      VectorPlot[{xVel3[t, x, y], yVel3[t, x, y]}, {x, 0.05, 2.15}, {y, 0.02, 0.4}, AspectRatio -> Automatic, Frame -> None]
      ,
      ParametricPlot[0.065 {r Cos[\[Tau]], r Sin[\[Tau]] Sin[\[Tau]/2]^1} + {0.415, 1/5}, {\[Tau], 0, 2 \[Pi]}, {r, 0, 1}, Frame -> None, Background -> None, Axes -> False, PlotStyle -> {Directive[White, Opacity[1]]}, Mesh -> None, Epilog -> {White, Thick, Line[{{0.4, 1/5}, {0.479, 1/5}}]}]
      ]
     }, ImageSize -> Large]
   , {t, 3, 11, 0.1}];
ListAnimate[p0]