1
$\begingroup$

This is a continuation of the post I've made Unable to solve Delay PDEs Error in Boussinesq Approximation. I apologise if I shouldn't have posted a seperate question for this but I think that the issue is different enough. I'm trying to solve the set of equations below describing the flow of a pot of water being heated slightly. The equations are 2D axisymmetric in nature. Using the suggestions from the linked post, I have updated my code as such:

x0 = 1;
y0 = 1;
\[CapitalOmega] = Rectangle[{0, 0}, {x0, y0}];
\[Mu]water = 8.9*10^-4;
\[Alpha]water = 2.1* 10^-4;
kwater = 0.6;
Cwater = 4184;
\[Rho]water = 1000;
g = 9.78;

eqn1=D[vr[r,z,t],t]+{vr[r,z,t],vz[r,z,t]}.Inactivate[Grad[vr[r,z,t],{r,z}],Grad]+D[p[r,z,t],r]/\[Rho]water-\[Mu]water(Inactivate[Div[Inactivate[Grad[vr[r, z, t], {r, z}], Grad], {r, z}], Div]+D[vr[r,z,t],r]/r-vr[r,z,t]/r^2)/\[Rho]water==0;
eqn2=D[vz[r,z,t],t]+{vr[r,z,t],vz[r,z,t]}.Inactivate[Grad[vz[r,z,t],{r,z}],Grad]+D[p[r,z,t],z]/\[Rho]water-\[Mu]water (Inactivate[Div[Inactivate[Grad[vz[r, z, t], {r, z}], Grad], {r, z}], Div]+D[vz[r,z,t],r]/r)/\[Rho]water  ==g \[Alpha]water T[r,z, t];
eqn3=D[vz[r,z,t],z]+D[r vr[r,z,t],r]/r==0;
eqn4=D[T[r,z,t],t]-kwater ((D[r*D[T[r,z,t],r],r]/r)+D[T[r,z,t],z,z])/(\[Rho]water Cwater)+{vr[r, z, t], vz[r, z, t]} . Inactive[Grad][T[r, z, t], {r, z}]==0;


wall = DirichletCondition[{vr[r, z, t] == 0, vz[r, z, t] == 0}, r == 1];
reference = DirichletCondition[p[r, z, t] == 0, r == 0 && z == 0];
temperatures = {DirichletCondition[T[r, z, t] == 1, r == 1]};
bcs = {wall, reference, temperatures};
ic = {vz[r, z, 0] == 0, vr[r, z, 0] == 0, p[r, z, 0] == 0, T[r, z, 0] == 0};
Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = NDSolveValue[{eqn1, eqn2, eqn3,eqn4, bcs, ic}, {vz, vr, p, T}, {r, z} \[Element] \[CapitalOmega], {t, 0, 200},Method -> {"PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement","MeshOptions" -> {"MaxCellMeasure" -> 0.0005}, "InterpolationOrder" -> {vz -> 2, vr -> 2, p -> 1, T -> 2}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];],currentTime]

Using EvaluationMonitor, I can see how much time has been computed by NDSolve. However, after running the code for a while I get that the time lingers at around 1.15s without any significant change. I'm not sure what the source of error is or how to fix the code. Could it be because only some parts of the equations are inactivated?

Improve this question
$\endgroup$
2
  • $\begingroup$ Unfortunately even in inactive form it runs up to t=0.0004 and then has a stop with the message 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 is turbulent convection. Maybe we need some advanced method to solve the problem. $\endgroup$
    – Alex Trounev
    Mar 19 at 6:30
  • $\begingroup$ To check is this is really turbulent, we'd need to estimate the Mach number u/c where u is the inflow velocity and c is the speed of sound for the fluid in question. Should the Mach number be larger then, say, 0.3 then, yes, we'd need to consider turbulence. $\endgroup$
    – user21
    Mar 20 at 7:51

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy