2
$\begingroup$

Bug introduced in 8.0.4 or earlier.


I try to solve the heat transfer equation with boundary condition that depends on time:

r0 = 0.75 10^-3;(*Beam spot size, m*)
ω = π ν;
ν = 1; (*pulse repetition rate, Hz*)

c = 1710;(*Heat capacity, W/(m·K)*)
ρ = 879; (*Density, kg/m^3*)
λ = 0.111;(*Heat conductivity, W/(m·K)*)
T0 = 300;(*initial temperature,K*)
T1 = 1000 - T0; (*Hot state temperature, K*)
Rm = 3 10^-3;(*Sample Radius, m*)
zm = 2 10^-3 ;(*Sample thickness, m*)   

 eq = D[T[R, z, t], t] == λ/(
    c ρ) (1/(R + 10^-20) D[R D[T[R, z, t], R], R] + 
      D[T[R, z, t], {z, 2}]);

init1 = T[R, z, 0] == T0;
bc1 = D[T[R, z, t], {R, 1}] == 0 /. R -> 0;
bc2 = D[T[R, z, t], {R, 1}] == 0 /. R -> Rm;
bc3 = T[R, z, t] == 
    Piecewise[{{T0 + T1 (1 - Abs@Sin[ω t])^500, 
       0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4}, 
  T[R, z, t], {R, 0, Rm}, {z, 0, zm}, {t, 0, 10},
  AccuracyGoal -> 30, MaxStepFraction -> 0.05]

The piecewise function defines the temperature at the left side of the z-interval as a short square pulses coming with frequency ω to central part of the disk.

However, solver returns the error-messages:

NDSolveValue::ibcinc: Warning: boundary and initial conditions are inconsistent. NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.

What's wrong? I've tried with FEM package but it produce even more error-messages :)

$\endgroup$
  • 1
    $\begingroup$ When t=0, init and bc3 are inconsistent, it makes T[R, 0, 0] both 300 and 1000 for 0 <= R <= 0.00075. $\endgroup$ – Feyre Feb 15 '17 at 9:43
  • $\begingroup$ @Feyre, ok, let's add small shift to t at bc3 making (t-5*10^-3) instead of t inside Sin. Unfortunately, it does not change principally anything. $\endgroup$ – Rom38 Feb 15 '17 at 10:42
  • $\begingroup$ Why don't change Sin to Cos in bc3? $\endgroup$ – zhk Feb 15 '17 at 10:51
  • $\begingroup$ @MMM, Is it principal? :) $\endgroup$ – Rom38 Feb 15 '17 at 10:53
  • $\begingroup$ @Rom38 I don't know about that. $\endgroup$ – zhk Feb 15 '17 at 10:54
4
$\begingroup$

Several issues here.

  1. As mentioned by Feyre, the i.c. and b.c. are inconsistent. I'll set "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100} in response. For more information check this post.

  2. ndnum is caused by a bug of NDSolve (yeah it's confirmed, I reported it WRI before), similar problem has been observed in this and this post. This can be resolved by adding Simplify`PWToUnitStep@ before Piecewise.

  3. The modeling for laser pulse is improper. Ideally DiracDelta is a possible choice, but currently NDSolve can't handle it properly, so we need to use an approximate one, for example:

    dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];
    

    Periodity can be achieved with Mod:

    Plot[dirac[Mod[t, 1, -1/2], 100], {t, 0, 10}, PlotRange -> All]
    

    Mathematica graphics

Here's the solution:

(* Previous code is not modified so I'd like to omit them in this post. *)
dirac[r_, a_] = Sqrt[a/Pi] Exp[-a r^2];

bc3 = T[R, z, t] == Simplify`PWToUnitStep@
         Piecewise[{{T0 + T1*dirac[Mod[t, Pi/ω, -2^(-1)], 100], 
        0 <= R <= r0}, {T0, True}}] /. z -> 0;
bc4 = T[R, z, t] == T0 /. z -> zm;

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

sol = NDSolveValue[{eq, init1, bc1, bc2, bc3, bc4}, 
    T, {R, 0, Rm}, {z, 0, zm}, {t, 0, 10}, 
    Method -> Union[mol[60, 2], mol[True, 100]]]; // AbsoluteTiming

Animate[Plot3D[sol[r, z, t], {r, 0, Rm}, {z, 0, zm}, PlotRange -> {0, 4000}, 
  PerformanceGoal -> "Quality"], {t, 0, 10}]

enter image description here

Remark

  1. NDSolve spits out eerr warning, but I think it's not a big problem since the error is small. Use more grid points will probably suppress the error, you can have a try if you have time.

  2. I think "FiniteElement" will do a better job on solving this problem.

  3. I took away AccuracyGoal because of the reason mentioned here.

$\endgroup$
  • $\begingroup$ Thanks a lot! I'm enough new in solution of DE in Mathematica.. $\endgroup$ – Rom38 Feb 16 '17 at 5:46

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.