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

$\endgroup$
  • 4
    $\begingroup$ Hard to say without information on what the BCs and ICs and Omega look like. $\endgroup$ – user21 Apr 29 at 4:45
  • $\begingroup$ also alpha,ta,wrc,ro,Qm,c ? $\endgroup$ – Mariusz Iwaniuk Apr 29 at 14:09
  • $\begingroup$ @MariuszIwaniuk Yes, they are constant coefficients $\endgroup$ – Hadi May 4 at 5:35
  • 1
    $\begingroup$ @Hadi Give a link to the article where this model is published. $\endgroup$ – Alex Trounev May 4 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 May 11 at 12:19
4
$\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

$\endgroup$

Your Answer

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

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