2
$\begingroup$

I am trying to solve a system of the partial differential equation with the help of NDSolve and method of lines.

Mathematica code for the above-described problem is

Clear[B, lm, Pr, pde1, pde2, pde3]
pde1 = D[u[t, x, r], x]+ D[v[t, x, r], r] + v[t, x, r]/r == 0;
pde2 = D[u[t, x, r], t]+u[t, x, r]*D[u[t, x, r], x]+v[t, x, r]*D[u[t, x, r], r]==T[t, x, r]+(1/(1+lm))*(D[u[t, x, r], r,r]+(1/r)*D[u[t, x, r], r]+B*(D[u[t, x, r], r,r,t]+D[v[t, x, r], r]*D[u[t, x, r], r,r]+v[t, x, r]*D[u[t, x, r], r,r,r]+D[u[t, x, r], r]*D[u[t, x, r], x,r]+u[t, x, r]*D[u[t, x, r], x,r,r]+(1/r)*D[u[t, x, r], r,t]+v[t, x, r]*D[u[t, x, r], r,r]/r+u[t, x, r]*D[u[t, x, r], x,r]/r));
pde3 = D[T[t, x, r], t]+u[t, x, r]*D[T[t, x, r], x]+v[t, x, r]*D[T[t, x, r], r]==(1/Pr)*(D[T[t, x, r], r,r]+D[T[t, x, r], r]/r);

(*with initial and boundary conditions:*)

ics = {u[0, x, r] == 0, v[0, x, r] == 0, T[0, x, r] == 0};
With[{lb=10}, bcs = {{u[t, 0, r] == 0, v[t, 0, r] == 0, T[t, 0, r] == 0}, {u[t, x, 1] == 0, v[t, x, 1] == 0, T[t, x, 1] == 1}, {u[t, x, lb] == 0, (D[u[t,x,r], r] /. r -> lb) == 0,T[t, x, lb] == 0}}];

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n,"MinPoints" -> n, "DifferenceOrder" -> o}} mol[tf : False | True, sf_: Automatic] := {"MethodOfLines","Differentiate\[Beta]oundaryConditions" -> {tf,"ScaleFactor" -> sf}}

Clear@solfunc With[{pts = 20, lb = 20}, solfunc[B_, Pr_,lm_, tend_] := NDSolveValue[{pde1, pde2[B,lm], pde3[Pr], bcs,ics}, {u, v, T}, {t, 0, tend}, {x, 0, 1}, {r, 0, lb}, Method -> Union[mol[pts, 4], mol[True, 100]]]] (sollst[#] = solfunc[#, 0.7, 0.5, 10]) & /@ {0.1, 0.2}; // Quiet

Plot[{sollst[#][[1]][1, 1, r] & /@ {0.1} // Evaluate, sollst[#][[1]][1, 1, r] & /@ {0.2} // Evaluate}, {r, 0, 10},  PlotRange -> All, AxesOrigin -> {0, 0}, Frame -> True]

But it doesn't work...

$\endgroup$
6
  • 2
    $\begingroup$ "But it doesn't work..." is really a lazy and useless problem description. This fits really well to the usage of Quiet. Please, start to read error messages. That will help you to figure what went wrong. Also sprinkle a couple of semicoli, in particular in the line starting with Clear. $\endgroup$ Mar 10 '20 at 7:09
  • 1
    $\begingroup$ You seem to use pde2 and pde2 as functions, but they are not defined as such. Something that you could have inferred from reading the error messages... $\endgroup$ Mar 10 '20 at 7:13
  • $\begingroup$ Side note: OP is refering to this post. $\endgroup$
    – xzczd
    Mar 10 '20 at 7:31
  • $\begingroup$ Sir when I change domain tp [0,1] in this link problem its also not working... u[t, x, 1] == 0, v[t, x, 1] == 0, T[t, x, 1] == 1 $\endgroup$ Mar 10 '20 at 7:55
  • 1
    $\begingroup$ Just remove Quiet. It's only purpose is to suppress error and warning messages. $\endgroup$ Mar 10 '20 at 8:29
2
$\begingroup$

In version 12, the solution method (DAE) is automatically selected

pde1 = D[u[t, x, r], x] + D[v[t, x, r], r] + v[t, x, r]/r == 0;
pde2 = (D[u[t, x, r], t] + u[t, x, r]*D[u[t, x, r], x] + 
     v[t, x, r]*D[u[t, x, r], r] == 
    T[t, x, r] + (1/(1 + lm))*(D[u[t, x, r], r, r] + (1/r)*
         D[u[t, x, r], r] + 
        B*(D[u[t, x, r], r, r, t] + 
           D[v[t, x, r], r]*D[u[t, x, r], r, r] + 
           v[t, x, r]*D[u[t, x, r], r, r, r] + 
           D[u[t, x, r], r]*D[u[t, x, r], x, r] + 
           u[t, x, r]*D[u[t, x, r], x, r, r] + (1/r)*
            D[u[t, x, r], r, t] + v[t, x, r]*D[u[t, x, r], r, r]/r + 
           u[t, x, r]*D[u[t, x, r], x, r]/r)));
pde3 = (D[T[t, x, r], t] + u[t, x, r]*D[T[t, x, r], x] + 
     v[t, x, r]*
      D[T[t, x, r], r] == (1/Pr)*(D[T[t, x, r], r, r] + 
       D[T[t, x, r], r]/r));

(*with initial and boundary conditions:*)
lb = 10; B = .5; Pr = .7; lm = .1;
ics = {u[0, x, r] == 0, v[0, x, r] == 0, T[0, x, r] == 0};
bcs = {u[t, 0, r] == 0, v[t, 0, r] == 0, T[t, 0, r] == 0, 
   u[t, x, 1] == 0, v[t, x, 1] == 0, T[t, x, 1] == 1 - Exp[-5 t], 
   u[t, x, lb] == 0, (D[u[t, x, r], r] /. r -> lb) == 0, 
   v[t, x, lb] == 0, T[t, x, lb] == 0};


{U, V, Ts} = 
  NDSolveValue[{pde1, pde2, pde3, bcs, ics}, {u, v, T}, {t, 0, 
     10}, {x, 0, 1}, {r, 1, lb}] // Quiet;

Visualisation

With[{t = 10}, {DensityPlot[U[t, x, r], {x, 0, 1}, {r, 1, 10}, 
   PlotRange -> All, ColorFunction -> "Rainbow", 
   PlotLegends -> Automatic, PlotLabel -> "u"], 
  DensityPlot[Ts[t, x, r], {x, 0, 1}, {r, 1, 10}, PlotRange -> All, 
   ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   PlotLabel -> "T"]}]

{Plot[U[t, 1, 1.19], {t, 0, 10}, AxesLabel -> Automatic, 
  PlotLabel -> u], 
 Plot[U[10, 1, r], {r, 1, 10}, AxesLabel -> Automatic, 
  PlotLabel -> u], 
 Plot[Ts[t, 1, 1.19], {t, 0, 10}, AxesLabel -> Automatic, 
  PlotLabel -> T], 
 Plot[Ts[10, 1, r], {r, 1, 10}, AxesLabel -> Automatic, 
  PlotLabel -> T]} 

Fig 1

$\endgroup$
4
  • $\begingroup$ thank you so much, dear Alex Trounev for your helpful comment. $\endgroup$ Mar 11 '20 at 4:55
  • $\begingroup$ Sir Why we take boundary condition T[t, x, 1] == 1 - Exp[-5 t]insted of T[t, x, 1] == 1 $\endgroup$ Mar 11 '20 at 5:10
  • $\begingroup$ @Mathematicain This is a damping factor for matching initial and boundary conditions. This is similar to parameter "ScaleFactor" in the"MethodOfLines". $\endgroup$ Mar 11 '20 at 11:31
  • $\begingroup$ Dear Alex Trounev, thank you so much. $\endgroup$ Mar 11 '20 at 12:08

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