How to solve a 3D permanent conjugate heat transfer problem?

I want to solve this conjugate heat transfer problem on NDSolve but Mathematica always quits the Kernel. Can you help me? I do not know exactly how to use the numerical methods that are NDSolve method options.
Besides, does it takes too long to solve a problem like that one? I want it in order to compare solutions.

k[x_, y_, z_]:=Piecewise[{{1, 0.25 <= x <= 0.75 ∧ 0.25 <= y <= 0.75 ∧ 0 <= z <= 5}}, 0.25]
w[x_, y_, z_]:=Piecewise[{{1, 0.25 <= x <= 0.75 ∧ 0.25 <= y <= 0.75 ∧ 0 <= z <= 5}}, 0.339266]
Bi = 1;

U[x_, y_] := Piecewise[{{(
     74076188 Sin[(491701844 (-(1/4) + x))/78256779] Sin[(
       491701844 (-(1/4) + y))/78256779])/31356257 + (
     67207613 Sin[(491701844 (-(1/4) + x))/26085593] Sin[(
       491701844 (-(1/4) + y))/78256779])/426732107 + (
     67207613 Sin[(491701844 (-(1/4) + x))/78256779] Sin[(
       491701844 (-(1/4) + y))/26085593])/426732107, 
    0.25 <= x <= 0.75 ∧ 0.25 <= y <= 0.75}}, 0]

solnum = NDSolve[{w[x, y, z] U[x, y] D[θ[x, y, z], z] == 
    Div[k[x, y, z] Grad[θ[x, y, z], {x, y, z}], {x, y, 
      z}], (D[θ[x, y, z], x] /. x -> 0) == 
    0, (D[θ[x, y, z], x] /. x -> 1) == 
    0, (-D[θ[x, y, z], y] /. y -> 0) == 
    1, ((D[θ[x, y, z], y] /. y -> 1) + Bi θ[x, 1, z]) ==
     0, θ[x, y, 0] == 1, (D[θ[x, y, z], z] /. z -> 5) ==
     0}, θ[x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 5}, 
  MaxStepSize -> 0.01, AccuracyGoal -> 3, PrecisionGoal -> 3]

Export["solnum.mx", solnum];

solnum = NDSolve[{w[x, y, z] U[x, y] D[[Theta][x, y, z], z] == 
    Inactive[Div][
      k[x, y, z] Inactive[Grad][[Theta][x, y, z], {x, y, z}], {x, y, 
       z}] + NeumannValue[1, y == 0] + 
     NeumannValue[Bi [Theta][x, y, z], y == 1],
   [Theta][x, y, 0] == 1
   }, [Theta][x, y, z], {x, 0, 1}, {y, 0, 1}, {z, 0, 5}
  ]