# Context

The following transient problem is the reciprocating (i.e., fully reversing) flow of a fluid $$0 over a thick heated block $$0 until the system reaches a cyclic steady-state (i.e., the system temperature oscillates around a mean when this point is reached). The given code simulates the system for a fluid velocity boundary condition $$u=3*\sin(2*\pi*1*t)$$, i.e., a frequency of $$1Hz$$ and a heat flux input of $$5000 W/m^2$$. tflow refers to the total time till which the flow oscillation occurs and t0 is the time-step. The code uses the flow solver developed by Alex Trounev and described in this page and as answer to one of my earlier question here.

## Code

Needs["NDSolveFEM"]
{f = 2;
L = 0.040, d = 0.003, e = 0.005, kf = 0.026499, ks = 16,
rho = 1.1492, rhos = 7860, mu = 18.923*10^-6, cp = 1.0069*10^3,
cps = 502.4}; u0 = 3; nu = mu/rho; om = 2 Pi f;
tflow = 1000;
t0 = .1;
NV = 2 f tflow;
nn = Round[NV \[Pi]/(om t0)]
Ti = 307; q = 5000/Ti;

reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 =
ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];

mesh = ToElementMesh[FullRegion[2], {{0, L}, {0, d}}];
mesh1 = ToElementMesh[FullRegion[2], {{0, L}, {-e, d}}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 307/Ti; appro =
With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cps rhos + (cp rho - cps rhos) appro[y]);
eqs = {Inactive[
u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] +
UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y] + (u[x, y] - UX[i - 1][x, y])/t0,
Inactive[
v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] +
UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], y] + (v[x, y] - VY[i - 1][x, y])/t0,
D[u[x, y], x] + D[v[x, y], y]};
bc[i_] := {DirichletCondition[{u[x, y] == u0*Sin[om*i*t0],
v[x, y] == 0}, x == L (1 - Sign[Sin[om*i*t0]])/2 && 0 < y < d],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 0 || y == d],
DirichletCondition[p[x, y] == 0,
x == L (1 + Sign[Sin[om*i*t0]])/2 && 0 < y < d]};

Monitor[Do[{UX[i], VY[i], P[i]} =
NDSolveValue[{eqs == {0, 0, 0} /. \[Mu] -> nu, bc[i]}, {u, v,
p}, {x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];, {i, 1,
nn}], ProgressIndicator[i, {1, nn}]]; // AbsoluteTiming

Monitor[Do[ux = If[y <= 0, 0, UX[i][x, y]];
vy = If[y <= 0, 0, VY[i][x, y]];
Tfs[i] =
NDSolveValue[{rde[
y] ((ux*D[T[x, y], x] +
vy*D[T[x, y], y]) + (T[x, y] - Tfs[i - 1][x, y])/t0) -
Inactive[Div][
NeumannValue[q, y == -e],
DirichletCondition[{T[x, y] == 1},
x == L (1 - Sign[Sin[om*i*t0]])/2 && 0 <= y <= d]},
T, {x, y} \[Element] mesh1,
Method -> {"FiniteElement",
"InterpolationOrder" -> {T -> 2}}] // Quiet;, {i, 1, nn}],
ProgressIndicator[i, {1, nn}]]; // AbsoluteTiming

ListLinePlot[
Table[{i t0, Tfs[i][0.5*L, -e/2]*Ti - 273.16}, {i, 0, nn}],
AxesLabel -> {"t(s)", "T(in deg. C)"}, PlotRange -> Full]

On evaluation, the temperature history of a point (0.5L,-0.5e) in the solid looks like:

## Question

In essence, the system temperature reaches a cyclic steady state (css) eventually, where it oscillates around a mean value (marked with purple color).

In modelling the system with the approach mentioned here, one has to march through the entire transient phase physically and wait till a css is reached (leading to large simulation times, as the time-step is restricted by the frequency of flow oscillation). My question is, if an alternative approach exists, through which a css can be directly calculated without iterating through the entire flow time? This will be especially useful for high f (which needs even smaller time-steps) and high heat input q (which push reaching css longer) Please note that I am not talking about a steady-state solution, as the flow direction reverses here after each half-cycle.

## Some resource

An approach to such direct calculation has been hinted in the Chapter 6-Periodic convection (pp. 259) of the book by Arpaci and Larsen called "Convection heat transfer" and Chapter 6 of the book "Conduction heat transfer" by the same authors. It is named as the method of Complex temperature. A copy of the former can be found here and for the latter here. The method assumes that the solution will oscillate at frequency equal to the imposed excitation frequency. Hence, in this case, the system temperature at css will oscillate with frequency f equal to that imposed in the velocity boundary condition. However, this is an analytical method and the problem at hand is conjugate (i.e., solid-fluid coupled), so its applicability in the present context is doubtful.

## Some comparisons

I am naming the original method as marching and the new proposed method by Alex as direct. Using these methods, I ran some comparisons for f=2, u0=3 and q=5000/Ti, 1000/Ti (i.e., different heat inputs). I compared the cyclic average (mean) solid temperature along the line y=-e/2.

## For q=1000/Ti

As evident, there is some considerable difference among the direct and marching solutions. This discrepancy leads to the question as to, which is more accurate in modelling this phenomena?

• Please, pay attention that we can include standard deviations in our analysis as follows Tlist = Table[Tfs[i][0.5*L, -e/2]*Ti - 273.16, {i, nn - 49, nn}]; Tmean = Mean[Tlist];dT = StandardDeviation[Tlist] , therefore we have $43.1296\pm 2.69356$ as expected temperature at q=1000/Ti. Also, you computation for Tfs is not complete. May be it goes down after reaching some maximum. Commented Dec 16, 2022 at 6:38

Since I am the author of several books and articles in the field of turbulence, I prefer to use statistical methods. For this case, this method can be explained as follows. First, we calculate the hydrodynamic part of the problem. We assume that the phenomenon is periodically repeated. As can be seen from the above graph, the temperature fluctuates relative to the average value. We choose one period and average temperature over that period of time. The simplest algorithm is the following.

Needs["NDSolveFEM"]
{f = 2;
L = 0.040, d = 0.003, e = 0.005, kf = 0.026499, ks = 16,
rho = 1.1492, rhos = 7860, mu = 18.923*10^-6, cp = 1.0069*10^3,
cps = 502.4}; u0 = 3; nu = mu/rho; om = 2 Pi f;
tflow = 1;
t0 = .01;
NV = 2 f tflow;
nn = Round[NV \[Pi]/(om t0)]
Ti = 307; q = 1000/Ti;

reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 =
ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];

mesh = ToElementMesh[FullRegion[2], {{0, L}, {0, d}},
MaxCellMeasure -> 10^-7];
mesh1 = ToElementMesh[FullRegion[2], {{0, L}, {-e, d}},
MaxCellMeasure -> 10^-7];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 307/Ti; appro =
With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cps rhos + (cp rho - cps rhos) appro[y]);
eqs = {Inactive[
Div][({{-\[Mu], 0}, {0, -\[Mu]}} .
Inactive[Grad][u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] +
UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y] + (u[x, y] - UX[i - 1][x, y])/t0,
Inactive[
Div][({{-\[Mu], 0}, {0, -\[Mu]}} .
Inactive[Grad][v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] +
UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], y] + (v[x, y] - VY[i - 1][x, y])/t0,
D[u[x, y], x] + D[v[x, y], y]};
bc[i_] := {DirichletCondition[{u[x, y] == u0*Sin[om*i*t0],
v[x, y] == 0}, x == L (1 - Sign[Sin[om*i*t0]])/2 && 0 < y < d],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 0 || y == d],
DirichletCondition[p[x, y] == 0,
x == L (1 + Sign[Sin[om*i*t0]])/2 && 0 < y < d]};

Monitor[Do[{UX[i], VY[i], P[i]} =
NDSolveValue[{eqs == {0, 0, 0} /. \[Mu] -> nu, bc[i]}, {u, v,
p}, {x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];, {i, 1,
nn}], ProgressIndicator[i, {1, nn}]]; // AbsoluteTiming

Visualization of the flow for one period

Table[DensityPlot[UX[i][x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotRange -> All,
AspectRatio -> Automatic, Frame -> False], {i, nn - 49, nn}]

Monitor[Do[ux = If[y <= 0, 0, UX[i][x, y]];
vy = If[y <= 0, 0, VY[i][x, y]];
Tfs[i] =
NDSolveValue[{rde[y] ((ux*D[T[x, y], x] + vy*D[T[x, y], y])) -
Inactive[Div][
NeumannValue[q, y == -e],
DirichletCondition[{T[x, y] == 1},
x == L (1 - Sign[Sin[om*i*t0]])/2 && 0 <= y <= d]},
T, {x, y} \[Element] mesh1,
Method -> {"FiniteElement",
"InterpolationOrder" -> {T -> 2}}] // Quiet;, {i, nn - 49,
nn}], ProgressIndicator[i, {nn - 49, nn}]]; // AbsoluteTiming

Visualization of the mean temperature Mean temperature in a fixed point

Sum[Tfs[i][0.5*L, -e/2]*Ti - 273.16, {i, nn - 49, nn}]/50

Out[]= 43.1296

• It does not matter since $\rho$ is a constant, and we use homogeneous bc as DirichletCondition[p[x, y] == 0, x == L (1 + Sign[Sin[om*i*t0]])/2 && 0 < y < d]. Commented Dec 28, 2022 at 8:30