# Implementation of FEM to steady-state coupled fluid flow and heat transfer

This is a steady-state conjugate heat transfer problem (the time-independent version of this problem). The problem is conjugate as the energy equation is being solved in thermally connected solid and fluid domain. It describes 2D-flow of fluid over a heated solid block. I presumed it would be straightforward with the NDSolve fem option, but alas I ran into errors:

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.05;(*inflow velocity*)
q = 5000;(*heat flux density*)

uavg = N[\!$$\*SubsuperscriptBox[\(∫$$, $$0$$, FractionBox[$$π$$, $$omega$$]]$$u0\ Sin[omega\ t] \[DifferentialD]t$$\)/(π/omega)];
re = d uavg/(nu);

(*Meshing*)
Nx = 70;(*number of elements in x-direction*)
NyF = 50;(*number of elements in y-direction in fluid*)
NyS = 40;(*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]}];
MeshTotal2 = MeshOrderAlteration[MeshTotal1, 2];

(*Fluid inflow profile*)
Clear[TopWall, BottomWall, reference, HeatInpBC, op, c, rampFunction,
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*)

(*Material property switching function*)
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]);

(*Stationary momentum, continuity and energy operator*)
c = If[ElementMarker == 0, 10^6,
0]; op = {{{u[x, y], v[x, y]}}.Inactive[Grad][u[x, y], {x, y}] +
Inactive[
u[x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(∂$$, $${x}$$]$$p[x, y]$$\) +
c u[x, y], {{u[x, y], v[x, y]}}.Inactive[Grad][v[x, y], {x, y}] +
Inactive[
v[x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(∂$$, $${y}$$]$$p[x, y]$$\) + c v[x, y], \!$$\*SubscriptBox[\(∂$$, $${x}$$]$$u[x, y]$$\) + \!$$\*SubscriptBox[\(∂$$, $${y}$$]$$v[x, y]$$\), {u[x, y],
v[x, y]}.Inactive[Grad][(u0 d) T[x, y], {x, y}] -
rde[y], {x, y}]};

(*Boundary conditions*)
TopWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 1];
BottomWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y <= 0];
HeatInpBC = NeumannValue[(ks)/(cs rhos), y == -(e/d)];
reference = DirichletCondition[p[x, y] == 0., x == l && y == 0];Clear[u, v, p, t, HeatDBC];
Clear[HeatDBC, Inflow, Outflow, bcs, UxFun, UyFun, pressure, TFun];
Inflow = DirichletCondition[{u[x, y] == 1, v[x, y] == 0},
x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{p[x, y] == 0, v[x, y] == 0},
x == l && y > 0 && y < 1];
HeatDBC = DirichletCondition[T[x, y] == 0, x == 0 && y >= 0 && y <= 1];
bcs = {TopWall, BottomWall, Inflow, Outflow, reference, HeatDBC};

(*Solution*)
{UxFun, UyFun, pressure, TFun} =
NDSolveValue[{op == {0, 0, 0, HeatInpBC}, bcs}, {u, v, p,
T}, {x, y} \[Element] MeshTotal2,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}},
InitialSeeding -> {u[x, y] == 0, v[x, y] == 0, T[x,y]==0}];


On solving this leads to the following error:

FindRoot::nosol: Linear equation encountered that has no solution. FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

As can be seen in code, I have also tried to provide an initial seed for the system InitialSeeding -> {u[x, y] == 1, v[x, y] == 0, T[x,y]==0}, but to no avail. Also, I have kept the fluid velocity low u0=0.05 (uavg is lower actually), which should be rendering the non-linearity to be low. What should be changed to find a stationary solution to this problem ?

Update

If the energy equation is non-dimensionalised as $$T=\frac{T^*}{\frac{qd}{k_s}}, u =\frac{u^*}{u_{avg}}, v =\frac{v^*}{u_{avg}}$$, it becomes:

$$u_{avg} d \rho c_p \bigg(u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y}\bigg) -k \nabla^2 T=0$$ with the boundary condition $$-\frac{\partial T}{\partial y}=1$$ at $$y=-e/d$$.

Then for MMA FEM, the boundary conditon should become NeumannValue[ks, y==-e/d]

After some recommendations in the comments, I have updated my approach. Please copy code till the definition of re from the already posted code.

Ti = 307.;
reg1 = ImplicitRegion[0 <= x <= L/d && 0 <= y <= 1, {x, y}];
reg2 = ImplicitRegion[0 <= x <= L/d && -e <= y <= 1, {x, y}];

mesh = ToElementMesh[FullRegion[2], {{0, L/d}, {0, 1}},
MaxCellMeasure -> 5 10^-4];
mesh1 = ToElementMesh[FullRegion[2], {{0, L/d}, {-e/d, 1}},
MaxCellMeasure -> 5 10^-4];

op = {{{u[x, y], v[x, y]}}.Inactive[Grad][u[x, y], {x, y}] +
Inactive[
u[x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(∂$$, $${x}$$]$$p[x, y]$$\), {{u[x, y], v[x, y]}}.Inactive[Grad][v[x, y], {x, y}] +
Inactive[
v[x, y], {x, y}]), {x, y}] + \!$$\*SubscriptBox[\(∂$$, $${y}$$]$$p[x, y]$$\), \!$$\*SubscriptBox[\(∂$$, $${x}$$]$$u[x, y]$$\) + \!$$\*SubscriptBox[\(∂$$, $${y}$$]$$v[x, y]$$\)};
TopWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 1];
BottomWall = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y <= 0];
reference = DirichletCondition[p[x, y] == 0., x == l && y == 0];

Clear[u, v, p, t];
Inflow = DirichletCondition[{u[x, y] == 1, v[x, y] == 0},
x == 0 && y > 0 && y < 1];
Outflow =
DirichletCondition[{p[x, y] == 0, v[x, y] == 0},
x == l && y > 0 && y < 1];
bcs = {TopWall, BottomWall, reference, Inflow, Outflow};

{UxFun, UyFun, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} ∈
mesh, Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}];
DensityPlot[UxFun[x, y], {x, 0, l}, {y, 0, 1},
ColorFunction -> "Rainbow", AspectRatio -> Automatic]
appro = With[{k = 2. 10^6}, ArcTan[#1 k]/\[Pi] + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])
rde[y_] := (cs rhos + (cf rhof - cs rhos) appro[y]);

ux = If[y <= 0, 0, UxFun[x, y]];
vy = If[y <= 0, 0, UyFun[x, y]];

T = NDSolveValue[{ uavg d rde[y] (ux \!$$\*SubscriptBox[\(\[PartialD]$$, $${x}$$]$$T[x, y]$$\) + vy \!$$\*SubscriptBox[\(\[PartialD]$$, $${y}$$]$$T[x, y]$$\)) -
y}] == NeumannValue[ks , y == (-e/d)],
DirichletCondition[{T[x, y] == 0}, x == 0 && 0 <= y <= 1]},
T, {x, y} \[Element] mesh1,
Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2}}];

DensityPlot[T[x, y], {x, 0, l}, {y, -e/d, 1},
AspectRatio -> Automatic, ColorFunction -> "Rainbow"]

DensityPlot[T[x, y], {x, 0, l}, {y, -e/d, 1},
AspectRatio -> Automatic, ColorFunction -> "Rainbow"]

Plot[{Ti + T[x, -e/d/2]*(q d/ks), Ti + T[x, 1/4]*(q d/ks)}, {x, 0, l}]


The solid (at $$y=\frac{-e/d}{2}$$) and fluid (at $$y=\frac{1}{4}$$) temperature profiles are:

Some issues:

1. While solving the energy equation, the following warnings are delivered:

Are these serious?

1. If I want to calculate the average solid temperature (line-integrated along its height, i.e., $$y-$$ direction), what should be the way? Basically, I want to calculate:

$$T_{s,avg}(x) = \frac{1}{e/d}\int_{-e/d}^0 T(x,y)\mathrm{d}y$$

Currently I tried using:

T1[x_] = (1/(e/d)) NIntegrate[Ti + T[x, y]*(q d/ks), {y, -e/d, 0}]

1. Let us assume from experiments, I know the length, width of a block to be $$L,W$$ (in m) to which some $$Q$$ (in W) heat is supplied. Now imagine, I am modelling this system in 2D (in Mathematica), where my model only has the length ($$L$$) feature and not the width ($$W$$). In this scenario, I have to emulate the heat input as a flux ($$q′′$$, q in this question) condition (Neumann value). How should I calculate this flux, to replicate the experiment ? Should it be $$q = \frac{experimetal \space Q}{L∗1}$$ or $$q = \frac{experimental \space Q}{L∗W}$$? Similar argument also follows for deciding the uavg value for the simulation as volumetric flow rate is known from experiments. Should it be $$\frac{experimental \space flowate}{d*1}$$ or $$\frac{experimental \space flowate}{d*W}$$ for the simulation.

Bounty

Let us look at the energy balance scenario. Let us assume that the width of the block is 1 (in $$z-$$ direction). In that case, analytical considerations say that the fluid outlet temperature can be calculated using Solve[q L 1 == rhof uavg cf d 1 (x - Ti), x], i.e., $$q*L*1 = \rho (d*1) u_{avg} c_p (T_{outlet} - T_i)$$, which results in

{{x -> 363.828}}


However, using this MMA model (in the Update section), we get

T2[x] is the average fluid temperature along the outlet boundary

T2[x_?NumericQ] := (1/1) NIntegrate[T[x, y], {y, 0, 1}, MaxRecursion -> 12];
Ti + Evaluate[T2[l]]*(q d/ks)

323.395


So there seems to be huge energy imbalance. Interstingly, I have seen that as the flow velocity, i.e., u0 is decreased, the amount of this discrepancy grows even further (hinting at diffusion being the main culprit). With increasing flow velcoity, this imbalance reduces.

• @Avrana Interpolating function warnings coming when we use mesh and exact coordinates together like in your code. Use T1[x_?NumericQ] := (1/(e/d)) NIntegrate[ Ti + T[x, y]*(q d/ks), {y, -e/d, 0}]; Plot[Evaluate[T1[x]], {x, 0, L/d}] Jan 15 at 7:54
• See this for verification. You should be using HeatTransferPDEComponent Jan 15 at 12:02
• The Interpolation issue is serious you should take care of this. Look at the message ref page for more information Jan 15 at 12:04
• For boundary integrals have a look at this. Jan 15 at 12:06
• @Avrana In a case of $q=q(x,y,z)$ with $0\le z \le W$ we can use, for instance, mean value $qm(x,y)=\int_0^W q(x,y,z)dz/W$. Jan 15 at 16:03

We can use this code for the case when inflow velocity is constant in time. The flow and temperature field in this system will to come to equilibrium after some period. The non-dimentional governing equations are solved here. Dimensional temperature and fluxes are than recalculated.

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 = 100;(*number of elements in x-direction *)
NyF = 20;(*number of elements in y-direction in fluid*)
NyS = 10;(*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*)
MeshTotal = MeshOrderAlteration[MergeMesh[MeshSolid, MeshFluid], 2];


Solution of PDE:

Clear[TopWall, BottomWall, reference, HeatInpBC, op, rampFunction, sf,
UinfProfile, Profile, x, y, t];
rampFunction[min_, max_, c_, r_] :=
Function[t,
If[t < 2 c, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t]),
1]]
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];
reference = DirichletCondition[p[t, x, y] == 0., x == 0 && y == 0];
HeatInpBC = NeumannValue[1, y == -e/d];
Inflow =
DirichletCondition[{u[t, x, y] == sf[t*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]*UinfProfile[y],
v[t, x, y] == 0}, x == l && y > 0 && y < 1];
HeatDBC =
DirichletCondition[T[t, x, y] == 0, x == 0 && y > 0 && y <= 1];

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

ti = 0;
tf = 10000 u0/d;
Clear[UxFun, UyFun, pressure, TFun];

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

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];


Visualization

Temperature hystory in point $$\{0.5L,0.5d\}$$ looks as follows:

pic1 = Plot[Ti + (d*q)/ks TFun[t*u0/d, l/2, 0.5], {t, 0, tf*d/u0},
PlotStyle -> {Thickness[0.003], RGBColor[0, 0, 0]},
PlotRange -> All, Frame -> True,
FrameLabel -> {"time [s]", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 400,
LabelStyle -> RGBColor[0, 0, 0]]


Temperature distributions along channel center line and in central cross section at equilibrium are as follows:

pic2 = Plot[Ti + (d*q)/ks TFun[tf, x/d, 0.5], {x, 0, L},
PlotStyle -> {Thickness[0.003], RGBColor[0, 0, 0]},
PlotRange -> All, Frame -> True,
FrameLabel -> {"x [m]", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 300,
LabelStyle -> RGBColor[0, 0, 0]]

pic3 = Plot[Ti + (d*q)/ks TFun[tf, 0.5, y/d], {y, -e, d},
PlotStyle -> {Thickness[0.003], RGBColor[0, 0, 0]},
PlotRange -> All, Frame -> True,
FrameLabel -> {"y [m]", "Temperature"},
FrameStyle -> RGBColor[0, 0, 0], BaseStyle -> 14, ImageSize -> 300,
LabelStyle -> RGBColor[0, 0, 0]]


Integral heat balance check

Let's write the integral heat balance for stationary sate:

$$$$c_f\rho_f\int_0^d\left(T(L,y)u_x(L,y)-T(0,y)u_x(0,y) \right)dy=\int_{\partial \Omega}k\frac{\partial T}{\partial n}ds$$$$

@user21 gave useful link where similar surface integrals are calculated. Let's take advantage of this approach.

bmesh = ToBoundaryMesh[MeshTotal];
tolerance = 10^-2 hy;
leftSurfFluid[x_, y_] := If[Abs[x] <= tolerance && y >= 0, 1, 0]
rightSurfFluid[x_, y_] := If[Abs[x - l] <= tolerance && y >= 0, 1, 0]
leftSurfSolid[x_, y_] := If[Abs[x] <= tolerance && y < 0, 1, 0]
rightSurfSolid[x_, y_] := If[Abs[x - l] <= tolerance && y < 0, 1, 0]
bottomSurf[x_, y_] := If[Abs[y + e/d] <= tolerance && y < 0, 1, 0]
topSurf[x_, y_] := If[Abs[y - 1] <= tolerance && y < 0, 1, 0]



Calculation of terms in right hand side of balance equation

Q1 = -q*d*(kf/ks)*
NIntegrate[
leftSurfFluid[x, y]*GradX[tf, x, y], {x, y} \[Element] bmesh]
Q2 = q*d*kf/
ks NIntegrate[
rightSurfFluid[x, y]*GradX[tf, x, y], {x, y} \[Element] bmesh]
Q3 = -q*d*
NIntegrate[
leftSurfSolid[x, y]*GradX[tf, x, y], {x, y} \[Element] bmesh]
Q4 = q*d*
NIntegrate[
rightSurfSolid[x, y]*GradX[tf, x, y], {x, y} \[Element] bmesh]
Q5 = -q*d*
NIntegrate[bottomSurf[x, y]*GradY[tf, x, y], {x, y} \[Element] bmesh]
Q6 = q*d*
NIntegrate[topSurf[x, y]*GradY[tf, x, y], {x, y} \[Element] bmesh]


gives

-2.49716

-0.0133357

-5.72554

-0.0324436

200.

0.


Right hand side equals to 191.732 W/m

Q1 + Q2 + Q3 + Q4 + Q5 + Q6


Calculation of left hand side gives 192.601 W/m:

deltaQ = (NIntegrate[
rightSurfFluid[x, y]*UxFun[tf, x, y]*
TFun[tf, x, y], {x, y} \[Element] bmesh] -
NIntegrate[
leftSurfFluid[x, y]*UxFun[tf, x, y]*
TFun[tf, x, y], {x, y} \[Element] bmesh])*(d^2*q*u0/ks)*rhof*cf


The exact value of supplied heat flux at $$y=-e$$ equals to $$q\cdot L$$=200 W/m. As we can see the heat balance discrepancy is not so large. In addition it depends on FE mesh used.

Update

By using previous FE mesh we obtained numerically heat flux Q3=$$\int_0^dk_s\frac{\partial T}{\partial x}|_{x=0}dy=-5.72$$ W/m which is not equals to 0. It's confuses a little bit. However in the vicinity of point {0,0} the large temperature gradients appear and more fine FE mesh is required in this region for better approximation. Let's try the refined mesh:

NyS = 10; (*number of elements in solid*)
NyF = 20; (*number of elements in fluid*)
NxS = 100; (*number of elements in x-direction*)
hyMin = 1/40; (*linear dimension of the smallest element*)
hxMin = hyMin;
hy = 1./NyF;(*linear dimension of element in fluid*)

meshX1d =
Line[{{0}, {l}}], <| "Alignment" -> "Left",
"ElementCount" -> NxS,  "MinimalDistance" -> hyMin|>];
meshY1d = MergeMesh[
Line[{{-e/d}, {0}}], <| "Alignment" -> "Right",
"ElementCount" -> NyS,  "MinimalDistance" -> hyMin|>],
Line[{{0}, {1}}], <| "Alignment" -> "Left",
"ElementCount" -> NyF,  "MinimalDistance" -> hyMin|>]
];
MeshTotal = ElementMeshRegionProduct[meshX1d, meshY1d];
MeshTotal["Wireframe"]



In this case we get Q3=-2.03978 W/m, deltaQ=193.113 W/m.

• Thanks for the detailed answer again, especially for showing the application of heat balance check. I was having a hard time understanding it from the tutorials. Jan 19 at 5:50
• Everything works fine, for the tests I have run. The net energy balance works good. However, there is one discrepancy which I would like to point out, i.e., the value of Q3, which comes out as significant. In your example, it is around -17.698, although, this left solid boundary is insulated. However, the net energy balance is within 1.27 %. Insulated is the default b.c. in NDSolve. Jan 20 at 6:21
• As to value of Q3. Yes, at first glance this value confuses a little. However we should don't forget that FEM is the approximate method. And the heat fluxes recalculated from temperature field will tend to Neumann b.c. used with decreasing of element size near boundary. In addition temperature gradients near point {0,0} are large and adaptive mesh should be used in this region for better approximation. Jan 21 at 8:34
• As to experimental values of heat flux $Q$ and flow rate $\frac{dV}{dt}$. For inflow velocity we have $u_0=\frac{dV}{dt}\frac{1}{Wd}$ and for heat flux density supplied $q=\frac{Q}{dW}$. Of course the last formula is correct if heat flux is distributed uniformly on surface $y=-e$ Jan 21 at 8:53
• In addition, the flow can be considered to be plane (2D) when $W>>d$. Otherwise 3D problem should be solved. And it makes the problem more complicated, and computations become more time consuming Jan 21 at 9:03