NDSolve Inconsistent equation dimensions

Question

I am getting an error while trying to solve for a differential equation. It is saying "NDSolveValue::femper: PDE parsing error... Inconsistant equation dimensions." I was wondering if anyone could help me out figuring what I did wrong. My domain is a prism, and I want to maintain a constant temperature on one of the faces, have the entire prism be that temperature at time = 0, and also has a symmetry boundary. The last faces I will put convective heat transfer through. Sorry if this is a bit hard to visualize through text.

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]

above is the PDE function, below is remainder of the code

length = 0.3;

plastic = 
  Prism[{{0.1335, 0, 0.1585}, {0.15, 0, 0.1415}, {0.15, 0, 
     0.1585}, {0.1335, length, 0.1585}, {0.15, length, 0.1415}, {0.15,
      length, 0.1585}}];

mesh = MeshRegion[plastic, PlotTheme -> "Lines"]; 
 GraphSurfaceMesh[{mesh}]

Subscript[T, hot] = 200;
h = 150;

Subscript[\[Rho], polystyrene] = 1045;
Subscript[Cp, polystyrene] = 1.25;
Subscript[k, polystyrene] = 0.14; 

(* boundary conditions *) 
Subscript[Γ, 
   temp] = {DirichletCondition[T[t, x, y, z] == Subscript[T, hot], 
    y >= length ]};
Subscript[Γ, symmetry] = {NeumannValue[0, x == 0.15]};
Subscript[Γ, 
   convective] = {NeumannValue[h*(Subscript[T, cold] - T[t, x, y, z]),
     z == 0.1585]};
Subscript[Γ, 
   convective1] = {NeumannValue[
    h*(Subscript[T, cold] - T[t, x, y, z]), 
    InfinitePlane[{0.1335, 0, 0.1585}, {0.15, length, 0.1415}, {0.15, 
      0, 0.1415}]]};

ic = {T[0, x, y, z] == Subscript[T, hot]};
parameters = {ρ -> Subscript[ρ, polystyrene], 
   Cp -> Subscript[Cp, polystyrene], k -> Subscript[k, polystyrene]};

tend = 30; (* s *)
pde = {TimeHeatTransferModel[T[t, x, y, z], t, {x, y, z}, k, ρ, 
     Cp, "NoFlow", "NoSource"] == 
    Subscript[Γ, symmetry] + 
     Subscript[Γ, convective] + 
     Subscript[Γ, convective1], 
   Subscript[Γ, temp], ic} /. parameters;
measure = 
  AbsoluteTiming[
   MaxMemoryUsed[
     Monitor[Tfun = 
       NDSolveValue[pde, T, {t, 0, tend}, {x, y, z} ∈ mesh, 
        EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], 
      monitor]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]

Answer 1

There were a couple missing/faulty definitions. Namely, $T_{cold}$ and the heat capacity of polystyrene should be 1000x larger at that density. Also note that the parentheses should have been removed from the NeumannValue specification and that a zero flux is the default setting.

After making those modification, you can use the BoundaryElementMarkerUnion property of the mesh to find the element and point markers assigned to the boundaries for easier boundary condition assignment.

Here is a possible workflow:

Needs["NDSolve`FEM`"]
length = 0.3
ρpolystyrene = 1045;
Cppolystyrene = 1250;
kpolystyrene = 0.14;
Thot = 200;
Tcold = 20;
h = 150;
tend = 100;
plastic = 
  Prism[{{0.1335, 0, 0.1585}, {0.15, 0, 0.1415}, {0.15, 0, 
     0.1585}, {0.1335, length, 0.1585}, {0.15, length, 0.1415}, {0.15,
      length, 0.1585}}];
mesh = ToElementMesh[plastic, "MaxBoundaryCellMeasure" -> 0.000002, 
   "MaxCellMeasure" -> 0.00008];
(* Visualize Boundary Markers for easier BC Assignment *)
groups = mesh["BoundaryElementMarkerUnion"]
temp = Most[Range[0, 1, 1/(Length[groups])]]
colors = ColorData["BrightBands"][#] & /@ temp
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors, 
  Axes -> True]]
(*boundary conditions*)
Γtemp = {DirichletCondition[T[t, x, y, z] == Thot, 
    ElementMarker == 2]};
(* Zero Flux is the default Neumann Condition *)
(*Subscript[Γ,symmetry]={NeumannValue[0,x\[Equal]0.15]};\
*)
Γconvective = 
  NeumannValue[h*(Tcold - T[t, x, y, z]), ElementMarker == 5];
Γconvective1 = 
  NeumannValue[h*(Tcold - T[t, x, y, z]), ElementMarker == 4];

ic = {T[0, x, y, z] == Thot};
parameters = {ρ -> ρpolystyrene, Cp -> Cppolystyrene, 
   k -> kpolystyrene};

(*s*)pde = {TimeHeatTransferModel[T[t, x, y, z], t, {x, y, z}, 
      k, ρ, Cp, "NoFlow", 
      "NoSource"] ==(*Subscript[Γ,
     symmetry]+*)Γconvective + \
Γconvective1, Γtemp, ic} /. parameters;
measure = 
  AbsoluteTiming[
   MaxMemoryUsed[
     Monitor[Tfun = 
       NDSolveValue[pde, T, {t, 0, tend}, {x, y, z} ∈ mesh, 
        EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], 
      monitor]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]

Visualization

One could use SliceContourPlot3D to view several clip planes of the solution as it evolves over time.

uRange = MinMax[Tfun["ValuesOnGrid"]];
frames = Table[
   Rasterize@
    SliceContourPlot3D[
     Tfun[t, x, y, z], {y == 0, y == length, x == 0.15, 
      z == 0.1585, {"YStackedPlanes", 2}}, {x, y, z} \[Element] mesh, 
     PlotRange -> uRange, 
     ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
     ContourStyle -> Opacity[0.5`], ColorFunctionScaling -> False, 
     Contours -> 15, Boxed -> False, Axes -> False, PlotPoints -> 50, 
     MaxRecursion -> 4, 
     ViewPoint -> {-0.36984446450781705`, 
       0.7256633889310892`, -3.284300186955811`}, 
     ViewVertical -> {-0.03520209233494262`, 
       0.9145566012946781`, -0.40292311391079266`}, 
     PlotLegends -> Automatic], {t, tend/10, tend, 2}];
ListAnimate[frames]