# FEM: modelling 2D flow under periodic forcing

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?

• How you know what period is for this flow? Periodic in space force produces some flow that may be not periodic or periodic with some unknown period. On the other hand periodic in time force generates periodic flow like in my post on mathematica.stackexchange.com/questions/246091/… Commented May 26, 2022 at 10:32
• You're right. I forgot to mention that I consider the Reynolds number to be very low, so the flow remains quasi-linear and periodic. Commented May 26, 2022 at 10:51
• The flow has a period of force or shorter but multiples of it.. Commented May 26, 2022 at 10:58
• Actually for l=10 the Reynolds number is about 2, and it is why on the plots Plot[xVel[tmax, x, .1], {x, 0, l}, PlotRange -> All] , Plot[pressure[tmax, x, .5], {x, 0, l}, PlotRange -> All] we can recognize anharmonic oscillations. What kind of solution do you try to reproduce with your bcD? Commented May 26, 2022 at 13:37
• I look for bcP definition which provide periodic solution so that fields (u,v,f) satisfy f(x+l,y)=f(x,y) and f'(x+l,y)=f'(x,y). If you want I can you analytical solution. Commented May 26, 2022 at 13:47