# Modified Heat Transfer in Fluid Flow

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["NDSolveFEM"]
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 ..."

• What actually do you try to solve? Commented Aug 30, 2021 at 3:32

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


• great! thanks a lot for your support! Commented Nov 2, 2021 at 12:18
• @ABCDEMMM, let me know if you can get a 2D version of this coupled to the Navier-Stokes equation. That would be a nice PDE model. Commented Nov 2, 2021 at 12:38
• in order to include the nonlinear term ("\phi=1/2 (1 + Tanh[(tr - temp)/r])"), must we use "HeatTransferPDEComponent" method in MMA 12.3? is it true in this case? Commented Dec 8, 2021 at 13:40
• @ABCDEMMM, I am not sure I understand the question. You can use HeatTransferPDEComponent but you certainly can also do without. HeatTransferPDEComponent is there to make the set up the equations easier, but if you set the equations up right, you do not need it. Commented Dec 8, 2021 at 14:29
• thanks a lot for your support again! Commented Dec 8, 2021 at 14:50

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.

• @user21this pde formulation can be found in: github.com/geo-fluid-dynamics/phaseflow-fenics/blob/master/… Commented Oct 28, 2021 at 20:04
• @ABCDEMMM I am on vacation now, but I'll have a look next week. The heat equation tutorial has a phase flow section. Maybe you can couple that to fluid flow? Commented Oct 29, 2021 at 5:11