1
$\begingroup$

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.

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]+vz[r,z,t]*D[vr[r,z,t],z]+vr[r,z,t]*D[vr[r,z,t],r]+D[p[r,z,t],r]/\[Rho]water-\[Mu]water (D[vr[r,z,t],z,z]+(D[r*D[vr[r,z,t],r],r]/r)-vr[r,z,t]/r^2)/\[Rho]water==0;
eqn2=D[vz[r,z,t],t]+vz[r,z,t]*D[vz[r,z,t],z]+vr[r,z,t]*D[vz[r,z,t],r]+D[p[r,z,t],z]/\[Rho]water-\[Mu]water (D[vz[r,z,t],z,z]+(D[r*D[vz[r,z,t],r],r]/r))/\[Rho]water==g \[Alpha]water T[x, y, 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)+vz[r,z,t] D[T[r,z,t],z]+vr[r,z,t] D[T[r,z,t],r]==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.0000625}, "InterpolationOrder" -> {vz -> 2, vr -> 2, p -> 1, T -> 2}}}}, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];],currentTime]

After running the code, I get the error that Delay PDEs are not supported by mathematica. When I search online, most of the methods to fix this error is to use the method of lines which I have used in NDSolve, but still have not been able to fix this. Anyone knows what's wrong?

Improve this question
$\endgroup$
10
  • $\begingroup$ @user21 I think this is another parsing bug? $\endgroup$
    – xzczd
    Mar 17 at 2:51
  • $\begingroup$ @xzczd I think the firs thing to do is to put it in inactive form. $\endgroup$
    – user21
    Mar 17 at 7:02
  • $\begingroup$ @user21 is it possible to put it into inactive form? I don’t know if I can specify grad and div in terms of cylindrical coordinates while using the axisymmetric condition. Also out of curiosity what’s the difference between the inactive and active forms when it comes to solving? $\endgroup$
    – Lucas
    Mar 17 at 7:08
  • $\begingroup$ Have a look at the DiffusionPDETerm (et al) ref pages, that have a "RegionSymmetry" -> "Axisymmetric" parameter. I have not tried this myself, you might need to fiddle a little with the equation too. Concerning the the inactive form, this specifies what the meaning of NeumannValue is, but also a few other things.... $\endgroup$
    – user21
    Mar 17 at 7:46
  • $\begingroup$ ...See here and here, Not quite what you are looking for but perhaps also useful is the axisymmetric Tubular Reactor example. $\endgroup$
    – user21
    Mar 17 at 7:46

0

Reset to default

Your Answer

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