1
$\begingroup$

I'm trying to numerically simulate heat flux in a 3d power cable model that i build. Essentially the equation is:

Temperature and Power equations

I can get an analytic solution from Mathematica easily and when i plot T(x) i get no trouble (putting the constraint T'(x)=0 and T'(l)=0):

enter image description here

second picture shows T'(x)=0

enter image description here

Because i want to add some nonlinear thermal behavior i had to use NDSlove, and first i try to solve the same equation that i wrote before to see if it works:

sol = NDSolve[{-T[x] + \[Tau] + 
 P[x]*(Subscript[\[Rho], t]) + (Subscript[\[Rho], 
   t])/(Subscript[\[Rho], c])*T''[x] == 0, T'[0] == 0,  T'[l] == 0}, T[x], {x, 0, l}]

and the numerical solution give me this:

enter image description here

i can't undertand what's the deal here and why it doesn't converge!!

Here complete Code:

P[x_] := 3*r*(\[Kappa] + \[Lambda]*((x - l/2)/(l/2))^2)^2
T0 = P[0]*Subscript[\[Rho], t]
\[Kappa] = 1597.14
\[Lambda] = 1788 - \[Kappa]
r = 0.0000133
\[Tau] = 20
Subscript[\[Rho], t] = 0.54877
Subscript[\[Rho], c] = 0.332502
l = 100000
\[Alpha] = 0
sol = NDSolve[{-h[x] + \[Tau] + 
     P[x]*(Subscript[\[Rho], t]) + (Subscript[\[Rho], 
       t])/(Subscript[\[Rho], c])*h''[x] == 0, h'[0] == 0, 
   h'[l] == 0}, h[x], {x, 0, l}]
$\endgroup$
  • $\begingroup$ Works for me with my chosen parameters: i.stack.imgur.com/BP3lR.png -- You'll need to give complete code to get much help. $\endgroup$ – Michael E2 Nov 23 '18 at 14:39
  • $\begingroup$ complete code added $\endgroup$ – Mattia Nov 23 '18 at 14:47
  • $\begingroup$ Parameters not defined. How can we verify the solution? $\endgroup$ – Alex Trounev Nov 23 '18 at 15:14
  • $\begingroup$ In your picture analytical solution does not satisfy boundary conditions T'[0] == 0, T'[l] == 0 $\endgroup$ – Alex Trounev Nov 23 '18 at 15:33
  • $\begingroup$ yes it's satisfied, slope of T(x) goes to zero at 0 and l but you can't see it from the picture because T''[x] is zero at 13 and graphic is from 0 to 100000, this is the reason why you can't see it. anyway i need help! $\endgroup$ – Mattia Nov 23 '18 at 15:49
1
$\begingroup$

1)Increase the order of the equation to solve the Dirichlet problem; 2) map the solution to the interval (0,1) by replacing x->l*x; 2) use the method of the false transient. Then the system of equations and the solution are

P[x_] := 3*r*(\[Kappa] + \[Lambda]*(2*x - 1)^2)^2

\[Kappa] = 1597.14;
\[Lambda] = 1788 - \[Kappa];
r = 0.0000133;
\[Tau] = 20;
rt = 0.54877;
rc = 0.332502;
l = 100000;
k = rt/rc/l^2;

eq = {-h[x, t] + D[P[x], x]*rt + k*D[h[x, t], x, x] - 
     D[h[x, t], t] == 0};
ic = h[x, 0] == 0;
bc = {DirichletCondition[h[x, t] == 0, x == 0], 
   DirichletCondition[h[x, t] == 0, x == 1]};

sol = NDSolveValue[{eq, ic, bc}, h, {x, 0, 1}, {t, 0, 20}] 
Plot3D[sol[x, t], {x, 0, 1}, {t, 1, 20}, PlotRange -> All, 
 Mesh -> None, ColorFunction -> Hue, 
 AxesLabel -> {"x", "t", "\!\(\*SubscriptBox[\(h\), \(x\)]\)"}]

Plot[Table[sol[x, t], {t, 5, 20, 5}], {x, 0, 1}, PlotRange -> All, 
 AxesLabel -> {"x", "\!\(\*SubscriptBox[\(h\), \(x\)]\)"}]

fig1

The solution of the problem reaches a stationary state already at t>1. Integrating over x, we find the solution to the original problem up to a constant

hS[x_] := NIntegrate[sol[x1, 20], {x1, 0, x}]

{Plot[hS[x], {x, 0, .01}], Plot[hS[x], {x, 1 - .01, 1}], 
 Plot[hS[x], {x, 0, 1},AxesLabel -> {"x", "T"}]}

fig2

We choose the integration constant = 90, then the solution is similar to the analytical one. Note that with numerical integration on the right-hand boundary, when x-> 1, a numerical instability arises, although the function sol[x1,20] is rather smooth. fig3

$\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.