# Coupled PDEs with second order spatial derivative

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:

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

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

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]}


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]}