# Modelling heat transfer in periodically reversing flow

This is a heat transfer problem, which involves reciprocating (fully-reversing) fluid flow over a heated block of solid. The objective is to determine the temperature field in the solid and the fluid as the system reaches a quasi-steady state (i.e., temperature oscillates around a mean). I have asked a version of this question before here and here, and I have received excellent answers by Alex and Oleksii. However, I had some problems with grid-independence tests and less than expected temperature values, so I decieded to go ahead from the scratch and non-dimensionalize the equations:

The domain is $$X \in [0, L], Y \in [-e, d]$$ with heat flux $$q$$ applied at $$Y=-e$$. The solid extends from $$Y \in [-e, 0]$$, while the fluid domain is $$Y \in [0,d]$$. The fluid oscillates with $$U = U_0 \sin(\omega t)$$, where $$\omega = 2 \pi f$$. The dimensional temperature $$T^*$$ is non-dimensionalised as:

$$T = \frac{T^* - T^*_{inlet}}{\alpha}$$ where $$\alpha =\frac{qd}{k_s}$$. It must be noted that fresh fluid at some temperature $$T_{inlet} = 0$$ enters the domain in each half-cycle. For some simplicity, this $$T_{inlet}$$ can be assumed to be equal to the initial temperature $$T_{initial} = 0$$ of the system.

Now I want to solve a different set of non-dimensional equations describing the same problem. The non-dimensional scheme I used is $$u=U/u_0, x=X/d, y=Y/d, p=\frac{P}{\rho u_0^{2}}, \tau = \omega t$$, $$Re=\frac{\rho U_0 d}{\mu}$$.

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}=0 \tag 1$$

$$\frac{\omega d}{u_0}\frac{\partial u}{\partial \tau} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}+\frac{\partial p}{\partial x}-\frac{1}{Re}\big(\nabla^2 u\big)=0 \tag 2$$

$$\frac{\omega d}{u_0}\frac{\partial v}{\partial \tau} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y}+\frac{\partial p}{\partial y}-\frac{1}{Re}\big(\nabla^2 v\big)=0 \tag 3$$

$$\omega d^2 \frac{\partial T}{\partial \tau}+u_0 d\big(u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y}\big)-\frac{k}{\rho c_p}\big(\nabla^2 T\big)=0 \tag 4$$

The b.c. becomes:

$$u(\tau)= \sin(\tau) \tag 5$$ at $$x=0$$ and $$-\frac{\partial T}{\partial y} = 1 \tag 6$$ at $$y=-e/d$$.

Using the previous answers I received, I have tried to solve the above set of equations using the following code. I acknowledge that this framework of implementation was proposed by Oleksii, which I have modified. I have also borrowed concepts from Alex's answers:

Needs["NDSolveFEM"]
Needs["MeshTools"]

L = 0.040 ;(*length of the channel*)
d = 0.003;(*depth of the fluid*)
e = 0.005;(*depth of the solid*)
l = L/d;(*dimensionless length*)
rhof = 1.1492;(*fluid density*)
rhos = 7860;(*density of solid*)
mu = 18.923*10^-6;(*dynamic viscosity*)
nu = mu/rhof;(*kinematic viscosity*)
ks = 16;(*conductivity of solid*)
kf = 0.026499;(*conductivity of liquid*)
cf = 1069;(*heat capacity of fluid*)
cs = 502.4;(*heat capacity of solid*)
AlphaF = kf/(cf*rhof);(*thermal diffusivity of fluid*)
AlphaS = ks/(cs*rhos);(*thermal diffusivity of solid*)
f = 1.0;(*flow oscillation frequency*)
period = 1/f;(*period*)
omega = 2*Pi/period;(*circular frequency*)
u0 = 0.5;(*inflow velocity*)
q = 5000;(*heat flux density*)
Ti = 307;
re = d u0/(nu);
Pr = nu/AlphaF;(*Pandtl number*)

gamma = If[ElementMarker == 0, AlphaF/AlphaS, 1];
sigma = kf/ks;

(*Meshing*)
Nx = 30;(*number of elements in x-direction*)
NyF = 15;(*number of elements in y-direction in fluid*)
NyS = 5;(*number of elements in y-direction in solid*)
hy = 1./NyF;(*linear dimension of element in fluid*)
raster = {{{0, 0}, {l, 0}}, {{0, 1}, {l, 1}}};
MeshFluid = StructuredMesh[raster, {Nx, NyF}];(*FE mesh in fluid*)

raster = {{{0, -e/d}, {l, -e/d}}, {{0, 0}, {l, 0}}};
MeshSolid = StructuredMesh[raster, {Nx, NyS}];(*FE mesh in solid*)

mesh = MergeMesh[MeshSolid, MeshFluid];
nodes = mesh["Coordinates"];
If[z < 0, 0, 1], {i, 1, Length[quads]}];
MeshTotal1 =
ToElementMesh["Coordinates" -> nodes,
MeshTotal2 = MeshOrderAlteration[MeshTotal1, 2];

(*Incident veolcity profile*)
Clear[TopWall, BottomWall, reference, HeatInpBC, op, c, rampFunction,
sf, UinfProfile, Profile];

rampFunction[min_, max_, c_, r_] :=
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]
sf = rampFunction[0, 1, 0.25, 100];

Profile =
Interpolation[{{0, 0}, {hy, 1}, {1 - hy, 1}, {1, 0}},
InterpolationOrder -> 1];
Uc = 1/NIntegrate[Profile[y], {y, 0, 1}];(*calibration coefficient*)
UinfProfile[y_] := Uc*Profile[y];(*inflow velocity profile*)

(*Functions defining thermo-physical properties of solid and fluid. This allows solving a single energy equation*)
appro = With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cs rhos + (cf rhof - cs rhos) appro[y]);

(*PDE operator definitions. Sink term added to momentum equations to make velocity zero in the solid domain, which is supplied to the energy equation*)
c = If[ElementMarker == 0, 10^6,
0]; op = {{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][
u[t, x, y], {x, y}] +
Inactive[
u[t, x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(\[PartialD]$$, $${x}$$]$$p[t, x, y]$$\) +
c u[t, x, y] + ((omega d) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$u[t, x, y]$$\))/
u0, {{u[t, x, y], v[t, x, y]}}.Inactive[Grad][v[t, x, y], {x, y}] +
Inactive[
v[t, x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(\[PartialD]$$, $${y}$$]$$p[t, x, y]$$\) +
c v[t, x, y] + ((omega d) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$v[t, x, y]$$\))/u0, \!$$\*SubscriptBox[\(\[PartialD]$$, $${x}$$]$$u[t, x, y]$$\) + \!$$\*SubscriptBox[\(\[PartialD]$$, $${y}$$]$$v[t, x, y]$$\), (u0 d) {u[t, x, y], v[t, x, y]}.Inactive[Grad][
T[t, x, y], {x, y}] -
rde[y], {x, y}] + (omega d^2) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$T[t, x, y]$$\)};

(*Boundary conditions*)
TopWall =
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y == 1];
BottomWall =
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y <= 0];
reference = DirichletCondition[p[t, x, y] == 0., x == 0 && y == 0];
HeatInpBC = NeumannValue[(q d)/(ks), y == -(e/d)]


Assuming an initial and fluid inlet temperature of $$0$$, the following solves for the velocity and temperature fields:

Clear[UxLast, UyLast, TLast, PLast];
UxLast[x_, y_] := 0;
UyLast[x_, y_] := 0;
TLast[x_, y_] := 0;
PLast[x_, y_] := 0;
SolutData = {};
SolutData1 = {};
SolutData2 = {};
K = 10;(*number of half-periods considered*)
Monitor[Do[Clear[u, v, p, t, HeatDBC];
ti = (k - 1)*Pi;
tf = ti + Pi;
Clear[HeatDBC, Inflow, Outflow, bcs, ic, UxFun, UyFun, pressure,
TFun];
If[k == 1,
Inflow =
DirichletCondition[{u[t, x, y] == sf[t]*Sin[t]*UinfProfile[y],
v[t, x, y] == 0}, x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{u[t, x, y] == sf[t]*Sin[t]*UinfProfile[y],
v[t, x, y] == 0}, x == l && y > 0 && y < 1],
Inflow =
DirichletCondition[{u[t, x, y] == Sin[t]*UinfProfile[y],
v[t, x, y] == 0}, x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{u[t, x, y] == Sin[t]*UinfProfile[y],
v[t, x, y] == 0}, x == l && y > 0 && y < 1]];
If[OddQ[k] == True,
HeatDBC =
DirichletCondition[T[t, x, y] == 0, x == 0 && y >= 0 && y <= 1],
HeatDBC =
DirichletCondition[T[t, x, y] == 0, x == l && y >= 0 && y <= 1]];

ic = {u[ti, x, y] == UxLast[x, y], v[ti, x, y] == UyLast[x, y],
p[ti, x, y] == PLast[x, y], T[ti, x, y] == TLast[x, y]};

bcs = {TopWall, BottomWall, Inflow, Outflow, reference, HeatDBC};

{UxFun, UyFun, pressure, TFun} =
NDSolveValue[{op == {0, 0, 0, HeatInpBC}, bcs, ic}, {u, v, p,
T}, {x, y} \[Element] MeshTotal2, {t, ti, tf},
MaxStepSize -> 1*10^-2,
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"TemporalVariable" -> t,
"SpatialDiscretization" -> {"FiniteElement",
"PDESolveOptions" -> {"LinearSolver" -> "Pardiso"},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}];

UxLast =
ElementMeshInterpolation[{MeshTotal2}, Last[UxFun["ValuesOnGrid"]]];
UyLast =
ElementMeshInterpolation[{MeshTotal2}, Last[UyFun["ValuesOnGrid"]]];
TLast =
ElementMeshInterpolation[{MeshTotal2}, Last[TFun["ValuesOnGrid"]]];
PLast =
ElementMeshInterpolation[{MeshTotal1},
Last[pressure["ValuesOnGrid"]]];
n = Length[TFun["ValuesOnGrid"]];
n1 = Length[UxFun["ValuesOnGrid"]];
n2 = Length[UyFun["ValuesOnGrid"]];
m = If[k < K, n - 1, n];
AppendTo[SolutData,
Take[Transpose[{TFun[][], TFun["ValuesOnGrid"]}], {1, m,
10}]];
m1 = If[k < K, n1 - 1, n1];
AppendTo[SolutData1,
Take[Transpose[{UxFun[][], UxFun["ValuesOnGrid"]}], {1, m1,
10}]];
m2 = If[k < K, n2 - 1, n2];
AppendTo[SolutData2,
Take[Transpose[{UyFun[][], UyFun["ValuesOnGrid"]}], {1, m2,
10}]];, {k, 1, K}], ProgressIndicator[k, {1, K}]]

(*Generating solution functions*)
Clear[TsolVec, TFun]
TsolVec =
Interpolation[Flatten[SolutData, 1], InterpolationOrder -> 1];
TFun[t_?NumericQ] :=
ElementMeshInterpolation[{MeshTotal2}, TsolVec[t]]

Clear[UxsolVec, UxFun]
UxsolVec =
Interpolation[Flatten[SolutData1, 1], InterpolationOrder -> 1];
UXFun[t_?NumericQ] :=
ElementMeshInterpolation[{MeshTotal2}, UxsolVec[t]]

Clear[UysolVec, UyFun]
UysolVec =
Interpolation[Flatten[SolutData2, 1], InterpolationOrder -> 1];
UYFun[t_?NumericQ] :=
ElementMeshInterpolation[{MeshTotal2}, UysolVec[t]]


I then plotted the temperature history at a point in the solid and temperature profile in the solid. These results look qualitatively correct but their magnitudes are far blown up:

1. Temperature history for 70 half-periods
Plot[(TFun[t][0.5 l, -e/(2 d)]), {t, 0, K*Pi}, GridLines -> Automatic, PlotRange -> Full] 1. Cyclic average temperature profile in the solid
Tsm[x_] = (2 π)^-1 (Sum[TFun[t][x, -e/d/2], {t, (K - 2) π, K π}]);
Plot[Tsm[x], {x, 0, L/d}, GridLines -> Automatic] 1. I plotted the time variation of the $$x-$$velocity at the channel mid and found no unreasonable values
Plot[{UXFun[t][l/2, 1/2]}, {t, (K - 2) Pi, K Pi}, GridLines -> Automatic] This implies that there must be something wrong with the way I am implementing the energy equation and its boundary conditions. However, I have not been able to figure out what.

• Code consists of several typos. I can't reproduce solution. But it looks like you need to make this correction HeatInpBC = NeumannValue[(q d)/(cs rhos), y == -(e/d)]; Jan 11 at 13:03
• Have edited code, you can try now. Will try this edit suggested by you. However, as per the derivation $-k_s \frac{\partial T}{\partial y} = q$, then how is the Neumann value $\frac{qd}{c_s \rho_s}$ ? Jan 11 at 13:32
• @Avrana I will try to figure out tomorrow in this code. It is blackout in my city now. It seems that something wrong in coefficients Jan 11 at 14:10
• @Avrana Please, look at your equation for T where you use Inactive[Div][(ade[y] Inactive[Grad][T[t, x, y], {x, y}])/rde[y], {x, y}]. In FEM algorithm it generates bc as -ade[-e]/rde[-e] D[T,y]=qd/ks, and it is why you have temperature about 10^6. Normally it should be -D[T,y]=qd/ks, To equalize sides you need to multiply on ade[-e]/rde[-e]=ks/(cs rhos). Jan 11 at 14:27
• If you suspect an issue with the heat equation you could use HeatTransferPDEComponent in stead or to double check what you have. But if the scale is wrong, does not not suggest an error in the non dimensionalization? Jan 11 at 21:10

Let's the velocity, temperature and pressure are measured in units $$u_0,\, dq/k_s,\, \rho u_0^2$$ respectively, time and space coordinates are measured in units $$d/u_0$$ and $$d$$. In this case the governing equations in dimensionless form are as follows:

Navier-Stokes equations:

$$\begin{equation} \frac{\partial \vec{V}}{\partial t}+ (\vec{V}\cdot\nabla)\vec{V}=-\nabla P+\frac{1}{Re}\Delta \vec{V}-C\cdot\vec{V} \end{equation}$$

$$\begin{equation} \nabla\cdot \vec{V}=0 \end{equation}$$

Energy conservation:

$$\begin{equation} \gamma Pe\left(\frac{\partial T}{\partial t}+(\vec{V}\cdot\nabla)T \right)=\nabla\cdot\left(\gamma\nabla T\right) \end{equation}$$

where $$Re=u_0d/\nu$$ is the Reinolds number, $$Pe$$ is the variable in space Peclet number $$\begin{equation} Pe=\begin{cases} u_0d/\alpha_s, & \{x,y\}\in solid \\ u_0d/\alpha_f, & \{x,y\}\in fluid \\ \end{cases} \end{equation}$$ Coefficient $$C$$ in penalty term is as follows: $$\begin{equation} C=\begin{cases} 10^6, & \{x,y\}\in solid \\ 0, & \{x,y\}\in fluid \\ \end{cases} \end{equation}$$ Coefficient $$\gamma$$:

$$\begin{equation} \gamma=\begin{cases} 1, & \{x,y\}\in solid \\ k_f/k_s, & \{x,y\}\in fluid \\ \end{cases} \end{equation}$$

Equations contain variable in space coefficients. Solution of such PDE are described here. In OP the temperature on inflow boundary is changed rapidly at the beginning of every half-period. The boundary conditions at this moment are not consistent with initial conditions here and calculated temperature on inlet can differ from $$0$$. I propose to change the temperature on inlet gradually up to 0 during the time which is small compared with period of oscillation.

Input parameters and mesh generation:

Needs["NDSolveFEM"]
Needs["MeshTools"]

L = 0.040;(*length of the channel*)
d = 0.003;(*depth of the fluid*)
e = 0.005;(*depth of the solid*)
l = L/d;(*dimensionless length*)
rhof = 1.1492;(*fluid density*)
rhos = 7860;(*density of solid*)
mu = 18.923*10^-6;(*dynamic viscosity*)
nu = mu/rhof;(*kinematic viscosity*)
ks = 16;(*conductivity of solid*)
kf = 0.026499;(*conductivity of liquid*)
cf = 1069;(*heat capacity of fluid*)
cs = 502.4;(*heat capacity of solid*)
AlphaF = kf/(cf*rhof);(*thermal diffusivity of fluid*)
AlphaS = ks/(cs*rhos);(*thermal diffusivity of solid*)
f = 1.0;(*flow oscillation frequency*)
period = 1/f;(*period*)
omega = 2*Pi/period;(*circular frequency*)
u0 = 0.5;(*inflow velocity*)
q = 5000;(*heat flux density*)
Ti = 307; (*inflow temperature*)
re = d *u0/nu; (*reinolds number*)
gamma = If[y < 0, 1, kf/ks]; (*relation of conductivities*)
Pe = If[y < 0, u0*d/AlphaS, u0*d/AlphaF];  (*Peclet number*)
c = If[y < 0, 10^6, 0];(*constant in momentum sink term*)

Nx = 50;(*number of elements in x-direction *)
NyF = 5;(*number of elements in y-direction in fluid*)
NyS = 5;(*number of elements in y-direction in solid*)
hy = 1./NyF;(*linear dimension of element in fluid*)
raster = {
{{0, 0}, {l, 0}},
{{0, 1}, {l, 1}}
};
MeshFluid = StructuredMesh[raster, {Nx, NyF}];(*FE mesh in fluid*)
raster = {
{{0, -e/d}, {l, -e/d}},
{{0, 0}, {l, 0}}
};
MeshSolid = StructuredMesh[raster, {Nx, NyS}];(*FE mesh in solid*)
MeshTotal1 = MergeMesh[MeshSolid, MeshFluid];
MeshTotal2 = MeshOrderAlteration[MeshTotal1, 2];
Show[MeshTotal2["Wireframe"], ImageSize -> 600]


Implementation of PDE and BC

Clear[TopWall, BottomWall, reference, HeatInpBC, op, rampFunction, sf,
UinfProfile, Profile, x, y, t];
rampFunction[min_, max_, c_, r_] :=
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]
sf = rampFunction[0, 1, 0.25, 50];
Profile =
Interpolation[{{0, 0}, {hy, 1}, {1 - hy, 1}, {1, 0}},
InterpolationOrder -> 1]
UinfProfile[y_] := Profile[y]/NIntegrate[Profile[y], {y, 0, 1}]

op = {
D[u[t, x, y], t] +
Inactive[
Div][({{-1/re, 0}, {0, -1/re}} .
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}] + D[p[t, x, y], x] +
c*u[t, x, y],
D[v[t, x, y], t] +
Inactive[
Div][({{-1/re, 0}, {0, -1/re}} .
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}] + D[p[t, x, y], y] +
c*v[t, x, y], D[u[t, x, y], x] + D[v[t, x, y], y],
Pe*gamma*D[T[t, x, y], t] +
Pe*gamma*{{u[t, x, y], v[t, x, y]}} .
Inactive[Grad][T[t, x, y], {x, y}] +
Inactive[Div][{{-gamma, 0}, {0, -gamma}} .
Inactive[Grad][T[t, x, y], {x, y}], {x, y}]

};

TopWall =
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y == 1];
BottomWall =
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, y <= 0];
(*setting the pressure value in single node*)
reference = DirichletCondition[p[t, x, y] == 0., x == 0 && y == 0];
HeatInpBC = NeumannValue[1, y == -e/d];


Solution of PDE

Here only 25 periods are considered

Clear[UxLast, UyLast, TLast, PLast];
UxLast[x_, y_] := 0;
UyLast[x_, y_] := 0;
TLast[x_, y_] := 0.;
PLast[x_, y_] := 0;
SolutData = {};
K = 50;(*number of half period considered*)

Do[
Clear[u, v, p, t, HeatDBC];
ti = (k - 1)*0.5 period*u0/d;
tf = ti + 0.5 period*u0/d;

Clear[HeatDBC, Inflow, Outflow, bcs, ic, UxFun, UyFun, pressure,
TFun];
If[k == 1,
Inflow =
DirichletCondition[{u[t, x, y] ==
sf[t*d/u0]*Sin[t*(omega*d)/u0]*UinfProfile[y], v[t, x, y] == 0},
x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{u[t, x, y] ==
sf[t*d/u0]*Sin[t*(omega*d)/u0]*UinfProfile[y], v[t, x, y] == 0},
x == l && y > 0 && y < 1],

Inflow =
DirichletCondition[{u[t, x, y] ==
Sin[t*(omega*d)/u0]*UinfProfile[y], v[t, x, y] == 0},
x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{u[t, x, y] ==
Sin[t*(omega*d)/u0]*UinfProfile[y], v[t, x, y] == 0},
x == l && y > 0 && y < 1]
];

(*temperature on inlet changes gradually up to 0 during dt*)
dt = 0.01*0.5 period*u0/d;
If[OddQ[k] == True,
HeatDBC =
DirichletCondition[
T[t, x, y] ==
If[t <= ti + dt, TLast[0, y] - TLast[0, y]*(t - ti)/dt, 0],
x == 0 && y > 0 && y <= 1],
HeatDBC =
DirichletCondition[
T[t, x, y] ==
If[t <= ti + dt, TLast[l, y] - TLast[l, y]*(t - ti)/dt, 0],
x == l && y > 0 && y <= 1]
];

ic = {u[ti, x, y] == UxLast[x, y], v[ti, x, y] == UyLast[x, y],
p[ti, x, y] == PLast[x, y], T[ti, x, y] == TLast[x, y]};
bcs = {TopWall, BottomWall, Inflow, Outflow, reference, HeatDBC};

Monitor[
{UxFun, UyFun, pressure, TFun} =
NDSolveValue[{op == {0, 0, 0, HeatInpBC}, bcs, ic}, {u, v, p,
T}, {x, y} \[Element] MeshTotal2, {t, ti, tf},

MaxStepSize -> 5*10^-3*0.5 period*u0/d,
Method -> {

"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2},

"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"PDESolveOptions" -> {"LinearSolver" -> "Pardiso"},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}
, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])]
, currentTime];

UxLast =
ElementMeshInterpolation[{MeshTotal2}, Last[UxFun["ValuesOnGrid"]] ];
UyLast =
ElementMeshInterpolation[{MeshTotal2},
Last[UyFun["ValuesOnGrid"]]];
TLast =
ElementMeshInterpolation[{MeshTotal2}, Last[TFun["ValuesOnGrid"]]  ];
PLast =
ElementMeshInterpolation[{MeshTotal1},
Last[pressure["ValuesOnGrid"]]  ];

n = Length[TFun["ValuesOnGrid"]];
m = If[k < K, n - 1, n];
AppendTo[SolutData,

Take[Transpose[{TFun[][], TFun["ValuesOnGrid"]}], {1, m, 10}]
]

, {k, 1, K}
]


Postprocessing

Clear[TsolVec, TFun]
TsolVec =
Interpolation[Flatten[SolutData, 1], InterpolationOrder -> 1];
TFun[t_?NumericQ] :=
ElementMeshInterpolation[{MeshTotal2}, TsolVec[t]]


Dynamics of temperature in point $$\{L/2,d/2\}$$ looks as follows:

pic1=Plot[Ti + (d*q)/ks TFun[t*u0/d][l/2, 0.5], {t, 0,K*0.5*period},
PlotStyle -> {Thickness[0.003], RGBColor[0, 0, 0]},
PlotRange -> All, Frame -> True,
FrameLabel -> {"time", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 500,
LabelStyle -> RGBColor[0, 0, 0]]; Let's look at temperature dynamics in points $$\{0,d/2\}$$ and $$\{L,d/2\}$$ during first 5 periods:

pic2 = Plot[{Ti + (d*q)/ks TFun[t*u0/d][0, 0.5],
Ti + (d*q)/ks TFun[t*u0/d][l, 0.5]}, {t, 0, 5*period},
PlotStyle -> {{Thickness[0.003],
RGBColor[0, 0, 0]}, {Thickness[0.003], RGBColor[1, 0, 0]}},
PlotRange -> All, Frame -> True,
FrameLabel -> {"time [s]", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 500,
LabelStyle -> RGBColor[0, 0, 0],
PlotLegends -> {"{0,0.5d}", "{L,0.5d}"}] • Thankyou for the well explained answer. I see that you have managed to remove the discrepancy regarding the Inflow b.c. in this version. I will go through the solution diligently and get back. Jan 14 at 2:55
• I have run the code and will be running some more tests. It certainly satisfies the b.c.(s) in each period. Meanwhile, I have some queries. $(1)$ I see that you have set up the non-dimensional equations using the $Pe$ number, which is treated as a space variable coefficient (SVC). As a result, you had to introduce $\gamma$, which again is a SVC. Why could we just not use $\frac{k}{\rho c_p}$ as the SVC ? Jan 14 at 9:10
• I have changed a little bit the code. Jan 18 at 16:42
• No, I have just simplified the code. I did not studied the optimal time step. The computations of course are time consuming Jan 18 at 16:54
• @Avrana Hello! I will look tomorrow Jun 6 at 19:10