I am trying to solve the convection-coupled melting benchmark in MMA 12.3 as presented in FEniCS/03-ConvectionCoupledMelting-MixedElement-AMR.ipynb. This model in MMA is based on the solution from user21 and energy-transport. In the simulation, "Matrix SparseArray is singular" problem cannot be solved. The final solution may look like this:
(see reference 1)
This notebook shows
Needs["NDSolve`FEM`"]
sizes = {length -> 1, hight -> 1};
\[CapitalOmega] = Rectangle[{0, 0}, {length, hight}] /. sizes;
\[Phi][temp_, tr_, r_] := 1/2 (1 + Tanh[(tr - temp)/r])
Plot[\[Phi][T, 0.01, 0.025], {T, -0.6, 0.6}, PlotRange -> All]
tvars = {T[t, x, y], t, {x, y}};
regParam = 0.025;
steNr = 0.045;
Tnr = 0.01;
pars = <|"MassDensity" -> 1, "SpecificHeatCapacity" -> 1,
"ThermalConductivity" -> 1/Pr*IdentityMatrix[2],
"HeatConvectionVelocity" -> {u[t, x, y], v[t, x, y]},
"HeatSource" -> (1/steNr)*\[Phi][T[t, x, y], Tnr , regParam]|>;
TransientHeatModel = HeatTransferPDEComponent[tvars, pars];
muL = 1.0;
muS = 10^8;
\[Mu][muL_, muS_, phi_] := muL + (muS - muL)*phi;
CFDModel = {
\!\(\*SuperscriptBox[\(u\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x,
y] + \[Mu][muL, muS, \[Phi][T[t, x, y], Tnr, regParam]]*
Inactive[Div][(-\[Nu] 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] + \[Mu][muL, muS, \[Phi][T[t, x, y], Tnr, regParam]]*
Inactive[Div][(-\[Nu] 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]};
parameters = {\[Nu] -> 1, Pr -> 56.2, Ra -> 327000.0};
tHot = 1.0;
tCold = -0.01;
ics = T[0, x, y] ==
With[{tHot = tHot, tCold = tCold, tMeltThickness = 0.1},
If[x < tMeltThickness, tHot, tCold]];
icuvp = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
Subscript[\[CapitalGamma],
temp] = {HeatTemperatureCondition[x == 0, tvars,
pars, <|"SurfaceTemperature" -> tHot|>],
HeatTemperatureCondition[x == 1, tvars,
pars, <|"SurfaceTemperature" -> tCold|>]};
Subscript[\[CapitalGamma], wall] =
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True];
Subscript[\[CapitalGamma], p] =
DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];
bcs = {Subscript[\[CapitalGamma], wall], Subscript[\[CapitalGamma],
p], Subscript[\[CapitalGamma], temp]};
pde =
{CFDModel == {0, 0, 0}, TransientHeatModel == 0, bcs, ics,
icuvp} /. parameters;
tEnd = 0.001;
{uVel, vVel, pressure, Tfun} =
NDSolveValue[
pde, {u, v, p, T}, {t, 0,
tEnd}, {x, y} \[Element] \[CapitalOmega], Method -> {
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}];
Your help would be much appreciated!
Remark:
1.) do we have time increments problem in this model? To debug the code, tEnd = 0.001 is already chosen in the original model (tEnd means the total time).
2.) energy-transport model has linear parameters in CFD model, in our case we have nonlinear term (\mu) in CFD model.
3.) (Very) fine mesh is used in this benchmark test, since the adaptive mesh is applied (Fenics solver).
[
] (see reference 1)
Updates:
1.) MeshOptions
If we test the model from AlexTrounev, we can find the following error in case we set: "MeshOptions" -> {"MaxCellMeasure" -> 0.001}
If we set "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}
Thus, a numerical issue still remains in this post, namely, why we cannot use very fine mesh (e.g. {"MaxCellMeasure" -> 0.0001}) for modeling convection-coupled melting in Mathematica 12.3.
Reference:
Zimmerman, Alexander & Kowalski, Julia. (2018). Monolithic Simulation of Convection-Coupled Phase-Change: Verification and Reproducibility: Proceedings of the 4th International Conference on Computational Engineering (ICCE 2017) in Darmstadt. 10.1007/978-3-319-93891-2_11.
