Modelling heat transfer in periodically reversing flow

Question

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["NDSolve`FEM`"]
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"];
quads = mesh["MeshElements"][[1]][[1]];
mark = Table[z = Mean[nodes[[quads[[i]]]]][[2]];
   If[z < 0, 0, 1], {i, 1, Length[quads]}];
MeshTotal1 = 
  ToElementMesh["Coordinates" -> nodes, 
   "MeshElements" -> {QuadElement[quads, mark]}];
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[
     Div][({{-(1/re), 0}, {0, -(1/re)}}.Inactive[Grad][
       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[
     Div][({{-(1/re), 0}, {0, -(1/re)}}.Inactive[Grad][
       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}] - 
   Inactive[Div][(ade[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[[3]][[1]], TFun["ValuesOnGrid"]}], {1, m, 
     10}]];
  m1 = If[k < K, n1 - 1, n1];
  AppendTo[SolutData1, 
   Take[Transpose[{UxFun[[3]][[1]], UxFun["ValuesOnGrid"]}], {1, m1, 
     10}]];
  m2 = If[k < K, n2 - 1, n2];
  AppendTo[SolutData2, 
   Take[Transpose[{UyFun[[3]][[1]], 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]

2. 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]

3. 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.

Answer 1

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["NDSolve`FEM`"]
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[[3]][[1]], 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}"}]