Modified Heat Transfer in Fluid Flow

Question

I am trying to simulate Modified Heat Transfer in Fluid Flow (based on Buoyancy-Driven Flow in a Square Cavity ).

The modified heat transfer takes the form:

with the solid volume fraction:

The whole model is defined as:

ClearAll["Global`*"]

Needs["NDSolve`FEM`"]
sizes = {length -> 1, hight -> 1};
\[CapitalOmega] = Rectangle[{0, 0}, {length, hight}] /. sizes;
em = ToElementMesh[\[CapitalOmega] , MaxCellMeasure -> .0001, 
   "MeshOrder" -> 2];

Pr = 50;
Ra = 2.27*10^5;
timetol = 0.00001;
Th = 1.0;
Tc = -0.01;


ClearAll[\[Nu]]
op = {
   
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] + 
    Inactive[Div][(- 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}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y], 
   
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] + 
    Inactive[Div][(-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}] + 
\!\(\*SuperscriptBox[\(p\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] - Ra/Pr*T[t, x, y], 
   
\!\(\*SuperscriptBox[\(u\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y] + 
\!\(\*SuperscriptBox[\(v\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, y],
   
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x, 
      y]*(1.0 + 1/0.045*20* Sech[40* (0.01 - T[t, x, y])]^2 ) + 
    Inactive[Div][(-Inactive[Grad][T[t, x, y], {x, y}]), {x, 
      y}] + {u[t, x, y], v[t, x, y]} . 
     Inactive[Grad][T[t, x, y], {x, y}]};


wall = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True];
reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];
temperatures = {DirichletCondition[T[t, x, y] == Th, x == 0], 
   DirichletCondition[T[t, x, y] == 0, x == length]};
bcs = {wall, reference, temperatures} /. sizes;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, 
   T[0, x, y] == 0};

Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = 
    NDSolveValue[{op == {0, 0, 0, 0}, bcs, ic}, {u, v, p, 
      T}, {x, y} \[Element] em, {t, 0, timetol}, 
     Method -> {"TimeIntegration" -> {"IDA", 
         "MaxDifferenceOrder" -> 2}, 
       "PDEDiscretization" -> {"MethodOfLines", 
         "SpatialDiscretization" -> {"FiniteElement", 
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, 
             T -> 2}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}])];], currentTime]

However, this modified simulation model produces some numerical errors:

How to set the solver parameters that we can fix the bug "Matrix SparseArray[<1395840>, {101404, 101404}] is singular ..."

Answer 1

Here is a version of the linked example, works just fine.

\[Phi][temp_, tr_, r_] := 1/2 (1 + Tanh[(tr - temp)/r])
regParam = 0.005;
Plot[\[Phi][T, 0, regParam], {T, -0.6, 0.6}, PlotRange -> All]

steNr = 0.045;
vars = {T[t, x], t, {x}};
pars = <|"DiffusionCoefficient" -> 1, 
   "HeatSource" -> steNr*\[Phi][T[t, x], 0, regParam]|>;
tHot = 1;
tCold = -0.01;
ics = T[0, x] == 
  With[{tHot = tHot, tCold = tCold, tMeltThickness = 0.1}, 
   If[x < tMeltThickness, tHot, tCold]
   ]

pde = {HeatTransferPDEComponent[vars, pars] == 0,
  HeatTemperatureCondition[x == 0, vars, 
   pars, <|"SurfaceTemperature" -> tHot|>], 
  HeatTemperatureCondition[x == 1, vars, 
   pars, <|"SurfaceTemperature" -> tCold|>]
  }

tEnd = 1;
fun = NDSolveValue[{pde, ics}, 
  T, {t, 0, tEnd}, {x} \[Element] Line[{{0}, {1}}], 
  Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}}]


Manipulate[
 Plot[fun[t, x], {x} \[Element] Line[{{0}, {1}}], 
  PlotRange -> {-0.1, 1.2}], {t, 0, tEnd}]

Manipulate[
 Plot[\[Phi][fun[t, x], 0, regParam], {x} \[Element] Line[{{0}, {1}}],
   PlotRange -> {-0.1, 1.2}], {t, 0, tEnd}]

Answer 2

This is not an answer but a tip how to change the implicit solver:

"TimeIntegration" -> {"IDA",(*"MaxDifferenceOrder"->2,*)
  "ImplicitSolver" -> {"Newton", "LinearSolveMethod" -> "Pardiso"}}

Sometimes this is useful when you get a message that a matrix is singular. In version 13 this will be the default. Unfortunately, this does not solve problem at hand; I suspect that the PDE model is not quite correct but I can not prove it.