# Meshing an irregular domain using quads to solve conjugate heat transfer problem

I am trying to mesh the following domain to solve a heat transfer + fluid flow problem:

The continuity+momentum equations are to be solved in $$ABGH$$, while the energy equation is to be solved across the entire domain, i.e., $$ABCDEFGH$$ (basically line $$CF$$ is the solid-fluid interface). To generate mesh for $$ABGH$$, I used the following:

Needs["MeshTools"]
L = 0.050, d = 0.003, e = 0.005, delta = 0.010, l = L/d + 2 delta/d
mesh = ToElementMesh[FullRegion[2], {{0, l}, {0, 1}}, MaxCellMeasure -> 8 10^-3];

(*ElementMesh[{{0., 23.3333}, {0., 1.}}, {QuadElement["<" 3132 ">"]}]*)


which produces the following:

Now, to mesh the entire domain, I use the tool StructuredMesh in the following code (about which I learned from the various answers on this post):

raster = {{{0, 0}, {delta/d, 0}, {delta/d, -e/d}, {delta/d + L/d, -e/d}, {delta/d + L/d, 0}, {l, 0}}, {{0, 1}, {delta/d, 1}, {delta + L/(3 d), 1}, {delta + (2 L/(3 d)), 1}, {delta/d + L/d, 1}, {l, 1}}}
mesh2 = StructuredMesh[raster, {200, 20}]
mesh2["Wireframe"]
(*ElementMesh[{{0., 23.3333}, {-1.66667, 1.}}, {QuadElement["<" 4000 ">"]}]*)


This produces the following:

To produce this, I have basically defined the domain using two lines, each having six points through the raster command. However, as visible the mesh quality is pretty bad.

How should the $$ABCDEFGH$$ domain be meshed to achieve a smooth quad mesh as $$ABGH$$? Also, it will be beneficial if I can have the same number of elements in the ABGH portion, while meshing ABCDEFGH as it had when it was individually meshed.

Also, this leads me to think, can both these sub-domains be meshed separetly and then joined?

Update 1 Upon recommendations in the comment, I tried the following:

MergeMesh for different size domains work on BoundaryMesh type object. Atleast, that is what I found in the documentation.

Needs["NDSolveFEM"]
Needs["MeshTools"]
{L = 0.050, d = 0.003, e = 0.005, delta = 0.010, l = L/d + 2 delta/d};
reg1 = ImplicitRegion[0 <= x <= l && 0 <= y <= 1, {x, y}];
reg2 = ImplicitRegion[
delta/d <= x <= l - delta/d && -e/d <= y <= 0, {x, y}];
reg3 = RegionUnion[reg1, reg2];

bm1 = ToBoundaryMesh[reg1]
bm2 = ToBoundaryMesh[reg2]
bm3 = MergeMesh[{bm1, bm2}]
meshfluid = ToElementMesh[bm1, MaxCellMeasure -> 0.01]["Wireframe"]
meshenergy = ToElementMesh[bm3, MaxCellMeasure -> 0.01]["Wireframe"]
(*meshsolid = ToElementMesh[bm2, MaxCellMeasure -> 0.01]["Wireframe"]*)



However, this is a triangular mesh. TriangleToQuadMesh does not work on meshenergy object. Also if you see the image above closely, the refinement at the interface (i.e., internal boundary) is localised to a few regions. Also MeshCellCount does not work when I try it on meshenergy, meshfluid or meshsolid.

Update 2

A comparison among tri and quad mesh results

I used the answer given here by Oleksii, to solve the problem of conjugate heat transfer (reciprocating flow over a heated block). Obviously, I have utlised the domain described in this question. Both the cases have been run for a flow time of 5 seconds, with similar MaxStepSize setting

Triangular mesh results

Created using the Update 1 in this question.

Please note that ToElementMesh creates a 2nd order mesh. To interpolate the p which is supplied with p->1 in NDSolve, I had to create a 1st order mesh using MeshOrderAlteration

1. Mesh statistics

1. Solid temperature variation at point ((delta+L)/d/2,-e/d/2)

1. Fluid temperature variation at both the inlets of the fluid domain

Mesh generated using Oleksii's answer in this question

Similar number of elements in y-direction in the fluid domain (compared to trimesh) could be achieved with the setting NyF. Hence, total number of elements is lesser than trimesh.

1. Mesh statistics

1. Solid temperature variation at point ((delta+L)/d/2,-e/d/2)

1. Fluid temperature variation at both the inlets of the fluid domain

In conclusion, the results seem fairly simialr. I will run these tests for longer flow times and report if there are any discrepancies that arise. However, it seems quadmesh allows more control over element size in such regular domains (thus leading to lesser total element count).

Addendum to Update 2 The following contour plots, show the difference in solution between the quad and triangle mesh (at the same flow time). Please note that the triangle mesh is generated using the answer of @user21.

Solid temperature

Fluid temperature Have replaced AspectRatio->Automatic to AspectRatio->1/2, as the fluid region is thin.

Note: Total elements for triangle mesh: 3246 Total elements for quad mesh: 1248

• You can try to build uniform meshes separately for fluid and solid domains and then to join them by MergeMesh from MeshTools Feb 4 at 12:11
• @OleksiiSemenov I gave it a try. See Update. Feb 4 at 13:08
• Great. A few comments. You should not use the mesh from update 1 - this already looks bad and I suspect it has holes or some such. For the triangular mesh use the mesh that I used for the quad mesh generation. Feb 10 at 8:31
• If you want to put in the work, what would be much more informative is to plot the difference between the solutions. But again, for this type of problem I would not expect much of a difference in the solution. Triangle elements are typically a bit more efficient to evaluate. Feb 10 at 8:33
• @user21 Thankyou for the comments. See "Addendum to Update 2". The fluid domain shows considerable difference it seems. Feb 11 at 3:36

Function MergeMesh from MeshTools package allows to join separate meshes. Let's take advantage of this useful function for this problem.

Needs["MeshTools"]
L = 0.050; d = 0.003; e = 0.005; delta = 0.010;
l = L/d + 2 delta/d;

NxS = 40;
NyS = 10;
MeshSolid =
RectangleMesh[{delta/d, -e/d}, {(L + delta)/d, 0}, {NxS, NyS}];

NxF = 10;
NyF = 10;
MeshFluid1 = RectangleMesh[{0, 0}, {delta/d, 1}, {NxF, NyF}];
MeshFluid2 = RectangleMesh[{(delta + L)/d, 0}, {l, 1}, {NxF, NyF}];
MeshFluid3 = RectangleMesh[{(delta)/d, 0}, {(L + delta)/d, 1}, {NxS, NyF}];
MeshFluid = MergeMesh[MergeMesh[MeshFluid1, MeshFluid3], MeshFluid2]

MeshTotal = MergeMesh[MeshFluid, MeshSolid];
MeshTotal["Wireframe"]


• Much much thanks. This is exactly what I wanted. I wasnt aware of the RectangleMesh function too. Keeping similar NxS probably helped. I would also like to ask, if such rectangle meshes are better for fluid flow problems ? Feb 4 at 19:20
• @Avrana, it would be interesting to hear what you find when you use this mesh in place of a triangulated region. I doubt it is going to make much of a difference, but I could be wrong. Feb 6 at 14:49
• @Avrana I will respond on your question about tri mesh vs quads mesh a little bit later. It is blackout in my city now and I am on smartphone Feb 6 at 15:23
• @Avrana about tri mesh vs quad mesh in cfd. It is popular debatable issue. And as far as I know the choice of proper mesh depends on particular problem in question. May be this work researchgate.net/publication/… will be interesting for you. As concerns the problem you are working on, I think you should carry out the computational experiments with different meshes and then make a choice. Interesting to look at these results. Feb 7 at 21:23
• @Avrana The points $\{\frac{\delta}{d},0\}$ and $\{\frac{L+\delta}{d},0\}$ are certainly the trouble spots where the large temperature gradients appears. And it would be better to make gradual mesh refinement there. Feb 7 at 21:28

A different approach might be to use ToQuadMesh:

Needs["FEMAddOns"]
L = 0.050; d = 0.003; e = 0.005; delta = 0.010; l = L/d + 2 delta/d;
bmesh = ToBoundaryMesh[
"Coordinates" -> {{0, 0}, {delta/d,
0}, {delta/d, -e/d}, {delta/d + L/d, -e/d}, {delta/d + L/d,
0}, {2*delta/d + L/d, 0}, {L/d + 2*delta/d, 1}, {0, 1}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
5}, {5, 6}, {6, 7}, {7, 8}, {8, 1}, {2, 5}}]}];
bmesh["Wireframe"]


mesh = ToElementMesh[bmesh];

Not, quite what you were looking for but perhaps an alternative option. Note, that the internal boundary is not respected by ToQuadMesh. If you specify region markers then it works:
mesh = ToElementMesh[bmesh,

• Thankyou for this. I will test this mesh too. I hope the RegionMarker can be used to assign separate properties to each zone? Feb 10 at 4:48
• The RegionMarker can indeed do this. The documentation gives some pretty good examples of use cases. Feb 10 at 5:33