This is a question similar to the previous one with small features. I'm solving this problem by modifying the method of @Alex. The code can run without error, but it was always stuck at ~t=2.328e-14 for a long time (no progress over 2 hrs). I guess it is due to the low mesh quality. To improve the mesh, I have also tried to increase the resolution using MaxCellMeasure-> 5*10^-5. I also added "MaxBoundaryCellMeasure" -> 0.01 and AccuracyGoal -> 5. But I still have no luck ...
The following is my full code. Though it looks long but it is very straightforward, and I also added some comments to make it clearer. Please help me out. And any suggestions are welcome. Thank you!
(*using scale L = max length to rescale all geometrical dimensions*)
Needs["NDSolve`FEM`"];
L = 2800;(*mm, length scale*)
sizes = {length -> 1900, height -> 2800, ductW -> 30, offsetx1 -> 5, offsety1 -> 1000, offsety2 -> 15, offsety3 -> 30, offsety4 -> 45, slab -> 10, r -> 5/2};(*dimensional sizes in mm*)
(*non-dimensionalization*)
sizes[[All, 2]] = sizes[[All, 2]]/L;
\[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes;
\[CapitalOmega]d1 = Rectangle[{0, offsety1}, {ductW, slab + offsety1}] /. sizes;
\[CapitalOmega]d2 = Rectangle[{0, offsety1 + slab + offsety4 + offsety2}, {ductW, offsety1 + 2*slab + offsety4 + offsety2}] /. sizes;
\[CapitalOmega]d = RegionUnion[\[CapitalOmega]d1, \[CapitalOmega]d2];
(*scaled heat source surfaces*)
HSsurf = First[{(100/L <= x <= 1850/L && y == 600/L) || (600/L <= y <= 750/L && x == 100/L) || (600/L <= y <= 750/L && x == 1850/L) || (100/L <= x <= 1850/L && y == 750/L)}];
(*scaled duct surfaces*)
ds = First@{(0 <= x <= ductW && y == offsety1) || (0 <= x <= ductW && y == offsety1 + slab) || (offsety1 <= y <= offsety1 + slab && x == ductW) || (offsety1 + slab + offsety4 + offsety2 <= y <= offsety1 + 2*slab + offsety4 + offsety2 && x == ductW) || (0 <= x <= ductW && y == offsety1 + slab + offsety4 + offsety2) || (0 <= x <= ductW && y == offsety1 + 2*slab + offsety4 + offsety2)} /. sizes;
tubes = {Disk[{offsetx1, offsety1 + slab + offsety2}, r], Disk[{offsetx1, offsety1 + slab + offsety3}, r], Disk[{offsetx1, offsety1 + slab + offsety4}, r], Disk[{3*offsetx1, offsety1 + slab + offsety2}, r], Disk[{3*offsetx1, offsety1 + slab + offsety3}, r], Disk[{3*offsetx1, offsety1 + slab + offsety4}, r], Disk[{5*offsetx1, offsety1 + slab + offsety2}, r], Disk[{5*offsetx1, offsety1 + slab + offsety3}, r], Disk[{5*offsetx1, offsety1 + slab + offsety4}, r]} /. sizes;
(*tube surfaces*)
t1 = Table[(x - n*offsetx1)^2 + (y - (offsety1 + slab + offsety2))^2 <= r^2, {n, {1, 3, 5}}] /. sizes;
t2 = Table[(x - n*offsetx1)^2 + (y - (offsety1 + slab + offsety3))^2 <= r^2, {n, {1, 3, 5}}] /. sizes;
t3 = Table[(x - n*offsetx1)^2 + (y - (offsety1 + slab + offsety4))^2 <= r^2, {n, {1, 3, 5}}] /. sizes;
(*put the surfaces into a form appropriate for setting the bcs*)
ts = t1[[1]] || t1[[2]] || t1[[3]] || t2[[1]] || t2[[2]] || t2[[3]] || t3[[1]] || t3[[2]] || t3[[3]];
\[CapitalOmega]3 = RegionDifference[\[CapitalOmega]1, RegionUnion[Flatten@tubes]];
HS = ImplicitRegion[100/L <= x <= 1850/L && 600/L <= y <= 750/L, {x, y}];(*heat source region*)
\[CapitalOmega]4 = RegionDifference[\[CapitalOmega]3, HS];
\[CapitalOmega] = RegionDifference[\[CapitalOmega]4, \[CapitalOmega]d];
(*observe the domain and small features*)
RegionPlot[\[CapitalOmega], PlotPoints -> 80, MaxRecursion -> 2, BoundaryStyle -> Black, PlotStyle -> LightBlue, AspectRatio -> Automatic]
smallReg = {{0, 1/56}, {7/20, 11/28}};
RegionPlot[\[CapitalOmega], PlotPoints -> 80, MaxRecursion -> 2, BoundaryStyle -> Black, PlotStyle -> LightBlue, AspectRatio -> Automatic, PlotRange -> smallReg]
mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 5*10^-5, "MaxBoundaryCellMeasure" -> 0.01, AccuracyGoal -> 5];
mesh["Wireframe"]
Show[mesh["Wireframe"["MeshElementStyle" -> {Directive[FaceForm[Blue]]}]], PlotRange -> smallReg]
parameters = {Pr -> 0.72, Ra -> 10^6};
op = {Derivative[1, 0, 0][u][t, x, y] + Inactive[Div][-Pr*Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + {u[t, x, y], v[t, x, y]} . Inactive[Grad][u[t, x, y], {x, y}] + Derivative[0, 1, 0][p][t, x, y],
Derivative[1, 0, 0][v][t, x, y] + Inactive[Div][-Pr*Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + {u[t, x, y], v[t, x, y]} . Inactive[Grad][v[t, x, y], {x, y}] + Derivative[0, 0, 1][p][t, x, y] - Ra*Pr*T[t, x, y],
Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y],
Derivative[1, 0, 0][T][t, x, y] + Inactive[Div][-Inactive[Grad][T[t, x, y], {x, y}], {x, y}] + {u[t, x, y], v[t, x, y]} . Inactive[Grad][T[t, x, y], {x, y}]} /. parameters;
(*BCs and ICs*)
walls = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, ts || x == length || y == 0 || y == height || ds || HSsurf]} /. sizes;
inletflow = DirichletCondition[{u[t, x, y] == 1, v[t, x, y] == 0}, x == 0 && offsety1 + slab <= y <= offsety1 + slab + offsety4 + offsety2] /. sizes;
Subscript[\[CapitalGamma], outlet] = DirichletCondition[p[t, x, y] == 0, (x == 0 && 0 <= y <= offsety1) || (x == 0 && offsety1 + 2*slab + offsety4 + offsety2 <= y <= height)] /.sizes;
reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];
bcsflow = {walls, inletflow, Subscript[\[CapitalGamma], outlet], reference};
bcstemperatures = {DirichletCondition[T[t, x, y] == 12/17, x == 0 && offsety1 + slab <= y <= offsety1 + slab + offsety4 + offsety2] /. sizes(*inlet temp*),
DirichletCondition[T[t, x, y] == 0, ts],
DirichletCondition[T[t, x, y] == 1, HSsurf]};
ic = {u[0, x, y] == 1, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 12/17};
tmax = 1;
Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} =
NDSolveValue[Flatten[{op == {0, 0, 0, 0}, bcsflow, bcstemperatures, ic}], {u, v, p, T}, Element[{x, y}, mesh], {t, 0, tmax}, Method -> {"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2}, "PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}, EvaluationMonitor :> (currentTime = Row[{"t=", CForm[t]}])];], currentTime]
Thank you for your time!
ToElementMeshhas the optionMeshRefinementFunctionfor this purpose. $\endgroup$