Message Text not Found in FEM Model

Question

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["NDSolve`FEM`"]
\[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!!

Answer 1

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]