3

Bug introduced in 5.0, persisting through 13.2.

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 :)

xzczd
  • 65,995
  • 9
  • 163
  • 468
Rom38
  • 5,129
  • 13
  • 28
  • 1
    When t=0, init and bc3 are inconsistent, it makes T[R, 0, 0] both 300 and 1000 for 0 <= R <= 0.00075. – Feyre Feb 15 '17 at 09:43
  • @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. – Rom38 Feb 15 '17 at 10:42
  • Why don't change Sin to Cos in bc3? – zhk Feb 15 '17 at 10:51
  • @MMM, Is it principal? :) – Rom38 Feb 15 '17 at 10:53
  • @Rom38 I don't know about that. – zhk Feb 15 '17 at 10:54
  • There're several issues here, including a possible bug. Before giving an answer, I want to ask: why does the (1 - Abs@Sin[ω t])^500 term exist? Are you approximate DiracDelta with it? If so, then it's incorrect because these peaks don't integrate to 1, if not, then this term can be simply taken away, because it'll almost have no effect in the solution. – xzczd Feb 15 '17 at 10:54
  • The Sin is defining the time-profile of the pulse. It can be different but principally the solution should be available by the same method. At least I guess so :\ – Rom38 Feb 15 '17 at 10:56
  • @xzczd, this term defines the periodical heat incomes at left boundary. I did not realize how to set the periodical peaks with DiracDelta. Do you have any ideas? The integration errors should lead to certain problems with convergence of the solution but in this case solver do not start at all.. Actually, these peaks can have any area from the common point - the heat incomes can has an arbitrary form. – Rom38 Feb 15 '17 at 11:01
  • I added a "bugs" tag to the question. Though the ndnum warning isn't the only issue behind OP's problem, I think it's good to have a question related to this bug tagged with "bugs", so future visitors can find it easier. – xzczd Feb 15 '17 at 12:32

1 Answers1

5

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.

Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468