Heat Flow in a Pipe

Question

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???

addendum #2(21.11.2019)

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!

Answer 1

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"]]

Answer 2

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["NDSolve`FEM`"]
ResourceFunction["FEMAddOnsInstall"][]
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[
     Div][(-{{k/r^2, 0}, {0, k}}.Inactive[Grad][
        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"]

Boundary Mesh with Heated Edge Marked

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[
     Div][(-{{alpha, 0}, {0, alpha}}.Inactive[Grad][
        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]];

Answer 3

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] == 
        Simplify`PWToUnitStep@
         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]