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

Question

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["NDSolve`FEM`"];

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

Answer 1

Posting an answer to consolidate some of the comments above.

• 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["NDSolve`FEM`"];

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.