# Help with 3D FEM calculation of a heat equation

I want to solve a heat transport problem in a long tube where 4 coolings rods are inserted. Fluid flows down axially, and there's radial heat conduction.

First, the shape is defined:

<< NDSolveFEM
shellID=0.6;
tubeLength=6;
rodOD=0.03;
rodCenters = shellID {{+0.25, -0.25}, {+0.25, +0.25}, {-0.25, -0.25}, {-0.25, +0.25}};
rodSection = RegionUnion @@ (Disk[#, rodOD] & /@ rodCenters);
tubeSection = RegionDifference[Disk[{0, 0}, shellID], rodSection];

mesh = ToElementMesh[tubeSection];


Then I set up the equation to be solved. Constants:

urcp=1500;
lamb=15;
u=200000;


Boundary conditions: the temperature is homogenous at the top, and the rods provide cooling.

init=DirichletCondition[t[z, x, y] == 240, ({x, y} ∈ mesh) && z == 0];

tubeBC = NeumannValue[(100 - t[z, x, y]), {z, x, y} ∈
RegionProduct[Line[{{0}, {tubeLength}}], rodSection]];


The convection/diffusion equation is:

equation = urcp*Derivative[1, 0, 0][t][z, x, y] == lamb (Derivative[0, 2, 0][t][z, x, y] + Derivative[0, 0, 2][t] t[z, x, y]) + u*tubeBC


FInally, the solution:

sol = NDSolveValue[{equation, init}, t, {z, 0, tubeLength}, {x, y} ∈ mesh]


This triggers error messages ('The ranges cannot be combined into a region'). How can I make this code work?

• On which border is the condition set NeumannValue[(100 - t[z, x, y])? – Alex Trounev Nov 15 '18 at 18:37
• You can use NeumannValue[(100 - t[z, x, y]), x^2 + y^2 < shellID^2] if you want the Robin condition on the rods only. – user21 Nov 16 '18 at 9:52
• The condition is set on the rods only. In this calculation the outer wall is adiabatic (no heat transfer). – Whelp Nov 16 '18 at 10:03
• According to the documentation, NeumannValue determines the coefficients of the flux based on the differential equation. How do we know for sure in this case that the flux is refered to the radial diffusion part lamb (Derivative[0, 2, 0][t][z, x, y] + Derivative[0, 0, 2][t][z, x, y]), and indeed does not take the convective part urcp*Derivative[1, 0, 0][t][z, x, y] into account ? – Whelp Nov 21 '18 at 9:53

Solution in the case when the Neumann condition is given at all boundaries where possible

shellID = 0.6;
tubeLength = 6;
rodOD = 0.03;
rodCenters =
shellID {{+0.25, -0.25}, {+0.25, +0.25}, {-0.25, -0.25}, {-0.25, \
+0.25}};
rodSection = RegionUnion @@ (Disk[#, rodOD] & /@ rodCenters);
tubeSection = RegionDifference[Disk[{0, 0}, shellID], rodSection];
urcp = 1500;
lamb = 15;
u = 200000;
init = t[0, x, y] == 240;

tubeBC = NeumannValue[(100 - t[z, x, y]), True];
equation =
urcp*Derivative[1, 0, 0][t][z, x, y] -
lamb (Derivative[0, 2, 0][t][z, x, y] +
Derivative[0, 0, 2][t][z, x, y]) == u*tubeBC;
sol = NDSolveValue[{equation, init},
t, {z, 0, tubeLength}, {x, y} \[Element] tubeSection]
Table[ContourPlot[sol[z, x, y], {x, -.6, .6}, {y, -.6, .6},
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> Row[{"z=", z}]], {z, 1, 6, 1}]


If we specify a different BC at the tube wall, for example, no heat transfer at the outer wall, then

tubeBC = NeumannValue[(100 - t[z, x, y]),
Norm[{x, y} - rodCenters[[1]]] == rodOD ||
Norm[{x, y} - rodCenters[[2]]] == rodOD ||
Norm[{x, y} - rodCenters[[3]]] == rodOD ||
Norm[{x, y} - rodCenters[[4]]] == rodOD];
bc1 = NeumannValue[0, x^2 + y^2 == shellID^2];
equation =
urcp*Derivative[1, 0, 0][t][z, x, y] -
lamb (Derivative[0, 2, 0][t][z, x, y] +
Derivative[0, 0, 2][t][z, x, y]) == u*tubeBC + bc1;
sol = NDSolveValue[{equation, init},
t, {z, 0, tubeLength}, {x, y} \[Element] tubeSection]
Table[ContourPlot[sol[z, x, y], {x, -.6, .6}, {y, -.6, .6},
ColorFunction -> "TemperatureMap", PlotLegends -> Automatic,
Contours -> 20, PlotLabel -> Row[{"z=", z}]], {z, 1, 6, 1}]


• Thanks. How can I specify a different BC at the tube wall, for example, no heat transfer at the outer wall? I've tried adding tubeBC2 = NeumannValue[ t[z, x, y], {z, x, y} [Element] RegionProduct[Line[{{0}, {tubeLength}}], Disk[{0, 0}, shellID]]] but without success. – Whelp Nov 16 '18 at 9:59
• So in this solution, no FEM is explicitely involved? I suppose NDSolve calls it internally nonetheless? – Whelp Nov 16 '18 at 10:05
• In version 11.3 which I run FEM is used automatically. I updated the code for adiabatic boundary conditions on the outer surface. – Alex Trounev Nov 16 '18 at 13:30