Efficient generation of a 3D mesh to solve a diffusion PDE

Question

I'm new to FEM and trying to figure out what is the most efficient way of generating a mesh and solving an equation. As a warm-up to the particular problem I'm interested in, I want to calculate the evolution of the temperature profile in a gas flow within a reactor shell with two cooling tubes maintained at a given temperature.

i defined the region like this :

<< NDSolve`FEM` 
{reactorLength,shellID, tubeOD}={6.0, 0.5, 0.03};
tubeRegion = 
 RegionUnion[
 Cylinder[{{0, -0.25*shellID, 
 0}, {0, -0.25*shellID, +reactorLength}}, 1 tubeOD], 
 Cylinder[{{0, 0.25*shellID, 0}, {0, 0.25*shellID, +reactorLength}}, 
 1 tubeOD]];
tubeSection = 
RegionUnion[Disk[{0, -0.25*shellID}, 1 tubeOD], 
Disk[{0, +0.25*shellID}, 1 tubeOD]];
region = RegionDifference[
Cylinder[{{0, 0, 0}, {0, 0, +reactorLength}}, shellID], tubeRegion];

So we have a large shell with two tubes in the middle. Next I generate a mesh:

regionmesh = 
ToElementMesh[region, "MaxBoundaryCellMeasure" -> 0.05/10]

And the equation to be solved :

equation = 
D[t[x, y, z], z] == 
D[t[x, y, z], {x, 2}] + D[t[x, y, z], {y, 2}] + 
NeumannValue[
10*(150 - t[x, y, z]), Element[{x, y, z}, tubeRegion]]

With the initial value:

dirichlet=DirichletCondition[
t[x, y, z] == 
170, Element[{x, y},
RegionDifference[Disk[{0, 0}, shellID], tubeSection] && z == 0]]

The equation is then solved:

solmesh = 
NDSolveValue[{equation, dirichlet}, 
t, Element[{x, y, z}, regionmesh]];

The procedure works but the mesh does not look like it accurately reproduces the tube region. Besides, I'm specifying the initial value of the temeprature in a way that seems inconsistent with the mesh. I think it would be better if I could directly specify it on the mesh itself rather than having to specify a geometric section. Any suggestions to make this code more efficient ?

Answer 1

One thing you could do is to generate a 2D mesh and the extrude the 3D version from it.

mr = BoundaryDiscretizeRegion[RegionDifference[Disk[], tubeSection]
  (*,AccuracyGoal -> 3,PrecisionGoal -> 4*)]

Make a region product with a line:

rp = RegionProduct[mr, Line[{{0}, {reactorLength}}]]

And then you can generate a FEM mesh from it if you want:

mesh = ToElementMesh[rp];
mesh["Wireframe"]

Note, that even though this is a second order mesh, the boundary faces will not be curved.