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:
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]
Activate[eqElectroCoupled]
gives you a hint. NamelyNonrectangular 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 $\endgroup$