1/0 infinity while solving heat equation in polar coordinates

Question

The ndsolve

    uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

NDSolve[{a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t], (D[W[r, phi, t], r] /. r -> 1) == ud[phi] u[t], W[r, phi,0] ==W0[r, phi]}, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T}]

gives the error

Power::infy: Infinite expression 1/0.^2 encountered.
NDSolveValue::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

No idea on how to proceed. Any hints or ideas?

Answer 1

Caution : There something strange in this solution. Not sure it's OK. The strange thing is that the temperature tends towards 0. This is unlikely, except if the Op's target is precisely to do that.
To be analysed.

Here is the solution :

I works only on Mathematica >= 11

uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

opts = Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};

res=NDSolve[{
a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t]+NeumannValue[ud[phi] u[t],r==1], 
 W[r, phi,0] ==W0[r, phi]
 }, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T},
 opts]

 ContourPlot[ res[[1,1,2]][Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], 10],{x,-1,1},  
 {y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99]]

My approach :

• As I was aware that the 1/0 problem was solved on V11 and fem (See the penultimate paragraph here.), I added the option FiniteElement
• Then I had a new problem with a error message that let me think that your condition (D[W[r, phi, t], r] /. r -> 1) == ud[phi] u[t], which is a Neumann boundary conditon, was interpreted as a Dirichlet boundary condition (I don't know why), so i tried the explicit Neumann syntax and it worked

The solution above corresponds to the case where there is a adiabatic wall at phi=0. Here is the solution without the wall :

Edit from the Author, 01/01/2020 The following is not correct : The wall is not completely removed. The problem is that the flux that enters the wall is not the same than the flux that get outside the wall. A problem of this kind is solved here
END EDIT

uo = 1; m = 1; n = 1;  a = 1; T = 20;  M = 10;
W0[r_, phi_] := Exp[-(1 - r)^2]
ud[phi_] := UnitStep[phi] - UnitStep[phi - Pi/36]
u[t_] := uo (1 - Cos[(2 Pi m)/T t])^n

opts = Method -> {"MethodOfLines", 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.01}}};

res=NDSolveValue[{
a Laplacian[W[r, phi, t], {r, phi},"Polar"]== D[W[r, phi, t],t]+NeumannValue[ud[phi] u[t],r==1], 
PeriodicBoundaryCondition[W[r,phi,t],phi==2 Pi,Function[v,v+{0,-2 Pi}]],
W[r, phi,0] ==W0[r, phi]
 }, W, {r, 0, 1}, {phi, 0, 2 Pi}, {t, 0, T},
 opts]

ContourPlot[ res[Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], 0.3],{x,-1,1},  
{y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99]]

An animation of the same thing :

Labeled[ContourPlot[ res[Sqrt[x^2+y^2],Pi+ArcTan[-x,-y], #],{x,-1,1},  
  {y,-1,1},RegionFunction->Function[{x,y},x^2+y^2<0.99],
ContourLabels->True,Contours-> {
.9,.8613,.861,.860,.858,.856,.855,.850,.845,.840,.835,.830,
.825,.8,.775,.750,.725,.7,.68,.65,.625,.6,.55,.525,
.5,.475,.45,.425,.4,.38,.36,.34,.32,.30,.26,.24,.20,
.175,.15,.125,.10,.075,.05,
.020,1. 10^-3,1. 10^-8},
ColorFunction-> "Temperature",
ColorFunctionScaling->False
],Row[{"t=",#}]]& /@ Range[0,16,0.1] //
Export["therm00.gif",#, "DisplayDurations" -> 0.3]&