# Reciprocating flow in a channel over a heated surface

The following is a coupled heat transfer and fluid flow problem.

A thick plane channel is being heated with a constant flux from the bottom (at $$y=-e$$) with a constant heat flux $$q$$ as shown in the attached figure. The heated channel is subjected to a reciprocating flow with velocity $$U(t)=u_{max}\sin(2\pi ft)$$. I must state here that the flow velocity has a mean of $$0$$ which means for the first half of the cycle it reaches $$u_{max}$$ and in the second half reaches $$-u_{max}$$. In the expression for velocity the term $$f$$ stands for freqency of the oscillating flow. Thus, the time period of the flow is $$\tau=1/f$$. This indicates that for $$0, the boundary at $${x=0, 0 acts as the inlet, while for the second half, i.e., $$\tau>t>\tau/2$$ the boundary at $${x=L, d>y>0}$$ is the inlet. The top plane ($$y=d, 0), the left solid face ($$x=0, 0>y>-e$$) and the right solid face ($$x=L, -e) are all insulated. The solid and fluid domains are coupled through the temperature and flux continuity at the interface ($$y=0$$).

In this scenario, after a certain time the entire system is supposed to reach a cyclic steady state. I will now mention the governing equations and boundary conditions:

Fluid

Hydrodynamic

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag1$$ $$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \frac{-1}{\rho} \frac{\partial p}{\partial x} + \mu (\nabla^2 u) \tag2$$ $$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \frac{-1}{\rho} \frac{\partial p}{\partial y} + \mu (\nabla^2 v) \tag3$$

No slip and No penetration condition on $$y=0$$ and $$y=d$$.

Thermal

$$\rho c_p \frac{\partial T}{\partial t} + \rho c_p u \frac{\partial T}{\partial x} + \rho c_p v \frac{\partial T}{\partial y}= k_f (\nabla^2 T) \tag4$$

The fluid has an inlet temperature $$T=T_i$$. So the boundary condition will be $$T(x=0)=T_i, \frac{\partial T}{\partial x} \vert_{x=L} = 0$$ for $$\tau/2>t>0$$ and $$T(x=L)=T_i, \frac{\partial T}{\partial x} \vert_{x=0} = 0$$ for $$\tau/2<0<\tau$$

Solid

Thermal

$$\rho_s c_{p,s} \frac{\partial T_s}{\partial t} = k_s \nabla^2 T_s \tag 5$$

The boundary conditions are: $$\frac{\partial T_s}{\partial t}\vert_{x=0}=\frac{\partial T_s}{\partial t}\vert_{x=L} = 0$$ and $$-k_s\frac{\partial T_s}{\partial t} = q$$.Here, $$k_s$$ is the solid thermal conductivity.

Coupling

At the interface between the solid and the fluid, the following holds which couples the problem:

Continuity of Temperature $$T(x,0)= T_s(x,0), 0 Continuity of Flux $$-k_s \frac{\partial T_s}{\partial y}\vert_{y=0,x} = -k_f \frac{\partial T}{\partial y}\vert_{y=0,x} \tag7$$

The objective is to solve for the velocity and the temperature fields in the solid and fluid. Also, the flux transfer from the solid to the fluid at the interface is of interst. I have found a Mathematica fluid solver here, but it only simulates isothermal flows:

(1) Are there non-isothermal flow mathematica solvers ?

(2) How should I model the coupling between the solid and fluid using (6) and (7) ? In a time-step what should be calculated first, the solid temperature or the fluid? I do understand that both the fields have to be calculated simultaneously.

(3) Finally, since the flow is reversing, how should I model the switching of boundary conditions during each half cycle ?

Some parameter values L=0.025, d=0.002, e=0.002, k_f=0.614, k_s=390, q=5000, rho=997, rhos=8960, mu=8.90*10^-4, cp =4178, cps =385. These parameters represent flow of water over copper. A typical flow velocity profile can be U=0.3*sin(2*pi*1*t), which is a 1Hz flow giving tau=1.

A typical CFD result obtained from COMSOL

For the same material combinations but with a e=4, f=0.5, umax=0.22875, Ti=288 and a heat input of $$1W$$ which translates to q=40000, I attach the average interface ($$y=0$$) temperature (line integrated along the length from $$x=0$$ to $$x=L$$) variation with time.

• Could you provide $c_p,c_{ps}$ as well? Feb 2, 2022 at 12:20
• Have you seen the examples from the documentation? For example the Heat Exchanger or the Buoyancy-Driven Flow Feb 2, 2022 at 14:02
• @AlexTrounev I have added the values Feb 2, 2022 at 14:27
• @Avrana I'am confused a little about BC for eq.(4) on $x=0$ and $x=L$. Such a conditions should lead to sharp temperature drop (growth) on these boundaries. How to imagine it it physically? Feb 3, 2022 at 14:13
• @Avrana See my answer. Please note, that we don't need any boundary condition for temperature on interface since we solve problem with FEM in the region $0\le x\le L, -e\le y\le d$. Feb 3, 2022 at 15:23

We can solve this problem with method proposed on my page.

Solution1. We use nondimensional form of equations with scale d and $$t_s = d^2/(k_f/(c_p \rho))$$. We define two regions reg1, reg2 to describe fluid flow and temperature consequently. This is code in a case of divergent form of temperature equation. We use scaled form of heat flux qn = q ts/(cp rho)/d, but temperature is unscaled

{f = 1; L = 0.025, d = 0.002, e = 0.002, kf = 0.614, ks = 390,
q = 5000, rho = 997, rhos = 8960, mu = 8.90*10^-4,
cp = 4.178*10^3 (*J/kg/\[Degree]K*), cps = 385}; Pr =
mu/rho/(kf/(cp rho));
Pr0 = Pr;   a = ks/kf; ts = d^2/(kf/(cp rho)); as =
ts ks/(cp rho)/d^2; rs = cps rhos/(cp rho); u0 = .25 ts/d; om =
2 Pi f ts; t0 = 1/om/10; nn = Round[4 Pi/(t0 om)]; qn =
q ts/(cp rho)/d;
Needs["NDSolveFEM"]

reg1 = ImplicitRegion[0 <= x <= L/d && 0 <= y <= 1, {x, y}]; reg2 =
ImplicitRegion[0 <= x <= L/d && -e/d <= y <= 1, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 0; appro =
With[{k = 2. 10^4}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (as + (1 - as) UnitStep[y] /. UnitStep -> appro);
rde[y_] := (rs + (1 - rs) UnitStep[y] /. UnitStep -> appro);

Do[
{UX[i], VY[i], P[i]} =
NDSolveValue[{{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]} == {0, 0, 0} /. \[Mu] ->
Pr0, {

DirichletCondition[{u[x, y] == u0*Sin[om*i*t0]*y*(1 - y),
v[x, y] == 0}, x == 0],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0},
y == 0 || y == 1]},
DirichletCondition[p[x, y] == 0, x == L/d]}, {u, v,
p}, {x, y} \[Element] reg1,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}];
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[qn, y == -e/d ],
DirichletCondition[{T[x, y] == 0},

x == 0 + L/d (1 - Sign[Sin[om*i*t0]])/2 && 0 <= y <= 1]},
T, {x, y} \[Element] reg2,
Method -> {"FiniteElement",
"InterpolationOrder" -> { T -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.01}}] //
Quiet;, {i, 1, nn}] // AbsoluteTiming


Visualization of temperature

Table[DensityPlot[Tfs[i][x, y], {x, y} \[Element] mesh1,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
AspectRatio -> 1/2, PlotLabel -> i ts t0, PlotRange -> All,
PlotPoints -> 50], {i, 5, nn, 10}]


Solution 2. We use unscaled version of code with given input parameters and boundary conditions as in a paper Oscillatory Heat Transfer in a Pipe Subjected to a Laminar Reciprocating Flow by T. S. Zhao & P. Cheng

{f = 1; L = 0.025, d = 0.002, e = 0.002, kf = 0.614, ks = 390,
q = 5000, rho = 997, rhos = 8960, mu = 8.90*10^-4,
cp = 4.178*10^3 (*J/kg/\[Degree]K*), cps = 385}; u0 = .3; nu =
mu/rho; om = 2 Pi f ; t0 = .1; nn = Round[10 Pi/(om t0)];

Needs["NDSolveFEM"]

reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 =
ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 0; appro =
With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) UnitStep[y] /. UnitStep -> appro);
rde[y_] := (cps rhos + (cp rho - cps rhos) UnitStep[y] /.
UnitStep -> appro);

Do[
{UX[i], VY[i], P[i]} =
NDSolveValue[{{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]} == {0, 0, 0} /. \[Mu] ->
nu, {

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]}, {u, v,
p}, {x, y} \[Element] reg1,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.00000005}}];
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] == 0},
x == L (1 - Sign[Sin[om*i*t0]])/2 && 0 <= y <= d]},
T, {x, y} \[Element] reg2,
Method -> {"FiniteElement",
"InterpolationOrder" -> { T -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}] //
Quiet;, {i, 1, nn}] // AbsoluteTiming



Visualization of temperature vs time at x=L/2, y=0 and x=L/2 for different t=1,2,3,4,5 s

ListLinePlot[Table[{i t0, Tfs[i][L/2, 0]}, {i, 0, nn}],
AxesLabel -> {"t, s", "T"}]

Plot[Evaluate[Table[Tfs[i][L/2, y], {i, 10, nn, 10}]], {y, -e, d},
PlotLegends -> Table[i t0, {i, 10, nn, 10}]]


Temperature and velocity distributions for different time shown above

Note, that temperature distribution is same for scaled and unscaled form of equation. Let consider test example posted by Avrana and solved with COMSOL. To get same average temperature we need to increase q up to q=120000, then we have

Needs["NDSolveFEM"]

{f = .5, L = 0.025, d = 0.001, e = 0.004, kf = 0.614, ks = 390,
rho = 997, rhos = 8960, mu = 8.90*10^-4,
cp = 4.178*10^3 (*J/kg/\[Degree]K*), cps = 385}; u0 = 0.22875; nu =
mu/rho; om = 2 Pi f ; t0 = .1; nn =
Round[40 Pi/(om t0)]; Ti = 288; q = 120000/Ti;

reg1 = ImplicitRegion[0 <= x <= L && 0 <= y <= d, {x, y}]; reg2 =
ImplicitRegion[0 <= x <= L && -e <= y <= d, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;
Tfs[0][x_, y_] := 293/Ti; appro =
With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) UnitStep[y] /. UnitStep -> appro);
rde[y_] := (cps rhos + (cp rho - cps rhos) UnitStep[y] /.
UnitStep -> appro);

Do[
{UX[i], VY[i], P[i]} =
NDSolveValue[{{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]} == {0, 0, 0} /. \[Mu] ->
nu, {

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]}, {u, v,
p}, {x, y} \[Element] reg1,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.00000005}}];
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] reg2,
Method -> {"FiniteElement",
"InterpolationOrder" -> { T -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0000001}}] //
Quiet;, {i, 1, nn}];


Average temperature

ff = Interpolation[
Table[{i t0,
Ti/L NIntegrate[Tfs[i][x, 0], {x, 0, L}, AccuracyGoal -> 4,
PrecisionGoal -> 4]}, {i, 0, nn}], InterpolationOrder -> 2]

Plot[ff[t], {t, 0, 40}, AxesLabel -> {"t, s", "T"}, PlotRange -> All]


• Appreciate the very detailed answer. I will now go through it and come back with any questions. Thanks again. Feb 3, 2022 at 16:47
• @Avrana Please, check that we can use DirichletCondition[{u[x, y] == u0*Sin[om*i*t0], v[x, y] == 0}, x == 0] as well with same result for temperature. Feb 4, 2022 at 6:02
• I am going through the solution till this time. Actually, i had non dimensionalized using a different scheme as described in a paper I had linked to in one of my comments to the question. So I am rewriting the equations now using your approach to non-dimensionalization. I will surely run with the condition you just mentioned and get back. Feb 4, 2022 at 6:13
• @Avrana See Update 2 with unscaled version of code and with bc from Zhao paper. Feb 11, 2022 at 8:59
• @Avrana We use scales $d$ and $t_s = d^2/(k_f/(c_p \rho))$. With this scaling we have \[Mu] -> Pr0. This scaling used in convection problem - see, for instance, Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264. Jan 4, 2023 at 14:32

It seems that the main challenge in this problem is Dirichlet BC which should be switched periodically on $$x=0$$ and $$x=L$$. I don't know whether it possible to switch BC inside NDSolveValue but we can run NDSolveValue every half period with new BC and take use solution obtained at last half period as a initial conditions. Velocity field is not influenced by temperature so that partitioned coupling scheme can be used i.e. at first stage velocity is calculated than temperature is obtained. But in this case interpolation functions for velocity are involved into calculations inside the solver that significantly decelerates calculation. I propose to use here monolithic approach which implies calculation both temperature and velocity in single code. We will solve NS equations in computational domain which includes solid and fluid. By introducing the momentum sink term $$-C\cdot \vec{V}$$ into the momentum conservation equations one can set to zero the velocity in solid phase. Here $$C$$ is a large number. In current simulation the value $$C=10^6$$ was used. Lets write the governing equations in dimensionless form as in the paper [Zhao1996] which is mentioned by @Avrana in comments. Velocity and pressure are measured in units $$u_0$$ and $$u_0^2\rho$$, spatial coordinates and time are dimentionelized by $$d$$, $$\omega^{-1}$$ respectively. Here $$\omega$$, $$u_0$$ are circular frequency and inflow velocity accordingly. Temperature is measured in units $$qd/k_f$$, where $$q$$ is a heat flux density supplied from bottom wall.

Navier-Stokes equations

$$$$\frac{\partial \vec{V}}{\partial t}+\frac{A_0}{2}\left[ (\vec{V}\cdot\nabla)\vec{V} +\nabla P\right ]=\frac{1}{Re_{\omega}}\Delta \vec{V}$$$$

$$$$\nabla\cdot \vec{V}=0$$$$

Energy conservation in fluid $$$$\frac{\partial T}{\partial t}+\frac{A_0}{2}(\vec{V}\cdot\nabla)T=\frac{1}{Re_{\omega}Pr}\Delta T$$$$

Energy conservation in solid $$$$\frac{\partial T}{\partial t}=\frac{1}{\Gamma Pr Re_{\omega}}\Delta T$$$$

where $$Re_{\omega}=\omega d^2/\nu$$, $$A_0=2u_0/(d\omega)$$, $$Pr=\nu/\alpha_f$$, $$\Gamma=\alpha_f/\alpha_s$$ are dimensionless parameters.

Input parameters

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

L = 0.025;(*length of the channel*)
d = 0.002;(*depth of the fluid*)
e = d;(*depth of the solid*)
l = L/d; (*dimensionless length*)
rhof = 997;(*fluid density*)
rhos = 8960;(*density of solid*)
mu = 8.9*10^-4;(*dynamic viscosity*)
nu = mu/rhof;(*kinematic viscosity*)
ks = 390;(*conductivity of solid*)
kf = 0.614;(*conductivity of liquid*)
cf = 4178;(*heat capacity of fluid*)
cs = 385;(*heat capacity of solid*)
AlphaF =kf/(cf*rhof); (*thermal diffusivity of fluid*)
AlphaS = ks/(cs*rhos); (*thermal diffusivity of solid*)
period = 1.;(*period*)
omega = 2*Pi/period;(*circular frequency*)
u0 = 0.3;(*inflow velocity*)
q = 5000;(*heat flux density*)

(*dimensionless model input parameters *)
A0 = 2*u0/(d*omega);
re = omega*d^2/(nu);
Pr = nu/AlphaF;(*Pandtl number*)
gamma=If[ElementMarker == 0, AlphaF/AlphaS, 1];
sigma = kf/ks;


FE mesh generation

It is convenient to use structured FE mesh for this particular case. All the manipulations with meshes were done here by means of utilities from MeshTools package.

Nx = 100;(*number of elements in x-direction *)
NyF = 20;(*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"];
(*ElementMarker=0 in soilid and 1 in fluid*)

If[z < 0, 0, 1], {i, 1, Length[quads]}];
(*1d order mesh in total domain*)

MeshTotal1 =
ToElementMesh["Coordinates" -> nodes,
(*2d order mesh in total domain*)

MeshTotal2 = MeshOrderAlteration[MeshTotal1, 2];
pic1 = Show[MeshTotal1["Wireframe"], ImageSize -> 600]
Export["pic1.jpeg", pic1, ImageResolution -> 300]


Solution procedure

NS solver used little differ from those in documentation. Function rampFunction introduced there helps to increase inflow velocity smoothly in time. Velocity inflow profile  UinfProfile[y] used here has a trapezoidal shape.

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*)


Define a PDE operator with boundary conditions

c = If[ElementMarker == 0, 10^6,
0];(*define the constant in momentum sink term*)
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}] +
0.5 A0*{{u[t, x, y], v[t, x, y]}} .
Inactive[Grad][u[t, x, y], {x, y}] + 0.5 A0*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}] +
0.5 A0*{{u[t, x, y], v[t, x, y]}} .
Inactive[Grad][v[t, x, y], {x, y}] + 0.5 A0*D[p[t, x, y], y] +
c*v[t, x, y],
D[u[t, x, y], x] + D[v[t, x, y], y],
D[T[t, x, y], t] +
Inactive[
Div][(-(1/(Pr*re*gamma))*
Inactive[Grad][T[t, x, y], {x, y}]), {x, y}] +
0.5*A0*{u[t, x, y], v[t, x, y]} . Inactive[Grad][T[t, 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 pressure value in single node*)

reference = DirichletCondition[p[t, x, y] == 0., x == 0 && y == 0];
HeatInpBC = NeumannValue[sigma*AlphaS/(AlphaF*Pr*re), y == -1];


Finally, the next code realizes the solution procedure for first twenty half-periods

Clear[UxLast, UyLast, TLast, PLast];
UxLast[x_, y_] := 0;
UyLast[x_, y_] := 0;
TLast[x_, y_] := 0.;
PLast[x_, y_] := 0;

SolutData = {};
K = 20;(*number of half-periods considered*)
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};

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 -> Pi*10^-3,

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

Construction of interpolation function for temperature solution

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


Visualization of thermal history in point $$\{0.5L,0\}$$

Plot[TFun[t][0.5 l, 0], {t, 0, K*Pi},
PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 0]},
PlotRange -> {0, 0.02}, Frame -> True,
FrameLabel -> {"time", "Temperature"},
FrameTicks -> {Table[(k - 1) 2 \[Pi], {k, 1, 6}], Automatic},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 500,
LabelStyle -> RGBColor[0, 0, 0],
PlotLabel -> "Thermal history in point {0.5l,0}"]


It is necessary to calculate at least 5 periods to reach quasi-stationary regime under the given input parameters (properties, geometry, inflow velocity) as we can see from the last pic. Distributions of temperature in cross section $$x=0.5L$$ at different times looks as follows

t1 = 2 Pi; t2 = 4 Pi;
t3 = 6 Pi; t4 = 8 Pi;
t5 = 10 Pi; t6 = 18 Pi;

Plot[{TFun[t1][0.5 l, y], TFun[t2][0.5 l, y],
TFun[t3][0.5 l, y], TFun[t4][0.5 l, y], TFun[t5][0.5 l, y],
TFun[t6][0.5 l, y]}, {y, -1, 1}, PlotStyle -> Thickness[0.005],
PlotRange -> All, Frame -> True, FrameLabel -> {"y", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 500,
PlotLegends -> {2 \[Pi], 4 \[Pi], 6 \[Pi], 8 \[Pi], 10 \[Pi],
18 \[Pi]},
PlotLabel ->
"Temperature distribution along the line x=0.5L at different \
times", LabelStyle -> RGBColor[0, 0, 0]]


Multiplication $$T\cdot qd/k_f$$ with $$q=5\cdot 10^3 W/m^2$$ gives overheating (in Kelvins) above initial temperature of supplied water $$T_0$$.

Plot[{(q*d)/kf*TFun[t1][0.5 l, y/d], (q*d)/kf*TFun[t2][0.5 l, y/d], (
q*d)/kf*TFun[t3][0.5 l, y/d], (q*d)/kf*TFun[t4][0.5 l, y/d], (q*d)/
kf*TFun[t5][0.5 l, y/d], (q*d)/kf*TFun[t6][0.5 l, y/d]}, {y, -e,
d}, PlotStyle -> Thickness[0.004], PlotRange -> All, Frame -> True,
FrameLabel -> {"y, m", "T-\!$$\*SubscriptBox[\(T$$, $$0$$]\), K"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 500,
PlotLegends -> {"t=1s", "t=2s", "t=3s", "t=4s", "t=5s", "t=9s"},
PlotLabel ->
"Temperature distribution along line x=0.5L ar different times",
LabelStyle -> RGBColor[0, 0, 0]]


It's interesting to analyze velocity profiles in different cross sections

VelocProfArr = Table[
ParametricPlot[
{
{UxFun[t, 0, y], y},
{UxFun[t, l/2, y], y},
{Sin[t]*6 (y - y^2), y}
},
{y, 0, 1}, AspectRatio -> 0.25, Frame -> True,
FrameStyle -> RGBColor[0, 0, 0], FrameLabel -> {"Velocity", "y"},
BaseStyle -> 14, PlotRange -> {{-1.8, 1.8}, {0, 1}},
LabelStyle -> Black, PlotLabel -> "time=" <> ToString[t],
PlotLegends -> {"x=0", "x=0.5L", "Poiseuille profile"},
ImageSize -> 500]
, {t, ti, tf, 0.01*(tf - ti)}];

ListAnimate@VelocProfArr


$$V_x$$ distribution in longitudinal direction at $$y=d/2$$ is as follows

VelocLongArr = Table[
Plot[
{UxFun[t, x, 0.5],
Sin[t]*6*(0.5 - 0.5^2)
},
{x, 0, l}, AspectRatio -> 0.25, Frame -> True,
FrameStyle -> RGBColor[0, 0, 0], FrameLabel -> {"x", "Velocity"},
BaseStyle -> 14, PlotRange -> {{0, l}, {-1.8, 1.8}},
LabelStyle -> Black, PlotLabel -> "time=" <> ToString[t],
PlotLegends -> {"y=0", "Poiseuille flow"}, ImageSize -> 500]
, {t, ti, tf, 0.01*(tf - ti)}];

ListAnimate@VelocLongArr


As we can see, the velocity field differ from Poiseuille profile. Thereby the flow in channel can not be considered to be fully developed.

• @Avrana Did you install MeshTools` package prior to simulation? Feb 10, 2022 at 16:28
• Do not forget to install it. I gave the link on this package in answer Feb 10, 2022 at 16:42
• @AlexTrounev No it is not a type. As I already said I use dimensionless temperature $\tilde{T}=k_fT/(qd)$ in calculations Feb 11, 2022 at 10:02
• @AlexTrounev When such dimensionless variables are used the BC on the bottom wall reads as follows $-\frac{\partial \tilde{T}}{\partial y}\Big |_{y=-1}=k_f/k_s$ Feb 11, 2022 at 10:06
• @Avrana By the way. My first estimates whether the flow can be considered to be fully developed were wrong. Analysis of velocity field (animations in my post) suggest otherwise. Velocity profile differ from Poiseuille profile appreciably. Feb 11, 2022 at 18:37