As mentioned in the comments, one can follow the Heat Transfer Tutorial to quickly set up a variety of heat transfer related problems.

We can copy the code to set up operators for steady-state and transient problems:

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, ρ_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k}, 
  V = If[Velocity === "NoFlow", 
    0, ρ*Cp*Velocity.Inactive[Grad][T, X]];
  Q = If[Source === "NoSource", 0, Source];
  If[FreeQ[a, _?VectorQ], a = a*IdentityMatrix[Length[X]]];
  If[VectorQ[a], a = DiagonalMatrix[a]];
  (*Note the-sign in the operator*)
  a = PiecewiseExpand[Piecewise[{{-a, True}}]];
  Inactive[Div][a.Inactive[Grad][T, X], X] + V - Q]
TimeHeatTransferModel[T_, TimeVar_, X_List, k_, ρ_, Cp_, 
  Velocity_, Source_] := ρ*Cp*D[T, {TimeVar, 1}] + 
  HeatTransferModel[T, X, k, ρ, Cp, Velocity, Source]

We can follow one of the many examples and adapt it to your problem. The default NeumannValue is zero flux, which is the case for insulated walls. I adapted the example from here for your problem.

Ω2D = Rectangle[{0, 0}, {10, 10}];
tend = 30;
parameters = {ρ -> 1, Cp -> 1, k -> 1};
ic = u[0, x, y] == x + y;
pde = {TimeHeatTransferModel[u[t, x, y], t, {x, y}, k, ρ, Cp, 
      "NoFlow", "NoSource"] == 0, ic} /. parameters;
ufun = NDSolveValue[pde, 
   u, {t, 0, tend}, {x, y} ∈ Ω2D];
pRange = MinMax[ufun["ValuesOnGrid"]];
legendBar = 
  BarLegend[{"BrightBands", pRange}, 
   Sequence[50, LegendLabel -> Style["[°C]", Opacity[0.6]]]];
options = {PlotRange -> pRange, 
   Sequence[ColorFunction -> ColorData[{"BrightBands", pRange}], 
    ContourStyle -> None, FrameLabel -> {"x", "y"}, 
    ColorFunctionScaling -> False, 
    PlotLabel -> Style["Transient Temperature Field: u(t,x,y)", 18], 
    AspectRatio -> Automatic, Contours -> 75, PlotPoints -> 41, 
    ImageSize -> Medium]};
nframes = 30;
frames = Table[
   Show[Legended[
     ContourPlot[ufun[t, x, y], {x, y} ∈ Ω2D, 
      Evaluate[options]], legendBar]], {t, 0, tend, tend/nframes}];
frames = Rasterize[#1, "Image", ImageResolution -> 80] & /@ frames;
ListAnimate[Sequence[frames, SaveDefinitions -> True], 
 ControlPlacement -> Top]

Updated Response to Edited Code

The OP edited the code to produce:

heqn = D[u[x, y, t], t] == Laplacian[u[x, y, t], {x, t}]

shape = Rectangle[{0, 0}, {10, 10}]

ic = u[x, y, 0] == x+y

bc = NeumannValue[0, True]

sol = NDSolveValue[{heqn == bc, ic},u, {t, 0, 30}, {x, y} ∈ shape]

First, you should have your independent variables have consistent ordering with the NDSolve specification. Namely, {t, 0, 30}, {x, y} ∈ shape means that the dependent variable should be defined with variables {t,x,y} or u[t,x,y].

Second, your Laplacian specification differentiated {x,t}, when it should be {x,y}.

Third, your heqn already contains ==. For the FEM method to work, you will want to express your equation in coefficient form, or:

$$frac{{{partial ^2}}}{{partial {t^2}}}u + dfrac{partial }{{partial t}}u + nabla cdotleft( { - cnabla u - alpha u + gamma } right) + beta cdotnabla u + au - f = 0$$

So, you pull all terms of the operator to the left hand side and replace the $0$ on the right hand side with a sum of NeumannValues.

The corrected code for your edit should look like this:

op = D[u[t, x, y], t] - Laplacian[u[t, x, y], {x, y}]

shape = Rectangle[{0, 0}, {10, 10}]

ic = u[0, x, y] == x + y

sol = NDSolveValue[{op == 0, ic}, 
  u, {t, 0, 30}, {x, y} ∈ shape]