5
$\begingroup$

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:

enter image description here

with the solid volume fraction:

enter image description here

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: enter image description here

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

$\endgroup$
1
  • 1
    $\begingroup$ What actually do you try to solve? $\endgroup$ Aug 30, 2021 at 3:32

2 Answers 2

4
$\begingroup$

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]

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
5
  • $\begingroup$ great! thanks a lot for your support! $\endgroup$
    – ABCDEMMM
    Nov 2, 2021 at 12:18
  • $\begingroup$ @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. $\endgroup$
    – user21
    Nov 2, 2021 at 12:38
  • $\begingroup$ 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? $\endgroup$
    – ABCDEMMM
    Dec 8, 2021 at 13:40
  • $\begingroup$ @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. $\endgroup$
    – user21
    Dec 8, 2021 at 14:29
  • $\begingroup$ thanks a lot for your support again! $\endgroup$
    – ABCDEMMM
    Dec 8, 2021 at 14:50
5
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ @user21this pde formulation can be found in: github.com/geo-fluid-dynamics/phaseflow-fenics/blob/master/… $\endgroup$
    – ABCDEMMM
    Oct 28, 2021 at 20:04
  • $\begingroup$ @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? $\endgroup$
    – user21
    Oct 29, 2021 at 5:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.