1
$\begingroup$

I am trying to use NDSolve to find the numerical solution to a system of 4 first order PDEs with 8 boundary conditions (I think they are Dirichlet type conditions). There are 4 dependant variables and 2 independant variables. I am using cylindrical coordinates (Rr, Ttheta, Zz) but Ttheta is not used at all and so is excluded (may be the cuase of the problems?). I cannot for the life of me get Mathematica to run the code, I must be doing something wrong but I can't seem to figure out what. If anyone is able to help in anyway, that would be much appreciated.

The error that this code creates changes quite often depending on little 'fixes' I have attempted to make, however I have been slowly learning that a lot of the error messages often don't point to the cause of the problem so I have left them out of this post. The relevant equations in the code are also in the tradition form bellow the code just for ease of reading.

(*Constants*)
mi = 1.67*10^-27;
Vn = 10^4 ;
T = 5;
K = 8.617*10^-5;
Z = 1;

EQ0 = Vr[Rr, Zz]*mi*D[Log[ni[Rr, Zz]], Rr] + 
   Vz[Rr, Zz]*mi*D[Log[ni[Rr, Zz]], Zz] + mi*D[Vr[Rr, Zz], Rr] + 
   Vr[Rr, Zz]*mi + mi*D[Vz[Rr, Zz], Zz] == 0

EQ1 = -Vt[Rr, Zz]^2 == -Rr*((Z*K)/mi)*(D[T, Rr] + T*D[Log[ni[Rr, Zz]], Rr]) - 
   Rr*Vn*Vr[Rr, Zz]

EQ2 = Vr[Rr, Zz]*D[Vt[Rr, Zz], Rr]*Rr + Vr[Rr, Zz]*Vt[Rr, Zz] == -Vn*Vt[Rr, Zz]*Rr

EQ3 = Vr[Rr, Zz]*D[Vz[Rr, Zz], Rr] + Vz[Rr, Zz]*D[Vz[Rr, Zz], Zz] == 
    -((Z*K)/mi)*(D[T, Zz] + T*D[Log[ni[Rr, Zz]], Zz]) - Vn*Vz[Rr, Zz]

sol = NDSolve[{EQ0, EQ1, EQ2, EQ3,
    Vz[Rr, 0] == 10^3,
    Vt[Rr, 0] == 10^5*Rr,
    Vr[Rr, 0] == 0,
    Log[ni[Rr, 0]] == -(Rr^2/0.028^2) + Log[8.5*10^8],
    Vz[0.028, Zz] == 0,
    Vt[0, Zz] == 0,
    Vr[0, Zz] == 0,
    Log[ni[0.028, Zz]] == 0},
    {Vz, Vt, Vr, ni},
    {Rr, 0, 0.028}, {Zz, 0, 0.6}, AccuracyGoal -> 2, PrecisionGoal -> 2]

enter image description here Thanks!

$\endgroup$
  • $\begingroup$ The code contains derivatives of T, which is a constant. Is that intentional? $\endgroup$ – bbgodfrey Jan 22 '18 at 0:50
  • $\begingroup$ Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$ – bbgodfrey Jan 22 '18 at 0:57
  • $\begingroup$ @bbgodfrey yeah sorry I should have explained, that is intentional as T can be a function, however at the moment this is simplified a bit for the purposes of getting some sort of result and understanding how to properly use NDSolve. $\endgroup$ – Michael z Jan 22 '18 at 1:01
  • $\begingroup$ NDSolve appears to want to solve the PDE system as an initial value problem, which it is not. It would be nice to solve it using FEM, but NDSolve definitely cannot do so for nonlinear equations. Bottom line, I am not sure that NDSolve can handle this problem. (@MichaelE2, do you have an opinion on this?) $\endgroup$ – bbgodfrey Jan 22 '18 at 1:58
  • 1
    $\begingroup$ @bbgodfrey, even if FEM could handle non-linear PDEs, the FEM would not be a good fit here since this is a convection dominated PDE and one most likely would have to use some sort of stabilization technique. $\endgroup$ – user21 Jan 22 '18 at 7:19

Your Answer

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

Browse other questions tagged or ask your own question.