2
$\begingroup$

I'm trying to solve a system of PDE using NDsolve but the following error accors.

NDSolveValue::femcmsd: The spatial derivative order of the PDE may not exceed two.

Equations: 1 2

Qm and Q are constant.

    pde1 = {
   D[qr[r, x, t], t] + tq*D[D[qr[r, x, t], t], t] ==  
    alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), r] + 
     alpha*ta*
      D[D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), r], t] + 
     alpha*wrc*D[Te[r, x, t], r] + alpha*wrc*ta*D[D[Te[r, x, t], t], r]
   };

    pde2 = {
   D[qx[r, x, t], t] + tq*D[D[qx[r, x, t], t], t] ==  
    alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x] + 
     alpha*ta*
      D[D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x], t] + 
     alpha*wrc*D[Te[r, x, t], x] + alpha*wrc*ta*D[D[Te[r, x, t], t], x]
   };

    pde3 = {
   ro*c*D[Te[r, x, t], 
      t] == -(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]) + 
     wrc*(Tb - Te[r, x, t]) + Qm
   };

    sol = NDSolveValue[{pde1, pde2, pde3, BC1, BC2, BC3, BC4, BC5, IC1, 
    IC2, IC3, IC4, IC5, IC6}, 
   Te[r, x, t], {t, 0, 5}, {r, x} ∈ Ω];

The BS's and IC's :

when[t <= tp] {
   BC1 = DirichletCondition[qx[r, x, t] == Qp, x == 0];};
when[t > tp] {
    BC2 = DirichletCondition[qx[r, x, t] == 0, x == 0];};

BC3 = DirichletCondition[qx[r, x, t] == 0, x == xo];
BC4 = DirichletCondition[qx[r, x, t] == 0, r == 0];
BC5 = DirichletCondition[qx[r, x, t] == 0, r == ro];

IC1 = qx[r, x, 0] == 0;
IC2 = qr[r, x, 0] == 0;
IC3 = Te[r, x, 0] == To;

IC4 = Derivative[0, 0, 1][qx][r, x, 0] == 0;
IC5 = Derivative[0, 0, 1][qr][r, x, 0] == 0;
IC6 = Derivative[0, 0, 1][Te][r, x, 0] == 0;

Subscript[r, 1] = 0;
Subscript[r, 2] = 0.05;
Ω = 
  ImplicitRegion[
   True, {{r, Subscript[r, 1], Subscript[r, 2]}, {x, 0, 0.05}}];

Is there any way to handle this errorr?

link of the article: https://www.sciencedirect.com/science/article/pii/S1290072908002937

Improve this question
$\endgroup$
7
  • 4
    $\begingroup$ Hard to say without information on what the BCs and ICs and Omega look like. $\endgroup$
    – user21
    Commented Apr 29, 2019 at 4:45
  • $\begingroup$ also alpha,ta,wrc,ro,Qm,c ? $\endgroup$
    – Mariusz Iwaniuk
    Commented Apr 29, 2019 at 14:09
  • $\begingroup$ @MariuszIwaniuk Yes, they are constant coefficients $\endgroup$
    – Hadi
    Commented May 4, 2019 at 5:35
  • 1
    $\begingroup$ @Hadi Give a link to the article where this model is published. $\endgroup$
    – Alex Trounev
    Commented May 4, 2019 at 15:38
  • 1
    $\begingroup$ @Hadi Add the parameters that are defined in the article: alpha=?; ta=?; wrc=?; ro=?; Qm=?;ta=?;tq?. $\endgroup$
    – Alex Trounev
    Commented May 11, 2019 at 12:19

1 Answer 1

Reset to default
5
$\begingroup$

I debugged the code using equations (5), (7), (8) from the article Jianhua Zhou, Yuwen Zhang, J.K. Chen, An axisymmetric dual-phase-lag bioheat model for laser heating of living tissues. Unfortunately, I was not able to achieve a match with the data shown in Fig. 3, since I could not find the parameters they used. I give the results as it is with the set of parameters that I found in the article (some mixture of tests 1 and 2).

x1 = 1; r0 = 0; r1 = 2; \[CapitalOmega] = 
 ImplicitRegion[0 <= x <= x1 && r0 <= r <= r1, {r, x}];
alpha = 1; T0 = 0; Tb[t_] := (1 - Exp[-100*t]); Qm[t_] = 0; 
Q[t_] = 0; tq = .05; ta = .001; c = 4187; rho = 1000; rhob = 
 1.06*10^3; cb = 3860; wb = 1.87*10^\[Minus]3 ; b = 
 tq*wb*rhob*cb/(rho*c); wrc = 
 wb*rhob*cb/(rho*
     c);(*\[Rho]=1000 kg/m3,k=0.628 W/(m K),c=4187 J/(kg K) for the \
thermophysical properties;\[Rho]b=1.06\[Times]10^3 kg/m3,cb=3860 \
J/(kg K),wb=1.87\[Times]10^\[Minus]3 m3/(m3 tissue s),Tb=37 \
\[EmptySmallCircle]C,Q m=1.19\[Times]10^3 W/m3*)
pde1 = qr1[r, x, t] + tq*D[qr1[r, x, t], t] == 
  alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), r] + 
   alpha*ta*
    D[(1/r)*(D[r*qr1[r, x, t], r] + D[r*qx1[r, x, t], x]), x] + 
   alpha*wrc*D[T[r, x, t], r] + alpha*wrc*ta*D[T1[r, x, t], r];

pde2 = qx1[r, x, t] + tq*D[qx1[r, x, t], t] == 
   alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x] + 
    alpha*ta*
     D[(1/r)*(D[r*qr1[r, x, t], r] + D[r*qx1[r, x, t], x]), x] + 
    alpha*wrc*D[T[r, x, t], x] + alpha*wrc*ta*D[T1[r, x, t], x];

pde3 = tq*D[T1[r, x, t], t] + (1 + b)*T1[r, x, t] == 
   alpha*(Laplacian[T[r, x, t], {r, x}] + D[T[r, x, t], r]/r) + 
    alpha*ta*(Laplacian[T1[r, x, t], {r, x}] + D[T1[r, x, t], r]/r) + 
    wrc*(T0 - T[r, x, t]) + b*D[Tb[t], t] + (Qm[t] + Q[t])/(rho*c) + 
    tq*D[Qm[t] + Q[t], t]/(rho*c);
pde4 = {qr1[r, x, t] == D[qr[r, x, t], t], 
   qx1[r, x, t] == D[qx[r, x, t], t], T1[r, x, t] == D[T[r, x, t], t]};

bc = {DirichletCondition[{qx[r, x, t] == 0, qr[r, x, t] == 0, 
     qx1[r, x, t] == 0, qr1[r, x, t] == 0, T[r, x, t] == 0, 
     T1[r, x, t] == 0}, x == x1 || r == r1], 
   DirichletCondition[T[r, x, t] == Tb[t]*If[r <= 1, 1, 0], x == 0], 
   DirichletCondition[{qx[r, x, t] == 0, qr[r, x, t] == 0, 
     qx1[r, x, t] == 0, qr1[r, x, t] == 0}, r == r0]};

ic = {qx[r, x, 0] == 0, qr[r, x, 0] == 0, T[r, x, 0] == 0, 
   qx1[r, x, 0] == 0, qr1[r, x, 0] == 0, T1[r, x, 0] == 0};

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{qrs, qr1s, qxs, qx1s, Ts, T1s} = 
   NDSolveValue[{pde1, pde2, pde3, pde4, bc, ic}, {qr, qr1, qx, qx1, 
     T, T1}, {r, r0, r1}, {x, 0, x1}, {t, 0, .05}, 
    Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> "TensorProductGrid"}}, 
    EvaluationMonitor :> (currentTime = t;)];]


{Plot3D[Ts[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All, AxesLabel -> {"r", "x", "T"}], 
 Plot3D[T1s[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All, AxesLabel -> {"r", "x", "T1"}]}
Plot[Ts[r0, x, .05], {x, 0, 1}, PlotRange -> All, 
 AxesLabel -> {"x", "T"}]

{Plot3D[qxs[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All], 
 Plot3D[qrs[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All], 
 Plot3D[qx1s[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All], 
 Plot3D[qr1s[r, x, .05], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All]}

fig1

I debugged OP's code using equations (4), (7), (8) from the article Jianhua Zhou, Yuwen Zhang, J.K. Chen, An axisymmetric dual-phase-lag bioheat model for laser heating of living tissues. Here is a test with a volumetric source.The model is unstable. Waves develop even with a constant source.

x1 = .05; r0 = 0; r1 = .05; \[CapitalOmega] = 
 ImplicitRegion[0 <= x <= x1 && r0 <= r <= r1, {r, x}];
alpha = 1; av = 0.01; t0 = .9; T0 = 37; 
qp[t_] := 0*150*(1 - Exp[-10*t]); mua = 100; 
Qm[r_, x_, t_] := 1.19*10^3; 
Q[r_, x_, t_] := 
 qp[t]*mua*Exp[-mua*x]; tq = 1; ta = .5; c = 4187; rho = 1000; rhob = 
 1.06*10^3; cb = 3860; wb = 1.87*10^\[Minus]3 ; b = 
 tq*wb*rhob*cb/(rho*c); wrc = 
 wb*rhob*cb;(*\[Rho]=1000 kg/m3,k=0.628 W/(m K),c=4187 J/(kg K) for \
the thermophysical properties;\[Rho]b=1.06\[Times]10^3 kg/m3,cb=3860 \
J/(kg K),wb=1.87\[Times]10^\[Minus]3 m3/(m3 tissue s),Tb=37 \
\[EmptySmallCircle]C,Q m=1.19\[Times]10^3 W/m3*)
pde1 = {D[qr[r, x, t], t] + tq*D[D[qr[r, x, t], t], t] == 
   alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x] + 
    alpha*ta*
     D[D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x], t] + 
    alpha*wrc*D[Te[r, x, t], r] + 
    alpha*wrc*ta*D[D[Te[r, x, t], t], r] - alpha*D[Q[r, x, t], r] - 
    alpha*ta*D[Q[r, x, t], r, t]};

pde2 = {D[qx[r, x, t], t] + tq*D[D[qx[r, x, t], t], t] == 
    alpha*D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x] + 
     alpha*ta*
      D[D[(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]), x], t] + 
     alpha*wrc*D[Te[r, x, t], x] + 
     alpha*wrc*ta*D[D[Te[r, x, t], t], x] - alpha*D[Q[r, x, t], x] - 
     alpha*ta*D[Q[r, x, t], x, t]};

pde3 = {rho*c*
     D[Te[r, x, t], 
      t] == -(1/r)*(D[r*qr[r, x, t], r] + D[r*qx[r, x, t], x]) + 
     wrc*(T0 - Te[r, x, t]) + Qm[r, x, t] + Q[r, x, t]};

bc = {DirichletCondition[{qx[r, x, t] == 0, qr[r, x, t] == 0}, 
    x == x1 || r == r1], 
   DirichletCondition[{qx[r, x, t] == 0, qr[r, x, t] == 0}, x == 0], 
   DirichletCondition[{qx[r, x, t] == 0, qr[r, x, t] == 0}, r == r0]};

ic = {qx[r, x, 0] == 0, qr[r, x, 0] == 0, 
   Derivative[0, 0, 1][qx][r, x, 0] == 0, 
   Derivative[0, 0, 1][qr][r, x, 0] == 0, Te[r, x, 0] == T0};

Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{qrs, qxs, Ts} = 
   NDSolveValue[{pde1, pde2, pde3, bc, ic}, {qr, qx, Te}, {r, r0, 
     r1}, {x, 0, x1}, {t, 0, t0}, 
    Method -> {"IndexReduction" -> Automatic, 
      "EquationSimplification" -> "Residual", 
      "PDEDiscretization" -> {"MethodOfLines", 
        "SpatialDiscretization" -> "TensorProductGrid"}}, 
    EvaluationMonitor :> (currentTime = t;)];]
{Plot3D[qxs[r, x, t0], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All, AxesLabel -> {"r", "x", "qx"}, Mesh -> None, 
  ColorFunction -> Hue], 
 Plot3D[qrs[r, x, t0], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All, AxesLabel -> {"r", "x", "qr"}, Mesh -> None, 
  ColorFunction -> Hue], 
 Plot3D[Ts[r, x, t0], {r, x} \[Element] \[CapitalOmega], 
  PlotRange -> All, AxesLabel -> {"r", "x", ""}, PlotLabel -> T, 
  Mesh -> None, ColorFunction -> Hue]}

fig2

Improve this answer
$\endgroup$

Your Answer

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

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