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]&

u
is bounded atr=0
(center of disk). This causes one of the solutions causing the blow up to be discarded. This is standard method when solving heat PDE on disk. It looks like ndsolve is complaining about this in your example. $\endgroup$ – Nasser Jan 21 '17 at 22:08