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

Edit from the author of the answer

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

Edit

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

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