FEM: modelling 2D flow under periodic forcing

Question

Despite many examples of periodic boundary conditions being presented, it remains difficult to obtain a periodic solution. Let us consider a simple example of a flow in a two-dimensional channel under the action of a vertical force Cos[\[Pi] x + \[Pi]/4], periodic boundary conditions in x direction and free-slip at bottom and top

l = 2; h = 1;
\[CapitalOmega] = Rectangle[{0, 0}, {l, h}];
op = {Derivative[1, 0, 0][u][t, x, 
      y] + \[Nu] Inactive[
        Div][-Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + {u[t, x, 
        y], v[t, x, y]} . Inactive[Grad][u[t, x, y], {x, y}] + 
     Derivative[0, 1, 0][p][t, x, y],
    Derivative[1, 0, 0][v][t, x, 
      y] + \[Nu] Inactive[
        Div][-Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + {u[t, x, 
        y], v[t, x, y]} . Inactive[Grad][v[t, x, y], {x, y}] + 
     Derivative[0, 0, 1][p][t, x, y] - 
     Cos[\[Pi] x + \[Pi]/4] (1 - Exp[-t 10^-1]), 
    Derivative[0, 1, 0][u][t, x, y] + 
     Derivative[0, 0, 1][v][t, x, y]} /. \[Nu] -> 10^-1;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
bcD = {DirichletCondition[
    v[t, x, y] == 0, (y == 0 || y == h) && 0 <= x <= l], 
   DirichletCondition[p[t, x, y] == 0, x == l && y == 0]};
bcN = {NeumannValue[0, (y == 0 || y == h) && 0 <= x <= l], 0, 0};
bcP = {PeriodicBoundaryCondition[u[t, x, y], x == 0 && 0 < y < h, 
    TranslationTransform[{l, 0}]], 
   PeriodicBoundaryCondition[v[t, x, y], x == 0 && 0 < y < h, 
    TranslationTransform[{l, 0}]], 
   PeriodicBoundaryCondition[p[t, x, y], x == 0 && 0 < y <= h, 
    TranslationTransform[{l, 0}]]};
tmax = 1; Monitor[
 AbsoluteTiming[{xVel, yVel, pressure} = 
    NDSolveValue[
     Flatten@{op == bcN, bcD, bcP, ic}, {u, v, 
      p}, {x, y} \[Element] \[CapitalOmega], {t, 0, tmax}, 
     Method -> {"PDEDiscretization" -> {"MethodOfLines", 
         "SpatialDiscretization" -> {"FiniteElement", 
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t=", CForm[t]}])];], currentTime];

The result does not look periodic

Row[{StreamDensityPlot[{xVel[tmax, x, y], 
    yVel[tmax, x, y]}, {x, y} \[Element] xVel["ElementMesh"], 
   AspectRatio -> Automatic, PlotRange -> All, StreamPoints -> Fine, 
   PlotLegends -> Placed[Automatic, Top], ImageSize -> 300], 
  ContourPlot[
   pressure[tmax, x, y], {x, y} \[Element] pressure["ElementMesh"], 
   AspectRatio -> Automatic, PlotRange -> All, 
   ColorFunction -> "TemperatureMap", 
   PlotLegends -> Placed[Automatic, Top], ImageSize -> 300, 
   Contours -> 10, PlotRange -> All]}]

I've tried to put PeriodicBoundaryCondition on both sides

bcP = {PeriodicBoundaryCondition[u[t, x, y], x == l && 0 < y < h, 
    TranslationTransform[{-l, 0}]], 
   PeriodicBoundaryCondition[v[t, x, y], x == l && 0 < y < h, 
    TranslationTransform[{-l, 0}]], 
   PeriodicBoundaryCondition[p[t, x, y], x == l && 0 < y <= h, 
    TranslationTransform[{-l, 0}]], 
   PeriodicBoundaryCondition[u[t, x, y], x == 0 && 0 < y < h, 
    TranslationTransform[{l, 0}]],
   PeriodicBoundaryCondition[v[t, x, y], x == 0 && 0 < y < h, 
    TranslationTransform[{l, 0}]], 
   PeriodicBoundaryCondition[p[t, x, y], x == 0 && 0 < y <= h, 
    TranslationTransform[{l, 0}]]};

and result is getting better

but still it is not satisfying. Also the trick with ghost layer did not wowrk. I can get almost periodic flow only if I will heavily extended area taking l=10. Then right flow looks like

So, can anyone who is familiar with FEM help me set up the periodic boundary correctly?