(*Element info shortcuts*)
ebi = ElementIncidents[#["BoundaryElements"]][[1]] &;
ebm = ElementMarkers[#["BoundaryElements"]][[1]] &;
ebn = #["BoundaryNormals"][[1]] &;
ei = ElementIncidents[#["MeshElements"]][[1]] &;
em = ElementMarkers[#["MeshElements"]][[1]] &;
epi = ElementIncidents[#["PointElements"]][[1]] &;
epm = Flatten@ElementMarkers[#["PointElements"]] &;
(*extract boundary mesh from element mesh*)
bmesh = ToBoundaryMesh[mesh];
bcrd = bmesh["Coordinates"];
bi = ebi[bmesh];(*boundary element incidents*)
bm = ebm[bmesh];(*boundary element markers*)
bn = ebn[bmesh];(*boundary normals*)
(*find markers corresponding to the spherical cap*)
mrk3pos = Flatten@Position[bm, 3, 1];
(*generate necessary info to estimate surface integral*)
bn3 = bn[[mrk3pos]];
polys = Map[Polygon,
GetElementCoordinates[bcrd, #] & /@ bi[[mrk3pos]]];
area3 = Area /@ polys;
center3 = Map[Mean, GetElementCoordinates[bcrd, #] & /@ bi[[mrk3pos]]];
(*find positions of left and right side of spherical cap*)
posXids = Position[center3Position[center3[[All, 1]], _?(#[[1]]# >= 0 &), 1] // Flatten;
negXids = Complement[Range[Length[mrk3pos]], posXids];
Show[{Graphics3D[{Red, polys[[posXids]]}],
Graphics3D[{Blue, polys[[negXids]]}]}]
So, as the Péclet number goes down the Advective component becomes less significant indicating that a "wind tunnel" is perhaps not the best model to study this problem. Although you may want to consider modifying the geometry, you can substantially mitigate the effects of low Péclet numbers by changing the default wall conditions to DirichletCondition's. For exampleHere is a logarithmic sweep from Pe numbers from 0.01 to 100 (note that this process is slow):
Pe = 0.001;
solpfun = NDSolveValue[ParametricNDSolveValue[{D[T[x, y, z], x] ==
1/PeP Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 3],
DirichletCondition[T[x, y, z] == 0., ElementMarker == 4 ||
Or @@ (ElementMarker == 6 ||
# & /@ ElementMarker{4, ==6, 7]7})]}},
T, {x, y, z} ∈ mesh];mesh, {P}];
surf = {{x^2 + y^2 + z^2 ==
1.001^2}, {"XStackedPlanes", {7.5}}, {"YStackedPlanes", {0}}, \
{"ZStackedPlanes", {-0.8}}, {"BackPlanes"}};
SliceContourPlot3D[sol[xframes = SliceContourPlot3D[pfun[#][x, y, z],
surf, {x, y, z} ∈ mesh,
Contours -> 11,
PlotPoints -> 100, BoxRatios -> Automatic,
ColorFunction -> "ThermometerColors", PlotRange -> {-0.001, 1},
PlotLegends -> Automatic]Automatic, PlotLabel -> N@#] & /@ (10^# & /@
Subdivide[-2, 2, 20]);
ListAnimate[frames]
So, as the Péclet number goes down the Advective component becomes less significant indicating that a "wind tunnel" is perhaps not the best model to study this problem. Although you may want to consider modifying the geometry, you can substantially mitigate the effects of low Péclet numbers by changing the default wall conditions to DirichletCondition's. For example:
Pe = 0.001;
sol = NDSolveValue[{D[T[x, y, z], x] ==
1/Pe Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 3],
DirichletCondition[T[x, y, z] == 0., ElementMarker == 4 ||
ElementMarker == 6 ||
ElementMarker == 7]}}, T, {x, y, z} ∈ mesh];
surf = {{x^2 + y^2 + z^2 ==
1.001^2}, {"XStackedPlanes", {7.5}}, {"YStackedPlanes", {0}}, \
{"ZStackedPlanes", {-0.8}}, {"BackPlanes"}};
SliceContourPlot3D[sol[x, y, z], surf, {x, y, z} ∈ mesh,
Contours -> 11, PlotPoints -> 100, BoxRatios -> Automatic,
ColorFunction -> "ThermometerColors", PlotRange -> {-0.001, 1},
PlotLegends -> Automatic]
So, as the Péclet number goes down the Advective component becomes less significant indicating that a "wind tunnel" is perhaps not the best model to study this problem. Although you may want to consider modifying the geometry, you can substantially mitigate the effects of low Péclet numbers by changing the default wall conditions to DirichletCondition's. Here is a logarithmic sweep from Pe numbers from 0.01 to 100 (note that this process is slow):
pfun = ParametricNDSolveValue[{D[T[x, y, z], x] ==
1/P Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 3],
DirichletCondition[T[x, y, z] == 0.,
Or @@ (ElementMarker == # & /@ {4, 6, 7})]}},
T, {x, y, z} ∈ mesh, {P}];
surf = {{x^2 + y^2 + z^2 ==
1.001^2}, {"XStackedPlanes", {7.5}}, {"YStackedPlanes", {0}}, \
{"ZStackedPlanes", {-0.8}}, {"BackPlanes"}};
frames = SliceContourPlot3D[pfun[#][x, y, z],
surf, {x, y, z} ∈ mesh, Contours -> 11,
PlotPoints -> 100, BoxRatios -> Automatic,
ColorFunction -> "ThermometerColors", PlotRange -> {-0.001, 1},
PlotLegends -> Automatic, PlotLabel -> N@#] & /@ (10^# & /@
Subdivide[-2, 2, 20]);
ListAnimate[frames]