Heat Conduction where thermal diffusivity depends on z

Question

I'd like to know why Mathematica is not solving this problem for me. Let's imagine a body in which the thermal diffusivity varies with z (we have a cyllindrical simmetry). I want to know what is the temperature distribution in a time t, at radius r and depth z.

Needs["NDSolve`FEM`"]

tdiff[z_] := 0.5 + 0.1 UnitStep[z - 1] + 0.2 UnitStep[z - 2];

eqn = 
  D[u[r, z, t], z, z] + (1/r) D[r D[u[r, z, t], r], r] - (tdiff[z]) D[u[r, z, t], t];

Ω = 
  ImplicitRegion[True, {{r, 0.001, 200}, {z, 0, 3}, {t, 0, 1000}}];

sol = 
  NDSolve[
    {eqn == NeumannValue[30 - u[r, z, t], z == 0], 
     DirichletCondition[ u[r, z, t] == 35, z == 3], 
     DirichletCondition[u[r, z, t] == 0, r == 200], 
     DirichletCondition[u[r, z, t] == 35 + 5 Exp[-r^2] Exp[-z], t == 0], 
     D[u[r, z, t], t] == 0 /. t -> 0}, 
    u, {r, z, t} ∈ Ω]

Mathematica takes forever and does not solve the problem. Why is that?

Answer 1

There are a few issues with this: You give {r,z,t} as one region, also the initial condition is given as a DirichletCondition this leads NDSolve to treat this as a 3D example; that however, leads to problems as then the PDE is convection dominated. I also changed the equation to what I think you want, but check that. So this will work:

tdiff[z_] := 0.5 + 0.1 UnitStep[z - 1] + 0.2 UnitStep[z - 2];

eqn = tdiff[z]*
    D[u[t, r, z], 
     t] - (D[u[t, r, z], z, z] + (1/r) D[r D[u[t, r, z], r], r]);
Ω = ImplicitRegion[True, {{r, 1/1000, 200}, {z, 0, 3}}];

sol = NDSolveValue[{eqn == NeumannValue[30 - u[t, r, z], z == 0],
   DirichletCondition[u[t, r, z] == 35, z == 3],
   DirichletCondition[u[t, r, z] == 0, r == 200],
   u[0, r, z] == 0
   }, u, {t, 0, 1000}, {r, z} ∈ Ω]