7
$\begingroup$

I am working on a multiphysics problem involving heat transfer and electrostatics. I have been messing around with the Joule Heating Tutorial Case and got stuck with the heat source term not working. I am running Mathematica 12.1 on Windows 10 and can not update Mathematica as I will add to some existing code. The tutorial (note, that the link above refers to a newer version of Mathematica, I used the documentation under "Help") is one-way-coupled, as the electric conductivitiy does not depend on the temperature, therefore the voltage field is solved first and then used to calculate the source term for the temperature field.

The source term Q in the tutorial is defined as the following:

...
Vfun = NDSolveValue[pde, V, {x, y, z} ∈ mesh];
...
Q = G * Norm[Grad[Vfun[x, y, z], {x, y, z}]]^2;
...

Now I want to solve the voltage and temperature field coupled (to use temperature depedent material properties) and I run into problems defining the source term Q as follows:

QCoupled = G *Norm[Grad[V[x, y, z], {x, y, z}]]^2;

So I set up my pde:

(*Set up coupled pde*)
QCoupled = Subscript[G, tungsten]*Norm[Grad[V[x, y, z], {x, y, z}]]^2;
pdeElectroCoupled = ElectrostaticsModel[V[x, y, z], {x, y, z}, G];
pdeHeatCoupled = 
HeatTransferModel[T[x, y, z], {x, y, z}, 
  k + T[x, y, z] * 10^(-10), ρ, Cp, "NoFlow", QCoupled]
pdeCoupled = {pdeElectroCoupled == 0, 
    pdeHeatCoupled == 
     Subscript[Γ, radiation] + 
      Subscript[Γ, convective], 
    Subscript[Γ, volt], Subscript[Γ, 
    temperature]} /. parameters;

And try running NDSolveValue:

(*Solve coupled stuff*)
measure = 
      AbsoluteTiming[
       MaxMemoryUsed[{VCoupled, TCoupled} = 
          NDSolveValue[
           pdeCoupled, {V, T}, {x, y, z} ∈ mesh]]/(1024.^2)];

It throws errors:

"... does not evaluate to a numeric matrix of dimensions {1,3} at the
coordinate ..."

"The linearization process in PDESolve failed."

So my question is, how do I define the source term for the coupled system?

Find my whole code below, if it helps you helping me.

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

(*Geometry*)
L = 0.1;
r = 0.0025;
R = 0.0075;
arc1 = Table[
   R*{Cos[θ], Sin[θ], 0.}, {θ, -π/2, π/
     2, π/30}];
arc2 = Reverse[-arc1];
arc3 = Table[
   R*{Cos[θ], Sin[θ], 0.}, {θ, π/2, 
    0, -π/30}];
pathPts = 
  Join[{{0.`, 0.`, 0.`}, {L, 0.`, 0.`}}, ({L, R, 0.`} + #1 &) /@ 
    arc1, {{L, 2 R, 0.`}, {R + r, 2 R, 
     0.`}}, ({R + r, 3 R, 0.`} + #1 &) /@ 
    arc2, {{R + r, 4 R, 0.`}, {L + R + r, 4 R, 
     0.`}}, ({L + R + r, 3 R, 0.`} + #1 &) /@ 
    arc3, {{L + 2 R + r, 3 R, 0.`}, {L + 2 R + r, -r, 0.`}}];
Graphics3D[{Thick, Blue, Line[pathPts]}, Sequence[
  PlotRange -> {All, All, {(-2) r, 2 r}}, Axes -> True, 
   AxesLabel -> {x, y, z}, PlotLabel -> "Centerline of the wire", 
   Ticks -> {Automatic, {0, 0.03}, {-0.005, 0.005}}]];
wire = Graphics3D[{CapForm["Butt"], Tube[pathPts, r]}];
bmeshRegion = 
  BoundaryDiscretizeGraphics[wire, MaxCellMeasure -> 2.5*10^-10];
mesh = ToElementMesh[bmeshRegion]

(*Physics models*)
ClearAll[ElectrostaticsModel]
ElectrostaticsModel[V_, X_List, G_] := Module[{a, factor},
  factor = -G*IdentityMatrix[Length[X]];
  a = PiecewiseExpand[Piecewise[{{factor, True}}]];
  Inactive[Div][a.Inactive[Grad][V, X], X]] 

ClearAll[HeatTransferModel]
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]

(*Boundaries*)
tolerance = 0.01 r;
leftEndQ[x_, y_, z_] := Abs[x] <= tolerance && y^2 + z^2 <= r^2
rightEndQ[x_, y_, z_] := 
 Abs[y + r] <= tolerance && (x - (L + 2 R + r))^2 + z^2 <= r^2

(*Parameters*)
Subscript[ρ, tungsten] = 1.93*10^4;
Subscript[Cp, tungsten] = 134;
Subscript[k, tungsten] = 175;
Subscript[G, tungsten] = 1.79*10^7;
parameters = {G -> Subscript[G, tungsten], ρ -> Subscript[ρ,
     tungsten], Cp -> Subscript[Cp, tungsten], 
   k -> Subscript[k, tungsten]};

(*Boundary electro*)
V0 = 0.2;
Subscript[Γ, volt] =
  {DirichletCondition[V[x, y, z] == V0, leftEndQ[x, y, z]],
   DirichletCondition[V[x, y, z] == 0, rightEndQ[x, y, z]]};

(*Boundary heat*)
Subscript[T, amb] = 300;
Subscript[Γ, temperature] = 
  DirichletCondition[T[x, y, z] == Subscript[T, amb], 
   leftEndQ[x, y, z] || rightEndQ[x, y, z]];

σ = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
ε = 1;
Subscript[Γ, radiation] = 
  NeumannValue[ε*σ*(Subscript[T, amb]^4 - 
      T[x, y, z]^4), Not[leftEndQ[x, y, z] || rightEndQ[x, y, z]]];

h = 20;
Subscript[Γ, convective] = 
  NeumannValue[h*(Subscript[T, amb] - T[x, y, z]), 
   Not[leftEndQ[x, y, z] || rightEndQ[x, y, z]]];

(*Set up one way coupled pde*)
QUncoupled = 
  Subscript[G, tungsten]*Norm[Grad[Vuncoupled[x, y, z], {x, y, z}]]^2;

ElectrostaticsUncoupled = 
  ElectrostaticsModel[V[x, y, z], {x, y, z}, G];
HeatUncoupled = 
  HeatTransferModel[T[x, y, z], {x, y, z}, 
   k + T[x, y, z] * 10^(-10), ρ, Cp, "NoFlow", QUncoupled];

pdeElectroUncoupled = {ElectrostaticsUncoupled == 0, 
    Subscript[Γ, volt]} /. parameters;
pdeHeatUncoupled = {HeatUncoupled == 
    Subscript[Γ, radiation] + Subscript[Γ,
      convective], Subscript[Γ, temperature]} /. 
  parameters
(*Solve one way coupled stuff*)
Vuncoupled = 
  NDSolveValue[pdeElectroUncoupled, V, {x, y, z} ∈ mesh];

measure = 
  AbsoluteTiming[
   MaxMemoryUsed[
     Tuncoupled = 
      NDSolveValue[pdeHeatUncoupled, 
       T, {x, y, z} ∈ mesh]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]

(*Set up coupled pde*)
QCoupled = Subscript[G, tungsten]*Norm[Grad[V[x, y, z], {x, y, z}]]^2;
Qtest = Subscript[G, tungsten] * V[x, y, z]/ (10^10) + 
  1000000; (*this one works just fine, though making no sense*)

pdeElectroCoupled = ElectrostaticsModel[V[x, y, z], {x, y, z}, G];
pdeHeatCoupled = 
 HeatTransferModel[T[x, y, z], {x, y, z}, 
  k + T[x, y, z] * 10^(-10), ρ, Cp, "NoFlow", QCoupled]

pdeCoupled = {pdeElectroCoupled == 0, 
    pdeHeatCoupled == 
     Subscript[Γ, radiation] + 
      Subscript[Γ, convective], 
    Subscript[Γ, volt], Subscript[Γ, 
    temperature]} /. parameters;

(*Solve coupled stuff*)
measure = 
  AbsoluteTiming[
   MaxMemoryUsed[{VCoupled, TCoupled} = 
      NDSolveValue[
       pdeCoupled, {V, T}, {x, y, z} ∈ mesh]]/(1024.^2)];
Print["Time -> ", measure[[1]], "\nMemory -> ", measure[[2]]]

(*Result visualization*)
VRange = MinMax[VCoupled["ValuesOnGrid"]];
legendBar = BarLegend[{"TemperatureMap", VRange}, LegendLabel -> Text[
Style["[V]", 
Opacity[0.6]]]];
options = {
Sequence[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, {-r, r}}]};
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 = {
Sequence[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, {-r, r}}]};
Legended[RegionPlot3D[mesh, 
  ColorFunction -> 
   Function[{x, y, z}, 
    ColorData[{"TemperatureMap", TRange}][TCoupled[x, y, z]]], 
  Evaluate[options]], legendBar]
Improve this question
$\endgroup$

2 Answers 2

Reset to default
7
$\begingroup$

If you look carefully at the message you see that there is an Abs'[]. It could not take the derivative of Abs.

Replace the Abs with RealAbs to fix this issue:

(*Set up coupled pde*)
QCoupled = 
 Subscript[G, tungsten]*Norm[Grad[V[x, y, z], {x, y, z}]]^2 /. 
  Abs -> RealAbs
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much! To be honest, I don't really get the difference between the two from the documentation. But the documentation says "RealAbs is continuous and differentiable everywhere except at the origin". Is there a possibility of running into errors due to this? $\endgroup$
    – Tobias
    Jul 30, 2021 at 9:28
  • 1
    $\begingroup$ @Tobias, you might run into problems when you solve PDEs in the complex plane. From the type questions you have asked here, it's noting to worry about. $\endgroup$
    – user21
    Jul 30, 2021 at 9:30
6
$\begingroup$

Simplified version of the code working fine

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]

(*Geometry*)
L = 0.1;
r = 0.0025;
R = 0.0075;
arc1 = Table[
   R*{Cos[\[Theta]], Sin[\[Theta]], 0.}, {\[Theta], -\[Pi]/2, \[Pi]/
     2, \[Pi]/30}];
arc2 = Reverse[-arc1];
arc3 = Table[
   R*{Cos[\[Theta]], Sin[\[Theta]], 0.}, {\[Theta], \[Pi]/2, 
    0, -\[Pi]/30}];
pathPts = 
  Join[{{0.`, 0.`, 0.`}, {L, 0.`, 0.`}}, ({L, R, 0.`} + #1 &) /@ 
    arc1, {{L, 2 R, 0.`}, {R + r, 2 R, 
     0.`}}, ({R + r, 3 R, 0.`} + #1 &) /@ 
    arc2, {{R + r, 4 R, 0.`}, {L + R + r, 4 R, 
     0.`}}, ({L + R + r, 3 R, 0.`} + #1 &) /@ 
    arc3, {{L + 2 R + r, 3 R, 0.`}, {L + 2 R + r, -r, 0.`}}];
Graphics3D[{Thick, Blue, Line[pathPts]}, 
  Sequence[PlotRange -> {All, All, {(-2) r, 2 r}}, Axes -> True, 
   AxesLabel -> {x, y, z}, PlotLabel -> "Centerline of the wire", 
   Ticks -> {Automatic, {0, 0.03}, {-0.005, 0.005}}]];
wire = Graphics3D[{CapForm["Butt"], Tube[pathPts, r]}];
bmeshRegion = 
  BoundaryDiscretizeGraphics[wire, MaxCellMeasure -> 2.5*10^-10];
mesh = ToElementMesh[bmeshRegion]
(*Physics models*)
ClearAll[ElectrostaticsModel]
ElectrostaticsModel[V_, X_List, G_] := 
 Module[{a, factor}, factor = -G*IdentityMatrix[Length[X]];
  a = PiecewiseExpand[Piecewise[{{factor, True}}]];
  Inactive[Div][a . Inactive[Grad][V, X], X]]

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]

(*Boundaries*)
tolerance = 0.01 r;
leftEndQ[x_, y_, z_] := Abs[x] <= tolerance && y^2 + z^2 <= r^2
rightEndQ[x_, y_, z_] := 
 Abs[y + r] <= tolerance && (x - (L + 2 R + r))^2 + z^2 <= r^2

(*Parameters*)
Subscript[\[Rho], tungsten] = 1.93*10^4;
Subscript[Cp, tungsten] = 134;
Subscript[k, tungsten] = 175;
Subscript[G, tungsten] = 1.79*10^7;
parameters = {G -> Subscript[G, tungsten], \[Rho] -> 
    Subscript[\[Rho], tungsten], Cp -> Subscript[Cp, tungsten], 
   k -> Subscript[k, tungsten]};

(*Boundary electro*)
V0 = 0.2;
Subscript[\[CapitalGamma], 
   volt] = {DirichletCondition[V[x, y, z] == V0, leftEndQ[x, y, z]], 
   DirichletCondition[V[x, y, z] == 0, rightEndQ[x, y, z]]};

(*Boundary heat*)
Subscript[T, amb] = 300;
Subscript[\[CapitalGamma], temperature] = 
  DirichletCondition[T[x, y, z] == Subscript[T, amb], 
   leftEndQ[x, y, z] || rightEndQ[x, y, z]];

\[Sigma] = First[UnitConvert[Quantity["StefanBoltzmannConstant"]]];
\[CurlyEpsilon] = 1;
Subscript[\[CapitalGamma], radiation] = 
  NeumannValue[\[CurlyEpsilon]*\[Sigma]*(Subscript[T, amb]^4 - 
      T[x, y, z]^4), True];

h = 20;
Subscript[\[CapitalGamma], convective] = 
  NeumannValue[h*(Subscript[T, amb] - T[x, y, z]), True];

(*Set up one way coupled pde*)
QUncoupled = 
  Subscript[G, tungsten]*Norm[Grad[Vuncoupled[x, y, z], {x, y, z}]]^2;

ElectrostaticsUncoupled = 
  ElectrostaticsModel[V[x, y, z], {x, y, z}, G];
HeatUncoupled = 
  HeatTransferModel[T[x, y, z], {x, y, z}, 
   k + T[x, y, z]*10^(-10), \[Rho], Cp, "NoFlow", QUncoupled];

pdeElectroUncoupled = {ElectrostaticsUncoupled == 0, 
    Subscript[\[CapitalGamma], volt]} /. parameters;
pdeHeatUncoupled = {HeatUncoupled == 
    Subscript[\[CapitalGamma], radiation] + 
     Subscript[\[CapitalGamma], convective], 
   Subscript[\[CapitalGamma], temperature]} /. parameters
(*Set up coupled pde*)
QCoupled = 
  Subscript[G, 
    tungsten]*(Grad[V[x, y, z], {x, y, z}] . 
     Grad[V[x, y, z], {x, y, z}]);
Qtest = Subscript[G, tungsten]*
  V[x, y, z]/(10^10);(*this one works just fine,though making no \
sense*)pdeElectroCoupled = 
 ElectrostaticsModel[V[x, y, z], {x, y, z}, G];
pdeHeatCoupled = 
 HeatTransferModel[T[x, y, z], {x, y, z}, 
  k + T[x, y, z]*10^(-10), \[Rho], Cp, "NoFlow", QCoupled]

pdeCoupled = {pdeElectroCoupled == 0, 
    pdeHeatCoupled == 
     Subscript[\[CapitalGamma], radiation] + 
      Subscript[\[CapitalGamma], convective], 
    Subscript[\[CapitalGamma], volt], 
    Subscript[\[CapitalGamma], temperature]} /. parameters;

(*Solve coupled stuff*)
{VCoupled, TCoupled} = 
  NDSolveValue[pdeCoupled, {V, T}, {x, y, z} \[Element] mesh];

Result visualization

VRange = MinMax[VCoupled["ValuesOnGrid"]];
legendBar = 
  BarLegend[{"TemperatureMap", VRange}, 
   LegendLabel -> Text[Style["[V]", Opacity[0.6]]]];
options = {Sequence[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, {-r, r}}]};
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 = {Sequence[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, {-r, r}}]};
Legended[RegionPlot3D[mesh, 
  ColorFunction -> 
   Function[{x, y, z}, 
    ColorData[{"TemperatureMap", TRange}][TCoupled[x, y, z]]], 
  Evaluate[options]], legendBar]

Figure 1

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much! To sum up your solution: substitute Norm[Grad[V[x, y, z], {x, y, z}]]^2 with (Grad[V[x, y, z], {x, y, z}] . Grad[V[x, y, z], {x, y, z}]) $\endgroup$
    – Tobias
    Jul 30, 2021 at 9:31
  • 1
    $\begingroup$ Yes, it is. Also we can use trick with RealAbs, as user21 recommended. $\endgroup$
    – Alex Trounev
    Jul 30, 2021 at 10:48

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.