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

Question

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["NDSolve`FEM`"]
{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[
      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] == 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][
          ade[y]*Inactive[Grad][T[x, y], {x, y}], {x, y}] == 
        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["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 = 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"];
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];
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[
     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/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[
     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/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[
     Div][((ade[y]/rde[y]) Inactive[Grad][T[t, x, y], {x, y}]), {x, 
     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.

Answer 1

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["NDSolve`FEM`"]
{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][
          ade[y]*Inactive[Grad][T[x, y], {x, y}], {x, y}] == 
        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]