I'm trying to solve for the heat transfer in 2D within a fluid. The geometry I'm using is a rectangle of water heated from the bottom. The code I'm using is as follows:

Needs["NDSolveFEM"]
\[Rho]water = 1000;
kwater = 0.598012;
\[Mu]water = 0.0010015961431205716;
Cwater = 4184;
\[Alpha]water = 0.00018;
g = 9.81;
Theat = 375;

HeatModel= {Cwater \[Rho]water D[T[x,y,t],t]+Cwater \[Rho]water u[x,y,t],Cwater \[Rho]water v[x,y,t]}. Grad[T[x, y, t],{x, y}]+Div[(-kwater Grad[T[x, y, t],{x, y}]),{x, y}];
CFDModel={\[Rho]water {u[x, y, t], v[x, y, t]} . Grad[u[x, y, t],{x, y}]+Div[(-\[Mu]water Grad[u[x, y, t],{x, y}]),{x, y}] +D[p[x,y,t],x]+D[u[x,y,t],t],\[Rho]water {u[x, y, t], v[x, y, t]} . Grad[v[x, y, t],{x, y}] + Div[(-\[Mu]water Grad[v[x, y, t],{x, y}]),{x,y}]+D[p[x,y,t],y]+D[v[x,y,t],t],D[u[x,y,t],x]+D[v[x,y,t],y]}
{Subscript[F, x], Subscript[F, y]} = {0, \[Alpha]water g \[Rho]water (T[x, y, t] - 300)};

r0 = 0.03;
h = 0.2;
\[CapitalOmega] = Rectangle[{0, 0}, {r0, h}];
mesh = ToElementMesh[\[CapitalOmega], {{0, r0}, {0, h}}, "MaxCellMeasure" -> {"Length" -> 20*10^-4}];

bcbase = {DirichletCondition[T[x, y, t] == Theat, y == 0]};
bcnoslip1 =DirichletCondition[{u[x, y, t] == 0, v[x, y, t] == 0}, y == 0 || y == h];
bcpressure = DirichletCondition[p[x, y, t] == 0, y == h];
bcnoslip2 =DirichletCondition[{u[x, y, t] == 0, v[x, y, t] == 0}, x == 0 || x == r0];
ics1 = {T[x, y, 0] == 300, (D[T[x, y, t] == 0, t]) /. t -> 0};
ic2 = {u[x, y, 0] == 0, v[x, y, 0] == 0, p[x, y, 0] == 0};
bcs = {bcbase, bcnoslip1, bcpressure, bcnoslip2, ics1, ic2};

pde = {CFDModel == {Subscript[F, x], Subscript[F, y], 0}, HeatModel == 0, bcs};
solveNS =NDSolveValue[pde, {u, v, p, T}, {x, y} \[Element] mesh, {t, 0, 10}, Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement","InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}]


When the code is ran, the only error message thrown is that -- Message Text not fonud --. Interesintgly NDSolve outputs a result but the time range in the output goes from t=0s to t=0s, while the time domain provided in NDSolve is from t=0s to t=10s. Does anyone know what's wrong? Thanks!!

• I can not reproduce this in V13.2. I get: LinearSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular. Mar 13, 2023 at 13:50
• With this parameters it should be turbulent convection. Mar 14, 2023 at 15:45

According to message LinearSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular this problem can't be solved with using nonlinear FEM. Nevertheless we can solve it with using linear FEM as described here. First, we transform system to dimensionless form as follows

\[Rho]water = 1000;
kwater = 0.598012;
\[Mu]water = 0.0010015961431205716;
Cwater = 4184;
\[Alpha]water = 0.00018;
g = 9.81;
Theat = 375; Pr0 = Cwater \[Mu]water/kwater (*7.007682559574845*);
r0 = 0.03; u0 = kwater/r0/\[Rho]water/Cwater(* velocity scale = 4.76427660930529*^-6*);
t0 = r0/u0 (*time scale = 6296.863608088132*);
R = 300 \[Alpha]water g r0/u0^2 (Theat/300 - 1)(* force scale = 1.75037*10^8*);


With this scaling the system has a form

tmax = 200; Pr0 = 7.0; Ra = 1.7503709383883822*^8/Pr0; R =
Ra*Pr0(*1.7503709383883822*^8*); t0 = 6296.863608088132; dt =
10^-5; a = 0; h = 20/3;
\[CapitalOmega] = ImplicitRegion[0 <= x <= 1 && 0 <= y <= h, {x, y}];
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
\[CapitalRho][0][x_, y_] := 0;
TK[0][x_, y_] := 0; TX[0][x_, y_] := 0;

Do[{UX[i], VY[i], \[CapitalRho][i], TK[i], TX[i]} =
NDSolveValue[{{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])/dt - R*T[x, y]*Sin[a],
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] - R*T[x, y]*Cos[a] + (v[x, y] - VY[i - 1][x, y])/dt, D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /.
\[Mu] -> Pr0, {DirichletCondition[{p[x, y] == 0}, y == h && x == 1], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, y == 0 || y == h]},
{UX[i - 1][x, y]*D[T[x, y], x] + VY[i - 1][x, y]*D[T[x, y], y] + (T[x, y] - TK[i - 1][x, y])/dt -
(D[T[x, y], x, x] + D[T[x, y], y, y]) == NeumannValue[0, x == 0 || x == 1 || y == h], DirichletCondition[{u[x, y] == 0, v[x, y] == 0},
x == 0], DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 1], DirichletCondition[T[x, y] == 1, y == 0]}},
{u, v, p, T, Derivative[1, 0][T]}, Element[{x, y}, \[CapitalOmega]],
Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}, "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}],
{i, 1, tmax}];


Note, that tmax dt t0 is about 12.5937s while in the question t<=10s. Visualization

mesh = UX[tmax]["ElementMesh"];

Table[StreamDensityPlot[{UX[i][x, y],
VY[i][x, y]}, {x, y} \[Element] \[CapitalOmega],
StreamPoints -> Fine, PlotLegends -> Automatic,
ColorFunction -> "TemperatureMap",
PlotLabel -> Row[{"t =", dt i t0}],
StreamColorFunction -> GrayLevel], {i, 50, 200, 50}]


Table[DensityPlot[TK[i][x, y], {x, y} \[Element] mesh,
PlotPoints -> 50, ColorFunction -> "TemperatureMap",
PlotRange -> All, AspectRatio -> Automatic, Frame -> False,
PlotLabel -> i dt t0], {i, 30, 200, 10}]


This is first stage of turbulent convection. We can compute up to tmax=400 to see this turbulent move, that we can animate as

lst = Table[
DensityPlot[TK[i][x, y], {x, y} \[Element] mesh, PlotPoints -> 50,
ColorFunction -> "TemperatureMap", PlotRange -> All,
AspectRatio -> Automatic, Frame -> False], {i, 30, 400, 2}];

ListAnimate[lst]


• Oh wow, thanks alot for this!! Mar 16, 2023 at 13:54
• @Lucas Yes, you are right, it looks very realistic. Mar 16, 2023 at 14:58