FEM nonlinear anisotropic heat transfer with element markers

Question

I'm trying to model a 2D nonlinear anisotropic heat transfer using FEM with element markers. I got stuck with errors and have no idea what to do to make it work...

I have simple toy-model: two connected squares with bottom side of the right heated through a constant heat flux and a left side of a left square is kept at constant temperature. Each of the squares has different material properties. From the FEM perspective they are distinguished through ElementMarkers. I've built this model using procedure described in the Heat Transfer tutorial (PDEModels/tutorial/HeatTransfer/HeatTransfer in Mathematica help) in "Modeling Heat Sources using Element Markers" section. Heat transfer FEM solving equations, as well as all the methods described below, come from the same Heat Transfer tutorial. I'm using Mathematica 12.1.

According to the Heat Transfer tutorial simulating anisotropic heat transfer is as easy as inputing k as a matrix rather than a single value. It works well with constants (you just need to remember to define k separately and not in "parameters" replace rules' set).

Again, according to Heat Transfer tutorial simulating nonlinear heat transfer (changing properties of a material with respect to its temperature) is as easy as defining the property as a function depending on temperature. And indeed it works well.

Finally, to give different properties to different parts of the model you just need a piecewise function with ElementMarkers. Again, on its own it works great.

The problem comes when combining all three methods. I tried combinations of two of the three methods and they all work - you can do nonlinear anisotropic without element markers (same properties for the whole simulation), nonlinear with markers, and anisotropic with markers. But when you try all three it gives errors.

Here is the code for the whole simulation and errors (errors as a picture as well as it is nicer to look at image than the code). Can you give me a clue what to do?

Needs["NDSolve`FEM`"]

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k},
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*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]

coordinates = {
   {0, 0}, {0, 1}, {1, 1}, {1, 0},
   {1, 0}, {1, 1}, {2, 1}, {2, 0}
   };
el = LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 1}, {5, 6}, {6, 7}, {7, 
     8}, {8, 5}}];
bmesh = ToBoundaryMesh["Coordinates" -> coordinates, 
   "BoundaryElements" -> {el}];
bmesh["Wireframe"]

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.5, 0.5}, 2}}];
mesh["Wireframe"[
  "MeshElementStyle" -> {Directive[FaceForm[Red]], 
    Directive[FaceForm[Green]]}]]

\[CapitalGamma]temp = DirichletCondition[T[x, y] == 300, x == 0];
\[CapitalGamma]heater = NeumannValue[1000, y == 0 && 1 <= x <= 2];

k2[T_] := 600 - T;

parameters = {
   \[Rho] -> 
    Piecewise[{{1, ElementMarker == 1}, {2, ElementMarker == 2}}],
   Cp -> Piecewise[{{3, ElementMarker == 1}, {4, ElementMarker == 2}}]
   };
kSet = Piecewise[{{{{5, 0}, {0, 50}}, 
     ElementMarker == 1}, {{{k2[T[x, y]], 0}, {0, 6}}, 
     ElementMarker == 2}}];

pde = {HeatTransferModel[T[x, y], {x, y}, kSet, \[Rho], Cp, "NoFlow", 
      "NoSource"] == \[CapitalGamma]heater, \[CapitalGamma]temp} /. 
   parameters;
Tfun = NDSolveValue[pde, T, {x, y} \[Element] mesh]

InitializePDECoefficients::femcnmd: The PDE coefficient {-(-(\[Piecewise]   {<<2>>} Equal[<<2>>]
{<<2>>} Equal[<<2>>]
0   True

)).{(T^(1,0))[x,y],(T^(0,1))[x,y]},-(-(\[Piecewise] {<<2>>} Equal[<<2>>]
{<<2>>} Equal[<<2>>]
0   True

)).{(T^(1,0))[x,y],(T^(0,1))[x,y]}} does not evaluate to a numeric matrix of dimensions {2,1} at the coordinate {0.295534,0.46115}; it evaluated to {{-0.5,-5.},{-0.5,-5.}} instead.

NDSolveValue::fempslf: The linearization process in PDESolve failed.

Answer 1

There seems to be a problem with the linearization of the PDE when it contains a Piecewise.

Try to use

ClearAll[HeatTransferModel]
HeatTransferModel[T_, X_List, k_, \[Rho]_, Cp_, Velocity_, Source_] :=
  Module[{V, Q, a = k}, 
  V = If[Velocity === "NoFlow", 
    0, \[Rho]*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]

with this kSet

kSet = {{If[ElementMarker == 1, 5, k2[T[x, y]] // Evaluate], 0}, {0, 
   If[ElementMarker == 1, 50, 6]}}

This works in V12.3 and V12.2 - I no longer have 12.1 installed.