# Heat Flow in a Pipe

I try to model the transient heat flow in a pipe, assuming that the temperature in radial direction doesn't change:

temperature u[t, \[CurlyPhi], z] , 0<z<10,0<\[CurlyPhi]<2Pi

The flux boundary conditions are described by NeumannValue. In direction of circumference two continuity conditions are considered:

Evaluating my model

NDSolveValue[{ 2.4 10^6   Derivative[1, 0, 0][u][t, \[CurlyPhi], z] == 172 (Derivative[0, 0, 2][u][t, \[CurlyPhi], z] +30.5 Derivative[0, 2, 0][u][t, \[CurlyPhi], z])
+ 1/138 1.85 10^7 NeumannValue[1, (-10 \[Degree] <= \[CurlyPhi]<= 10 \[Degree]) &&z == 10    ]
+ NeumannValue[0, (- 10 \[Degree]  <= \[CurlyPhi] <= 10 \[Degree]) &&z ==  0  ]
, u[0, \[CurlyPhi], z] == 0
, u[t, 0, z] == u[t, 2 Pi, z]
, Derivative[0, 1, 0][u][t, 0, z] ==Derivative[0, 1, 0][u][t, 2 Pi, z]},
u,{t, 0, 1}, {\[CurlyPhi], 0, 2 Pi}, {z, 0, 10}
, Method -> {"MethodOfLines", "TemporalVariable" -> t,"SpatialDiscretization" -> {"TensorProductGrid" }}]

(*NDSolveValue::bcart: Warning: an insufficient number of
boundary conditions have been specified for the direction of
independent variable z. Artificial boundary effects may be present in the solution.*)

gives an error message "insufficient boundary conditions in z-direction", inspite of two NeumannValue's, which I don't understand

What's wrong with my code? Thanks!

addendum #1(21.11.2019) (thanks for the several contributions! )

Substituting NeumannValue by explicit Derivative-boundaries the modified simulation

NDSolveValue[{ 2.4 10^6 (Derivative[1, 0, 0][u][t, \[CurlyPhi], z] ) ==172 (Derivative[0, 0, 2][u][t, \[CurlyPhi], z] + 30.5 Derivative[0, 2, 0][u][t, \[CurlyPhi], z])
, u[0, \[CurlyPhi], z] == 0
, u[t, 0, z] == u[t, 2 Pi, z]
, Derivative[0, 1, 0][u][t, 0, z] ==Derivative[0, 1, 0][u][t, 2 Pi, z]
, Derivative[0, 0, 1][u][t, \[CurlyPhi], 10] == If[-10 \[Degree] <= \[CurlyPhi]<= 10 \[Degree],1/138 1.85 10^7 , 0]
, Derivative[0, 0, 1][u][t, \[CurlyPhi], 0] == 0},
u,{t, 0, 1}, {\[CurlyPhi], 0, 2 Pi}, {z, 0, 10}, Method -> {"MethodOfLines", "TemporalVariable" -> t,"SpatialDiscretization" -> {"TensorProductGrid" }}

runs without error but only evaluates the solution u==0, that means the flux boundary at z==0 doesn't contribute to the solution???

Learning&Knowing that flux boundaries might be more easier formulated in FEM I tried to find a executable FEM-model:

NDSolveValue[{2.4 10^6 Derivative[1, 0, 0][u][t, \[CurlyPhi], z] ==172 (Derivative[0, 0, 2][u][t, \[CurlyPhi], z] + 30.5 Derivative[0, 2, 0][u][t, \[CurlyPhi], z]) +
1/138 1.85 10^7 NeumannValue[1, (-10 \[Degree] <= \[CurlyPhi] <= 10 \[Degree]) && z ==  10] +
NeumannValue[0, (-10 \[Degree] <= \[CurlyPhi] <= 10 \[Degree]) && z ==  0], u[0, \[CurlyPhi], z] == 0
, PeriodicBoundaryCondition[u[t, \[CurlyPhi], z], \[CurlyPhi] ==  \[Pi],Function[\[Phi], \[Phi] - 2 \[Pi]]]},
u, {t, 0, 1}, {\[CurlyPhi], -Pi, Pi}, {z, 0, 10},
Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement"}}]

But Mathematica gives error message NDSolveValue::bcnop: No places were found on the boundary where -10 \[Degree]<=\[CurlyPhi]<=10 \[Degree]&&z==0 was True, so NeumannValue[0,-10 \[Degree]<=\[CurlyPhi]<=10 \[Degree]&&z==0] will effectively be ignored.

What's wrong here? Thanks!

• I think the periodic boundary conditions should be formulated with, well, PeriodicBoundaryCondition. E.g. PeriodicBoundaryCondition[u[t, \[Phi], z], \[Phi] == 2 \[Pi], Function[\[Phi], \[Phi] - 2 \[Pi]]] and PeriodicBoundaryCondition[ Derivative[0, 1, 0][u][t, \[Phi], z], \[Phi] == 2 \[Pi], Function[\[Phi], \[Phi] - 2 \[Pi]]]. However, Mathematica complains about the latter. =| – Henrik Schumacher Nov 19 '19 at 21:17
• @HenrikSchumacher Thanks, it's dodgy... – Ulrich Neumann Nov 19 '19 at 21:40
• @Henrik Schumacher, what leads you to belive that? Maybe I need to improve the documentation? – user21 Nov 20 '19 at 4:41
• One obvious mistake is, you're using TensorProductGrid together with NeumannValue, but NeumannValue is only for FiniteElement, at least now. – xzczd Nov 20 '19 at 5:26
• @xzczd Thanks. In the documentation of NeumannValue I didn't find that restriction. It would imply that no flux -conditions could be imposed with TensorPorductGridand difference-methods? – Ulrich Neumann Nov 20 '19 at 7:32

There are certain triggers that force the use of the finite element method. These are listed here and I am working on adding this to the FEM documentation in product. So the use of NeumannValue is going to force the the FEM - but that contradicts with the TensorProductGrid method you requested. So either use the finite element method (then you need to get rid of the periodic derivatives - see this ) or get rid of the NeumannValue to use the "TensorProductGrid" (TPG) method, which I think you want to do here. To convert the NeuammValue to a Derivative it is mandatory to understand their relation which is discussed here. One thing that could be done though is to make NDSolve give a message in that case.

You could also solve this on an Annulus (if I understand correctly what you want to do) but in either case you'd need something that drives the temperature. As a template you could use:

fun = NDSolveValue[{Derivative[1, 0, 0][u][t, x, y] ==
10 (Laplacian[u[t, x, y], {x, y}]) +
NeumannValue[
1, (x^+y^2 >= 2) && (x > 1 && (-1 >= y && y <= 1))],
u[0, x, y] == 0
}, u, {t, 0, 1}, Element[{x, y}, Annulus[{0, 0}, {1, 2}]],
Method -> {"MethodOfLines", "TemporalVariable" -> t}];
ContourPlot[fun[1, x, y], {x, y} \[Element] fun["ElementMesh"]]

• Thanks for your hints, what is "TPG method"? My understanding of NeumannValue is the normal derivative of the unknown function on the boundary surface. In the documentation NeumannValue special type of pde is shown, somewhat confusing for me... – Ulrich Neumann Nov 20 '19 at 9:18
• @UlrichNeumann, TPG = TensorProductGrid method. Concerning the NeumannValue, that's not quite correct, though in some cases it can be used for that. The NeumannValues is tied to the c values in Div[c Grad[u[x,y],{x,y}],{x,y}] – user21 Nov 20 '19 at 9:47
• Thank you for your comment and the template! The purpose of my question is as follows: I wrote a huge simulation model (3D hollow cylinder, cartesian coordinates , transient,>10000elements) which works quite well, but gives error message peclet-number. Increasing element number doesn't succeed. That's why I try to reduce this modell. – Ulrich Neumann Nov 20 '19 at 9:58
• I tried to substitute the NeumannValue in your example by DirichletCondition[Grad[u[t, x, y], {x, y}].{x, y} ==1, (x^2 + y^2 >= 2) && (x > 1 && (-1 >= y && y <= 1 ))] but TensorProductGriddoesn't work as expected. Any idea how to make it work? Thanks! – Ulrich Neumann Nov 21 '19 at 9:22
• @UlrichNeumann, Neither, DirichletCondition nor NeumannValue nor PeriodicBoundaryCondition can be used with the TensorProductGrid method and more in the finite element method you can not use DirichletCondition[Grad... either. To make this work with the TPG method you need to eliminate DirichletCondtion, NeumannValue and PeriodicBoundaryCondition from your code. Maybe you could have a second post with your actual 3D example? – user21 Nov 21 '19 at 9:45

Your problem and @user21's answer suggests that you have half symmetry at $$\pi$$. You could use FEMAddOns to build a FEM mesh with a heated boundary $$0-5^\circ$$ and solve on the domain $$0-180^\circ$$ like so:

Needs["NDSolveFEM"]
Needs["FEMAddOns"]
ht = 10;
len = 5 \[Degree];
top = ht;
bot = 0;
left = 0;
right = len;
left2 = right;
right2 = 180 \[Degree];
bounds = <|wall -> 1, hot -> 2|>;
bmeshheated =
ToBoundaryMesh[
"Coordinates" -> {{left, bot}(*1*), {right, bot}(*2*), {right,
top}(*3*), {left, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*), {3, 4}(*4*)}, {1, 1, 1,
bounds[hot]}]}];
bmeshheated[
"Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]];

bmeshinsulated =
ToBoundaryMesh[
"Coordinates" -> {{left2, bot}(*1*), {right2, bot}(*2*), {right2,
top}(*3*), {left2, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*), {3, 4}(*4*)}, {1, 1, 1,
1}]}];
bmeshinsulated[
"Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]];
bmesh = BoundaryElementMeshJoin[bmeshheated, bmeshinsulated];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]];
bmesh["Wireframe"];
mesh = ToElementMesh[bmesh]
mesh["Wireframe"]

Once you have an well defined FEM mesh, it is relative straightforward to pose the FEM problem using coefficient form.

rhocp = 2.4 10^6;
k = 172.;
r = Sqrt[1/30.5];
q = 1/138 1.85 10^7/2.;
nv = NeumannValue[q, ElementMarker == bounds[hot]];
op = rhocp \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[t, x, y]$$\) -
Inactive[
u[t, x, y], {x, y}]), {x, y}];
uif = NDSolveValue[{op == nv , u[0, x, y] == 0},
u, {t, 0, 1}, {x, y} \[Element] mesh];
imgs = Plot3D[
If[x < 0, uif[#, -x, z], uif[#, x, z]], {x, -Pi, Pi}, {z, 9.5,
10}, PlotRange -> {-0.1, 1.4}, ColorFunction -> "DarkBands"] & /@
Subdivide[0, 1, 60];
ListAnimate@imgs

The mesh will need some refinement, but if you have the interpolation function, then you should be able to transform it to cylindrical coordinates.

# Update To Include Scaling and Elevated Surface

The spatial and temporal domain was very large relative to where all the action was occurring. In these cases, non-dimensionalizing the PDE will help size the domains relative to characteristic time and length scales. In this case, I used the concept of a penetration depth in infinite plate heating. The penetration depth is given by:

$${y_c} = 4\sqrt {\alpha t}$$

After non-dimensionalizing the equations, we can create a domain where the flux will penetrate to the symmetry condition and to the insulate wall condition at $$z=0$$. Here is the code to create a mesh on the new domain with refinement near the heater:

r = 1;
ht = 4 r;
len = 5 \[Degree];
top = ht;
bot = 0;
left = 0;
right = len;
left2 = right;
right2 = 180 \[Degree];
bounds = <|wall -> 1, hot -> 2|>;
bmeshheated =
ToBoundaryMesh[
"Coordinates" -> {{left, bot}(*1*), {right, bot}(*2*), {right,
top}(*3*), {left, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*), {3, 4}(*4*)}, {1, 1, 1,
bounds[hot]}]}];
bmeshheated[
"Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]];

bmeshinsulated =
ToBoundaryMesh[
"Coordinates" -> {{left2, bot}(*1*), {right2, bot}(*2*), {right2,
top}(*3*), {left2, top}(*4*)},
"BoundaryElements" -> {LineElement[{{1, 2}(*bottom edge*)(*1*), {4,
1}(*left edge*)(*2*), {2, 3}(*3*), {3, 4}(*4*)}, {1, 1, 1,
1}]}];
bmeshinsulated[
"Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]];
bmesh = BoundaryElementMeshJoin[bmeshheated, bmeshinsulated];
bmesh["Wireframe"["MeshElementMarkerStyle" -> Blue,
"MeshElementStyle" -> {Green, Red, Black}, ImageSize -> Large]]
mesh = ToElementMesh[bmesh,
MeshRefinementFunction ->
Function[{vertices, area},
area > 0.0002 (1 + 10 Norm[Mean[vertices] - {0, ht}])]];
mesh["Wireframe"]

## Element Mesh Showing Refinement

We will solve the non-dimensional PDE and visualize the results with an elevated surface using ParametricPlot3D.

ct = CoordinateTransformData["Cylindrical" -> "Cartesian", "Mapping"];
rhocp = 16/Pi^2;
alpha = 1;
q = 1;
nv = NeumannValue[q, ElementMarker == bounds[hot]];
op = rhocp \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]$$u[t, x, y]$$\) -
Inactive[
u[t, x, y], {x, y}]), {x, y}];
uif = NDSolveValue[{op == nv , u[0, x, y] == 0},
u, {t, 0, 10}, {x, y} \[Element] mesh];
ppfn =  ParametricPlot3D[{ct[{r + 10 uif[#, x, z], x, z}],
ct[{r + 10 uif[#, x, z], -x, z}]}, {x, z} \[Element] mesh,
ColorFunction ->
Function[{x, y, z, u, v},
ColorData["SunsetColors"][10 uif[#, u, v]]],
ColorFunctionScaling -> False,
MeshFunctions -> {Function[{x, y, z}, Norm@{x, y}]}, Mesh -> 25,
PlotRange -> {{-2, 4}, {-2, 2}, {0, ht}}, Boxed -> False,
Axes -> False,
ViewPoint -> {1.6822406480044492, -1.3581436516840697, \
-2.6029814105352025},
ViewVertical -> {0.9575175915664486,
0.060531137127969, -0.2819504269881192},
Background -> GrayLevel[.1]] &;
imgs = ppfn[#] & /@ Subdivide[0, 10, 80];
ListAnimate@imgs

The results seem reasonable. The temperature is highest at the flux boundary and the pipe begins to "swell" at the zero flux conditions.

# Time Varying Heater

Something a little more interesting/mesmerizing would be to add a time varying heat flux and watch the waves propagate down the pipe, which is easily done by replacing two lines and re-running the simulation.

q[t_] = 1 + Sin[Pi*t] + Cos[Pi*t] + Sin[2*Pi*t] + Cos[2*Pi*t];
nv = NeumannValue[q[t], ElementMarker == bounds[hot]];

• Thanks for your help. Symmetry is only fullfilled in my simplified example. Additionally I need the flux boundaries. – Ulrich Neumann Nov 20 '19 at 13:46
• The default boundary condition is zero flux, so a flux BC has been specified on all boundaries. – Tim Laska Nov 20 '19 at 14:17
• No, the flux boundary isn't zero at z==10! – Ulrich Neumann Nov 21 '19 at 8:04
• Thanks for the impressing simulation! – Ulrich Neumann Nov 21 '19 at 8:05
• You are welcome. I updated the post to show the heated edge boundary marker (2) for clarity. It was set to be non-zero by nv = NeumannValue[q, ElementMarker == bounds[hot]]. All other boundaries are set to zero and internal boundaries are ignored. – Tim Laska Nov 21 '19 at 12:45

I'll fix the code using TensorProductGrid in #1. There're at least 5 issues:

1. Considering your b.c., the domain for $$φ$$ should probably be {φ, -Pi, Pi}.

2. The default grid is too coase (15 points in each dimension) to capture the narrow flux at $$z=0$$.

3. The i.c. and b.c. are inconsistent, which leads to the problem mentioned here.

4. Setting a non-zero ScaleFactor as suggested above triggers a bug mentioned here.

5. The method for ODE solving needs to be adjusted a bit, or NDSolve will automatically choose a rather slow method when the grid is dense enough. (It's unclear what the method is. )

After resolving all these, we obtain:

test = NDSolveValue[
With[{u = u[t, φ, z]},
{2.4 10^6 D[u, t] == 172 (D[u, z, z] + 30.5 D[u, φ, φ]),
u == 0 /. t -> 0,
(u /. φ -> -Pi) == (u /. φ -> Pi),
D[u, z] ==
SimplifyPWToUnitStep@
PiecewiseExpand@
If[-10 ° <= φ <= 10 °, 1/138 1.85 10^7, 0] /. z -> 10,
D[u, z] == 0 /. z -> 0}], u, {t, 0, 1}, {φ, -Pi, Pi}, {z, 0, 10},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 10},
"SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 810,
"MinPoints" -> 810, "DifferenceOrder" -> 4}, Method -> Adams}]; // AbsoluteTiming
(* ibcinc warning *)
(* {86.9492, Null} *)

Plot[test[1, phi, 10], {phi, -Pi, Pi}, PlotRange -> All]

• Thank you very much, I have to elaborate some time. Testing your code I increased the simulation time, but the plot doesn't change, though there is no heat conduction in the pipe??? – Ulrich Neumann Nov 21 '19 at 10:21
• @UlrichNeumann I guess you forgot to change the time inside Plot? – xzczd Nov 21 '19 at 11:37
• Unfortunately not, I tried t=3 (simulation and plot). – Ulrich Neumann Nov 21 '19 at 11:45
• @UlrichNeumann Please show me your code. – xzczd Nov 21 '19 at 11:47
• Here is my code Show[{Plot[test[3, phi, 10], {phi, -Pi, Pi}, PlotRange ->All,GridLines -> {{-10 \[Degree], 10 \[Degree]}, None}],Plot[test[1, phi, 10], {phi, -Pi, Pi}, PlotRange -> All,GridLines -> {{-10 \[Degree], 10 \[Degree]}, None}]}]. It shows a small conduction effect! – Ulrich Neumann Nov 21 '19 at 11:50