Coupled pde for multiphysics joule heating with temperature dependent properties - error inconsistent dimensions

I want to simulate 3D joule heating and followed the wolfram multiphysics tutorial.

This multiphysics problem is comprised of the equation for the voltage field:

and the heat equation:

with the source term of the joul heating for the heat equation:

While the tutorial is one way coupled (solving voltagefield first and using the solution to calculate the source term of the heat equation), I want to use temperature dependent material properties, meaning G=f(T) and k=f(T). But when I try solving the coupled system (with the temperature dependent electr. conductivity G) I get an error.

I'm sorry, if this has been answered elsewhere, but I couldn't find anything that seemed applicable to my case.

Let's jump into my example case. Note, the interesting part is the creation the coupled pde, I guess.

Creating a mesh:

Needs["NDSolveFEM"];

(*geometry constants*)
xmin = 0; xmax = 0.1;
ymin = 0; ymax = 0.2;
l = 1;

(*mesh generation*)
cubi = Cuboid[{xmin, ymin, 0}, {xmax, ymax, l}];
mesh = ToElementMesh[cubi, MaxCellMeasure -> 0.00003]

Show[mesh[
"Wireframe"["MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm[LightBlue], Axes -> True,
AxesLabel -> {"x", "y", "z"},
AxesStyle -> {RGBColor[0, 0, 0], Opacity[1]}]]]


Defining the temperature dependent material properties (these will be interpolating functions loaded from a file):

(*material properties*)
k = 175; (* W/(m*K) *)
k1[temp_] := k*(1 + (210 - temp)*0.003)
Plot[k1[t], {t, 0, 300}]

G = 1.79*10^7; (* A/(V*m) *)
G1[temp_] := G*(1 + (210 - temp)*0.003)
Plot[G1[t], {t, 0, 300}]


Defining boundaries and boundary conditions (fixed temperatures at the ends, fixed voltage at the back and current density at the front):

(*boundary definitions*)
tol = 10^-6;
boundaryFront[x_, y_, z_] := Abs[z - 0] <= tol
boundaryBack[x_, y_, z_] := Abs[z - l] <= tol

(*Boundary Electro*)
bcElectroBack = DirichletCondition[V[x, y, z] == 0, boundaryBack[x, y, z]];
bcElectroFrontNeumann = NeumannValue[10^6, boundaryFront[x, y, z]];

(*Boundary heat*)
bcTempFront = DirichletCondition[T[x, y, z] == 200, boundaryFront[x, y, z]];
bcTempBack = DirichletCondition[T[x, y, z] == 220, boundaryBack[x, y, z]];


And here is where it fails, setting up the pde. Setting up the heat equation works fine, even with the temperature dependent conductivity. The voltage field equation with the constant material properties (with G) works, while the equation with the temperature dependent conductivity (G1) fails.

(*setting up coupled pde*)
QJouleCoupled = G1[T[x, y, z]]*Norm[Grad[V[x, y, z], {x, y, z}]]^2 /.Abs -> RealAbs;(*this works*)

eqHeatCoupled =
Inactive[Div][{{-k1[T[x, y, z]], 0, 0}, {0, -k1[T[x, y, z]], 0}, {0, 0, -k1[T[x, y, z]]}}.Inactive[Grad]
[T[x, y, z], {x, y, z}], {x, y, z}] - QJouleCoupled;

(*voltage field equation with temperature dependet conductivity does not work*)
eqElectroCoupled =
Inactive[Div][{{-G1[T[x, y, z]], 0, 0}, {0, -G1[T[x, y, z]], 0}, {0,0, -G1[T[x, y, z]]}}.Inactive[Grad]
[V[x, y, z], {x, y, z}], {x, y, z}];

(*voltage field equation with constant conductivity works*)
eqElectroCoupled =
Inactive[Div][{{-G, 0, 0}, {0, -G, 0}, {0, 0, -G}}.Inactive[Grad]
[V[x, y, z], {x, y, z}], {x, y, z}];

pdeCoupled = {{eqElectroCoupled == 0 + bcElectroFrontNeumann,eqHeatCoupled == 0},
{bcElectroBack, bcTempFront, bcTempBack}}


Solving the equations:

(*Solve coupled pde*)
measure = AbsoluteTiming[MaxMemoryUsed[{VCoupled, TCoupled} =
NDSolveValue[pdeCoupled,
{V, T},
{x, y, z} \[Element] mesh,
InitialSeeding -> {T[x, y, z] == 210, V[x, y, z] == 0.025},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {V -> 2, T -> 2}}}
]
]/(1024.^2)
];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]


Which throws the error when using the temperature dependent electr. conductivity. It says "PDE parsing error of ... Inconsistent equation dimensions" and "Lists ... are not of the same shape:

(*Visualize solution*)
VRange = MinMax[VCoupled["ValuesOnGrid"]];
legendBar = BarLegend[{"TemperatureMap", VRange}, LegendLabel -> Text[
Style["[V]",
Opacity[0.6]]]];
options = {AspectRatio -> Automatic, PerformanceGoal -> "Quality",
PlotPoints -> 50, Mesh -> None, PlotTheme -> "Detailed",
PlotLegends -> None, AxesLabel -> {x, y, z},
ColorFunctionScaling -> False, ImageSize -> Medium,
PlotLabel -> Style["Electric Potential Field: V(x,y,z)", 18],
Ticks -> {Automatic, Automatic, Automatic}};
Legended[RegionPlot3D[mesh,
ColorFunction ->
Function[{x, y, z},
ColorData[{"TemperatureMap", VRange}][VCoupled[x, y, z]]],
Evaluate[options]], legendBar]

TRange = MinMax[TCoupled["ValuesOnGrid"]];
legendBar =
BarLegend[{"TemperatureMap", TRange},
LegendLabel -> Text[Style["[K]", Opacity[0.6]]]];
options = {AspectRatio -> Automatic, PerformanceGoal -> "Quality",
PlotPoints -> 50, Mesh -> None, PlotTheme -> "Detailed",
PlotLegends -> None, AxesLabel -> {x, y, z},
ColorFunctionScaling -> False, ImageSize -> Medium,
PlotLabel -> Style["Temperature Field: T(x,y,z)", 18],
Ticks -> {Automatic, Automatic, Automatic}};
Legended[RegionPlot3D[mesh,
ColorFunction ->
Function[{x, y, z},
ColorData[{"TemperatureMap", TRange}][TCoupled[x, y, z]]],
Evaluate[options]], legendBar]


• Activate[eqElectroCoupled] gives you a hint. Namely Nonrectangular tensor encountered. You had a typo in in the last row of your matrix, {0,0-G1[...]} note the missing comma before -G1. Fixing this seems to work Commented Oct 18, 2021 at 16:35
• Actually the comma just went missing when creating the post. There still is a problem Commented Oct 18, 2021 at 16:55
• I am running Mathematica 12.1 by the way Commented Oct 18, 2021 at 16:55
• I only have access to a machine with 12.3 now, but copy-pasting the (now-corrected) code runs on 12.3. Perhaps this is related to this? Commented Oct 18, 2021 at 17:17
• Perhaps, I get the error reported in the link you send. Did you make sure you are using the temperature dependent definition of the electro equation eqElectroCoupled? I will try it on a newer version tomorrow. Commented Oct 18, 2021 at 18:05

• The OP's code runs on 12.3, so perhaps this is related to the parsing bug described here.

• Here's a slightly cleaned-up version to allow easier copy-pasting and verify indeed the correct definition of eqElectroCoupled was used

Needs["NDSolveFEM"];

xmin = 0; xmax = 0.1;
ymin = 0; ymax = 0.2;
l = 1;

cubi = Cuboid[{xmin, ymin, 0}, {xmax, ymax, l}];
mesh = ToElementMesh[cubi, MaxCellMeasure -> 0.00003];

k = 175;
k1[temp_] := k*(1 + (210 - temp)*0.003)

G = 1.79*10^7;
G1[temp_] := G*(1 + (210 - temp)*0.003)

bcElectroBack = DirichletCondition[V[x, y, z] == 0, z == l];
bcElectroFrontNeumann = NeumannValue[10^6, z == 0];

bcTempFront = DirichletCondition[T[x, y, z] == 200, z == 0];
bcTempBack = DirichletCondition[T[x, y, z] == 220, z == l];

QJouleCoupled =
G1[T[x, y, z]]*Norm[Inactive[Grad][V[x, y, z], {x, y, z}]]^2 /.
Abs -> RealAbs;

eqHeatCoupled =
Inactive[Div][{{-k1[T[x, y, z]], 0, 0}, {0, -k1[T[x, y, z]], 0}, {0,
0, -k1[T[x, y, z]]}} .
Inactive[Grad][T[x, y, z], {x, y, z}], {x, y, z}] - QJouleCoupled;

eqElectroCoupled =
Inactive[Div][{{-G1[T[x, y, z]], 0, 0}, {0, -G1[T[x, y, z]], 0}, {0,
0, -G1[T[x, y, z]]}} .
Inactive[Grad][V[x, y, z], {x, y, z}], {x, y, z}];

pdeCoupled = {{eqElectroCoupled == 0 + bcElectroFrontNeumann,
eqHeatCoupled == 0}, {bcElectroBack, bcTempFront, bcTempBack}};

measure =
AbsoluteTiming[
MaxMemoryUsed[{VCoupled, TCoupled} =
NDSolveValue[pdeCoupled, {V, T}, {x, y, z} \[Element] mesh,
InitialSeeding -> {T[x, y, z] == 210, V[x, y, z] == 0.025},
Method -> {"PDEDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {V -> 2, T -> 2}}}]]/(1024.^2)];


Note two things:

1. Used the simpler Dirichlet conditions z==0 and z==l which seem to work fine (does not affect solution)
2. Wrapped Grad in the definition of QJouleCoupled inside an Inactive (this does affect the Temperature field solution!)

Here are the resulting plots:

Without the use of Inactive the Temperature field solution is similar to OP's.

• Thank you very much! Commented Oct 18, 2021 at 18:45
• What is the benefit of wrapping the Grad inside an Inactive`? Commented Oct 19, 2021 at 6:59