# Refining mesh size leads to absurd results for a coupled heat transfer FEM model

I have been recently solving a conjugate heat transfer problem, which involves fully-reversing or reciprocating flow of fluid over a heated block of solid. The problem is 2D and the temperature field is being evaluated for both the solid and the fluid. The following code (courtesy Alex) uses the finite-element approach:

Needs["NDSolveFEM"]
{f = 8;
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.069*10^3,
cps = 502.4}; u0 = 3.0; nu = mu/rho; om = 2 Pi f;
tflow = 0.125;
t0 = tflow/100;
NV = 2 f tflow;
nn = Round[NV \[Pi]/(om t0)]
Ti = 307.0; q = 5000.0/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[
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])) -
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

Tsm[x_] = nn^-1 Sum[Tfs[i][x, -e/2]*Ti - 273.16, {i, 1, nn}];
Plot[Tsm[x], {x, 0, L}, PlotRange -> Full, GridLines -> Automatic]

Tfm[y_] = (1/nn) Sum[Tfs[i][L/2, y]*Ti - 273.16, {i, 1, nn}];
Plot[{ Tfm[y]}, {y, -e, d}, PlotRange -> Full, GridLines -> Automatic]


Tsm[x] and Tfm[x] are the cyclic average temperature profile in the solid (along the line y=-e/2, i.e., streamwise direction) and the cyclic average temperature profile at a particular cross-section (x=L/2), respectively. The above code works properly for the given mesh settings in the above code i.e., MaxCellMeasure -> 10^-7. However, when I attempt a grid-independence test, i.e., run the simulation for a finer mesh like MaxCellMeasure -> 10^-8, the results become absurd. Following are the plots for Tsm[x] and Tfm[x] at the finer mesh settings:

As can be seen, the values are absurdly high. For the coarser mesh, we get reasonable plots, like the following:

I think, that with a finer mesh, the results should stay same or get more accurate. I have tried with a smaller time step, i.e., t0=tflow/500 with MaxCellMeasure -> 10^-8, but that did not help. How can this problem be resolved ?

For Bounty

The singularities have been removed by the accepted answer, by creating mesh in the non-dimensional space. However, problems arise when I conduct a grid-independence test, even using the modified code. The solutions Tsm[x] and Tfm[x] seem to converge to a particular solution up until MaxCellMeasure -> 10^-3. However, when the mesh size is reduced further using MaxCellMeasure -> 10^-4, they again diverge (although not absurdly).

These results have been summarized in excel sheets here. In these sheets Grid 1-4 refer to MaxCellMeasure -> 5*10^-3, 10^-3, 5*10^-4 and 10^-4, respectively. I do understand now that some numerical methods diverge when $$h \rightarrow 0$$, but converge until $$h \rightarrow h_{min}$$. My question is, if I would have started my calculation directly using MaxCellMeasure -> 10^-4, and then refined further, I would not be able to identify the optimum Mesh size where the solution is diverging. How can in such a scenario, one be assured that the solution is mesh independent ?

I tried accommodating the suggestions of @user21. I have used the following non-dimensionalisation scheme $$u=U/u_0, x=X/d, y=Y/d, p=\frac{P}{\rho u_0^{2}}, \tau = \omega t$$. $$\omega=2\pi f$$ is the flow reciprocation angular frequency, where $$f$$ is the frequency in Hz. Using this non-dimensional time allows us to march in time by $$\pi$$ intervals (at pi the b.c. are to be switched). With this scheme I have the following equations:

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

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

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

$$\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$$

The temperature is kept in the dimensional form.

In the non-dimensional space, the initial condition then becomes $$T(0,x,y)=0, u(0,x,y)=0, v(0,x,y)=0, p(0,x,y)=0$$, while the boundary condition becomes, $$u(\tau, 0, y)=\sin(t)$$ and $$-\frac{\partial T}{\partial y}=\frac{qd}{k_s}$$ at $$y=-e/d$$. The fluid thermal boundary conditions then become $$T(\tau, 0, y)=0, T(\tau+\pi, L/d, y)=0$$. Utilizing the contributions from Oleksii's answer (which was solved for a different form of non-dimensionalized equation of the same problem), I have utilized the monolithic approach of adding a momentum sink term to the momentum equations and modelling both the solid and fluid domains in one go. Following is the modified code:

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 = 2.0;(*flow oscillation frequency*)
period = 1/f;(*period*)
omega = 2*Pi/period;(*circular frequency*)
u0 = 1.;(*inflow velocity*)
q = 5000;(*heat flux density*)
Ti = 307;
re = d u0/(nu);
Pr = nu/AlphaF;
gamma = If[ElementMarker == 0, AlphaF/AlphaS, 1];
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}];

raster = {{{0, -e/d}, {l, -e/d}}, {{0, 0}, {l, 0}}};
MeshSolid = StructuredMesh[raster, {Nx, NyS}];

mesh = MergeMesh[MeshSolid, MeshFluid];
nodes = mesh["Coordinates"];

If[z < 0, 0, 1], {i, 1, Length[quads]}];
MeshTotal1 =
ToElementMesh["Coordinates" -> nodes,
MeshTotal2 = MeshOrderAlteration[MeshTotal1, 2];
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*)

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

c = If[ElementMarker == 0, 10^6,
0];(*define the constant in momentum sink term*)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/u0) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$u[t, x, y]$$\), {{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/u0) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$v[t, x, y]$$\), \!$$\*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[
y}] + (omega d^2) \!$$\*SubscriptBox[\(\[PartialD]$$, $${t}$$]$$T[t, 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[(d q)/(ks Ti), y == -e/d];

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 -> 8.37*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}]]

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

Plot[TFun[t][0.5 l, -e/(2 d)]*Ti, {t, 0, K*Pi},
GridLines -> Automatic, PlotRange -> Full]


The last line plots the solid temperature (dimensional) at a point for the simulated non-dimensional time.

However, as can be seen the temperature values are considerably blown up.

• What makes you believe that the messages you hit with Quiet is not the root of the issue? Commented Jan 5, 2023 at 7:00
• You never needs to specify Interpolation order for single dependent variable equations. Commented Jan 5, 2023 at 7:00
• Why do you use two different meshes? Commented Jan 5, 2023 at 7:02
• I would avoid using the iterative approach for the Navier-Stokes equation but solve the equation in one go. That will be more efficient. If you need to ramp up the inflow velocity, you should use ParametricNDSolve for that. With your code you have made some choices that miss an explanation in the text. To really answer your second question you'd need to get this FEM solver in a more standard form to perform experiments, I think. Commented Jan 5, 2023 at 7:07
• @user21 Thankyou for the comments. If I understand correctly, I should be solving the Navier-Stokes equation without the (u[x, y] - UX[i - 1][x, y])/t0 but instead use D[u[t,x,y],t] and leave the time-stepping to NDSolve by using the MaxStepSize command ? In that case, the iterator i should only control the switching of the boundary conditions at the end of each half-period to execute the flow reciprocation, right ? Also, the energy equation still needs to be solved separately or that too should be in conjunction with the momentum and continuity ? Commented Jan 5, 2023 at 10:07

This is a typical numerical instability due to FEM matrix inversion. On the fine mesh there are big numbers in Navier-Stokes equation of about $$1/h^2$$, therefore at $$h\rightarrow 0$$ we have singular matrix to be inverted. If we plot velocity

DensityPlot[UX[41][x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", PlotRange -> All,
AspectRatio -> Automatic, Frame -> False, PlotLegends -> Automatic]


we will see singular like point with value of about $$\pm 10^6$$. To avoid singularity we can increase $$h$$ by mapping reg1, reg2 to

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


Then option MaxCellMeasure -> 10^-8 with QuadElement number of about 12000 in mesh can be replaced by option MaxCellMeasure -> 10^-3 with QuadElement number of 13504. Hence we can make research with new numerical model as follows

Needs["NDSolveFEM"]
{f = 8;
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.069*10^3,
cps = 502.4}; u0 = 3.0; nu = mu/rho; om = 2 Pi f;
tflow = 0.125;
t0 = tflow/100;
NV = 2 f tflow;
nn = Round[NV \[Pi]/(om t0)]
Ti = 307.0; q = 5000.0/Ti;

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

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

re = d u0/nu; re1 = 1/re

UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P[0][x_, y_] := 0;

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] == Sin[om*i*t0], v[x, y] == 0},
x == L /d (1 - Sign[Sin[om*i*t0]])/2 && 0 < y < 1],
DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 0 || y == 1],
DirichletCondition[p[x, y] == 0,
x == L /d (1 + Sign[Sin[om*i*t0]])/2 && 0 < y < 1]};

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


Let check velocity in the central point of channel

ListPlot[Table[UX[i][L/d/2, 1/2], {i, 100}], PlotRange -> All]


It looks fine, therefore we can compute temperature as well

Tfs[0][x_, y_] := 307/Ti; appro =
With[{k = 2. 10^6}, ArcTan[k #]/Pi + 1/2 &];
ade[y_] := (ks + (kf - ks) appro[y])/(d u0)
rde[y_] := (cps rhos + (cp rho - cps rhos) appro[y]);

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/(u0 ), y == -e/d],
DirichletCondition[{T[x, y] == 1},
x == L/d (1 - Sign[Sin[om*i*t0]])/2 && 0 <= y <= 1]},
T, {x, y} \[Element] mesh1,
Method -> {"FiniteElement", "InterpolationOrder" -> {T -> 2}}] //
Quiet;, {i, 1, nn}],
ProgressIndicator[i, {1, nn}]] // AbsoluteTiming


Visualization

Tsm[x_] = nn^-1 Sum[Tfs[i][x, -e/d/2]*Ti - 273.16, {i, 1, nn}];
Plot[Tsm[x], {x, 0, L/d}, PlotRange -> Full, GridLines -> Automatic]

Tfm[y_] = (1/nn) Sum[Tfs[i][L/d/2, y]*Ti - 273.16, {i, 1, nn}];
Plot[{Tfm[y]}, {y, -e/d, 1}, PlotRange -> Full, GridLines -> Automatic]


• Really appreciate this detailed analysis. So you have now created the mesh in non-dimensional space, which allows for finer grid but with larger numbers, right ? In this scheme, the non-dimensionalisation is done as $u=u'/u0$, so the flux is $q/u0$ ? Commented Dec 26, 2022 at 6:41
• @Avrana Yes it is, in this case we use nondimensional equations for velocity with scale u0, d. Please note, that re is the Reynolds number, and 1/re= 0.00182958. Commented Dec 26, 2022 at 11:19
• I tried this new modified code and ran a mesh independence test. I started from MaxCellMeasure -> 5*10^-3 and went down till MaxCellMeasure -> 5*10^-4. The results of Tsm[x] and Tfm[y] converging down till this refinement. However, when I tried MaxCellMeasure -> 10^-4, the profiles showed a different behaviour (although not unphysical like the old code). Since, these calculations take a long time, I am attaching my results in the form of an excel sheet here. Commented Dec 27, 2022 at 3:13
• @Avrana This is FEM matrix inversion problem. Note, that problem coming from large value rde[y]. We can avoid this problem just using appropriate normalization, for example, ade[y_] := (ks + (kf - ks) appro[y])/(d u0 cps rhos ); rde[y_] := (1 + (cp rho /(cps rhos-1) appro[y]); and q/(u0 cps rhos). Commented Jan 6, 2023 at 17:07
• @Avrana It could be better to ask a new question in a new topic. Commented Jan 11, 2023 at 4:24