3
$\begingroup$

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

$\endgroup$
9
  • $\begingroup$ $length$ is not defined. I could not find GraphSurfaceMesh in the documentation. You probably should be using ToElementMesh to creates meshes to be solved by NDSolve. $\endgroup$
    – Tim Laska
    Aug 12, 2020 at 19:27
  • $\begingroup$ You are also missing the polystyrene parameters. $\endgroup$
    – Tim Laska
    Aug 12, 2020 at 19:35
  • 1
    $\begingroup$ Perhaps h=length. I even have no found GraphSurfaceMesh on the internet. InfinitePlane wants the triple as a List. The list {x,y,z} in the pde causes problems. The indexed variables in the boundary condition do not work that way. NDSolve requires [values, Reals. Values for the plastics are missing. mathematica.stackexchange.com/questions/55546/… might be helpful. $\endgroup$ Aug 12, 2020 at 19:50
  • $\begingroup$ sorry guys, I forgot to copy and paste those parts of the code, give me a second to edit the post $\endgroup$
    – kmulc
    Aug 12, 2020 at 19:55
  • $\begingroup$ ok I have added them. Also I did try to use ToELementMesh but it was not showing up, which is why I used GraphSurfaceMesh. I found that here under applications: reference.wolfram.com/language/ref/MeshRegion.html $\endgroup$
    – kmulc
    Aug 12, 2020 at 20:02

1 Answer 1

3
$\begingroup$

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

Mesh

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]

Solution animation

$\endgroup$
3
  • $\begingroup$ Thanks for your help. I tried your code but I am getting this error "NDSolveValue::femper: PDE parsing error of {{1306.25 T$188889-NeumannValue[150 (20+Times[<<2>>]),ElementMarker==4]-NeumannValue[150 (20+Times[<<2>>]),ElementMarker==5]}}. Inconsistent equation dimensions. " I would like to get it fixed so that I can visualize it like you said $\endgroup$
    – kmulc
    Aug 13, 2020 at 14:57
  • $\begingroup$ @kmulc I started a new Mathematica session and copied the code for HeatTransferModel from your question and the code from my answer and it runs. I am using "12.1.1 for Microsoft Windows (64-bit) (June 19, 2020)". You may want to Clear variables or close out and restart Mathematica. $\endgroup$
    – Tim Laska
    Aug 13, 2020 at 18:38
  • $\begingroup$ That was the problem. Thank you so much! $\endgroup$
    – kmulc
    Aug 13, 2020 at 21:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.