Speed up finite element method solver with small features

Question

I post a question on making mesh for small features. @user21 gave a useful suggestion which resolves the small features. To improve the mesh, I have also tried to increase the resolution near the boundaries using MaxBoundaryCellMeasure. Also, I need a higher resolution in the center than that was given with default options (but can be coarser than those near the boundaries), thus adding MaxCellMeasure -> 100.

 mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 100, "MaxBoundaryCellMeasure" -> 5(*1*), AccuracyGoal -> 3]; (*{mesh["Wireframe"[PlotRange -> {{0, 20}, {240, 280}}]], mesh["Wireframe"[PlotRange -> All]]}*)

Now, I follow the answer provided by @Alex Trounev to solve the same set of equations. The cooling is from an array of small disks with a constant temperature and heating by a heat source with constant heat flux. For completeness, I post the codes together.

Needs["NDSolve`FEM`"]; ResourceFunction["FEMAddOnsInstall"][]; Needs["FEMAddOns`"] (*geometry*) sizes = {length -> 300, height -> 280, ductW -> 10, offsetx1 -> 5, offsetx2 -> 130, offsety1 -> 30, offsety2 -> 1, offsety3 -> 4, slab -> 1, r -> 8/10}; \[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes; \[CapitalOmega]d1 = Rectangle[{ductW, height - offsety1 - slab}, {ductW + 3*slab + offsetx1, height - offsety1}] /. sizes; \[CapitalOmega]d2 = Rectangle[{ductW, height/2}, {ductW + slab, height - offsety1}] /. sizes; \[CapitalOmega]d = RegionUnion[\[CapitalOmega]d1, \[CapitalOmega]d2]; (*duct surfaces*) ds = First@{(ductW <= x <= ductW + slab && y == height/2) || (ductW <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1) || (height/2 <= y <= height - offsety1 && x == ductW) || (height/2 <= y <= height - offsety1 - slab && x == ductW + slab) || (height - offsety1 - slab <= y <= height - offsety1 && x == ductW + 3*slab + offsetx1) || (ductW + slab <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1 - slab)} /. sizes tubes = Table[{Disk[{ductW + slab, height - offsety2 - n*offsety3}, r], Disk[{ductW + slab + offsetx1, height - offsety2 - n*offsety3}, r]}, {n, 0, 7}] /. sizes; (*tube surfaces*) t1 = Table[(x - (ductW + slab))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; t2 = Table[(x - (ductW + slab + offsetx1))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; (*put into a form appropriate for bcs*) ts = t1[[1]] || t1[[2]] || t1[[3]] || t1[[4]] || t1[[5]] || t1[[6]] || t1[[7]] || t1[[8]] || t2[[1]] || t2[[2]] || t2[[3]] || t2[[4]] || t2[[5]] || t2[[6]] || t2[[7]] || t2[[8]]; \[CapitalOmega]3 = RegionDifference[\[CapitalOmega]1, RegionUnion[Flatten@tubes]]; \[CapitalOmega] = RegionDifference[\[CapitalOmega]3, \[CapitalOmega]d]; (*mesh*) mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 100, "MaxBoundaryCellMeasure" -> 5(*1*), AccuracyGoal -> 3]; (*heat source surfaces*) HSsurf = First@{(130 <= x <= 170 && y == 170) || (0 <= y <= 170 && x == 130) || (0 <= y <= 170 && x == 170)}; HS = With[{q = 10}, If[130 <= x <= 170 && 0 <= y <= 170, q, 0]]; (*main*) parameters = {Pr -> 0.7, Ra -> 10^6}; op = { \!\(\*SuperscriptBox[\(u\), TagBox[ RowBox[{"(", RowBox[{"1", ",", "0", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[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}] + \!\(\*SuperscriptBox[\(p\), TagBox[ RowBox[{"(", RowBox[{"0", ",", "1", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[t, x, y], \!\(\*SuperscriptBox[\(v\), TagBox[ RowBox[{"(", RowBox[{"1", ",", "0", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[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}] + \!\(\*SuperscriptBox[\(p\), TagBox[ RowBox[{"(", RowBox[{"0", ",", "0", ",", "1"}], ")"}], Derivative], MultilineFunction->None]\)[t, x, y] - Ra*Pr*T[t, x, y], \!\(\*SuperscriptBox[\(u\), TagBox[ RowBox[{"(", RowBox[{"0", ",", "1", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[t, x, y] + \!\(\*SuperscriptBox[\(v\), TagBox[ RowBox[{"(", RowBox[{"0", ",", "0", ",", "1"}], ")"}], Derivative], MultilineFunction->None]\)[t, x, y], \!\(\*SuperscriptBox[\(T\), TagBox[ RowBox[{"(", RowBox[{"1", ",", "0", ",", "0"}], ")"}], Derivative], MultilineFunction->None]\)[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 == 0 || x == length || y == 0 || y == height || ds || HSsurf]} /. sizes; reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]; bcsflow = {walls, reference}; bcstemperatures = {DirichletCondition[T[t, x, y] == 0, ts]}; ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0}; tmax = 1; Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = NDSolveValue[Flatten@{op == {0, 0, 0, HS}, bcsflow, bcstemperatures, ic}, {u, v, p, T}, {x, y} \[Element] 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]

All boundary conditions looked fine. The code can run without error but is extremely slow after t=3.81..e^-10.

Questions:

Please help to speed up the code. I'm unsure if I used the correct method to solve this problem. For testing, we can reduce the number of small disks.
Because most of the domain features straight boundaries, should we discretize the domain into a few sub-domains with some regions having quadrilateral meshes? Would such a mesh be more accurate for FEM simulation?

I appreciate any help you can provide.

This is not a valid boundary condition: DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, {x, y}]. Also this is problematic: DirichletCondition[T[t, x, y] == 0, {x, y} [Element] RegionUnion[Flatten@tubes]]. Have a look at the reference page of DirichletCondition or the PDEModels and have a look there at what are good ways to model Dirichlet BCs. — user21
– user21, Commented Jul 17, 2024 at 8:10
I have voted to close this question, as one of your BCs contains an invalid predicate and another uses a numerically unstable method for the predicate. — user21
– user21, Commented Jul 18, 2024 at 8:28
@user21, thank you for the suggestion. I read the relevant doc and modified the the Dirichlet bc. Now the code can run without any error but extremely slow. — user55777
– user55777, Commented Jul 20, 2024 at 12:41
@user55777 Region $\Omega $ looks very simply, while it defined by very complicated way. Also, please note, that in convection problem we use length to define Ra. Therefore, region in nondimensional form has length=1. Since you take length -> 300, effectively in your problem Ra= 10^6 length^3=2.7 10^13. No one can solve this problem:). — Alex Trounev
– Alex Trounev, Commented Jul 21, 2024 at 10:18
@AlexTrounev, (1) The region is simple but with some small features. To add these features, the definition looks complicated, please refer to my previous question. We may reduce the number of disks to, say, 6. (2) Yes, Ra is proportional to l^3. Here, length -> 300 is the width of the domain, but not necessary to be the characteristic length. Could we take a smaller l to be a characteristic length, say the width of the disk array, offsetx1->5? — user55777
– user55777, Commented Jul 21, 2024 at 10:31

Alex Trounev · Accepted Answer · 2024-07-24 04:56:23Z

To solve this problem, we put characteristic scale L=300 and then rescale all geometrical parameters on L, as result we have

Needs["NDSolve`FEM`"]; (*geometry*)L = 300; sizes = {length -> 300, height -> 280, ductW -> 10, offsetx1 -> 5, offsetx2 -> 130, offsety1 -> 30, offsety2 -> 1, offsety3 -> 4, slab -> 1, r -> 8/10}; sizes[[All, 2]] = sizes[[All, 2]]/L; sizes; \[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes; \[CapitalOmega]d1 = Rectangle[{ductW, height - offsety1 - slab}, {ductW + 3*slab + offsetx1, height - offsety1}] /. sizes; \[CapitalOmega]d2 = Rectangle[{ductW, height/2}, {ductW + slab, height - offsety1}] /. sizes; \[CapitalOmega]d = RegionUnion[\[CapitalOmega]d1, \[CapitalOmega]d2]; (*duct surfaces*) ds = First@{(ductW <= x <= ductW + slab && y == height/2) || (ductW <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1) || (height/2 <= y <= height - offsety1 && x == ductW) || (height/2 <= y <= height - offsety1 - slab && x == ductW + slab) || (height - offsety1 - slab <= y <= height - offsety1 && x == ductW + 3*slab + offsetx1) || (ductW + slab <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1 - slab)} /. sizes; tubes = Table[{Disk[{ductW + slab, height - offsety2 - n*offsety3}, r], Disk[{ductW + slab + offsetx1, height - offsety2 - n*offsety3}, r]}, {n, 0, 7}] /. sizes; (*tube surfaces*) t1 = Table[(x - (ductW + slab))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; t2 = Table[(x - (ductW + slab + offsetx1))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; (*put into a form appropriate for bcs*) ts = t1[[1]] || t1[[2]] || t1[[3]] || t1[[4]] || t1[[5]] || t1[[6]] || t1[[7]] || t1[[8]] || t2[[1]] || t2[[2]] || t2[[3]] || t2[[4]] || t2[[5]] || t2[[6]] || t2[[7]] || t2[[8]]; \[CapitalOmega]3 = RegionDifference[\[CapitalOmega]1, RegionUnion[Flatten@tubes]]; \[CapitalOmega] = RegionDifference[\[CapitalOmega]3, \[CapitalOmega]d]; (*mesh*) mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> .01, "MaxBoundaryCellMeasure" -> .01(*1*), AccuracyGoal -> 5]; mesh["Wireframe"]

HSsurf = First[{(130/L <= x <= 170/L && y == 170/L) || (0 <= y <= 170/L && x == 130/L) || (0 <= y <= 170/L && x == 170/L)}]; HS = With[{q = 10}, If[130/L <= x <= 170/L && 0 <= y <= 170/L, q, 0]]; parameters = {Pr -> 0.7, 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; walls = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, ts || x == 0 || x == length || y == 0 || y == height || ds || HSsurf]} /. sizes; reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]; bcsflow = {walls, reference}; bcstemperatures = {DirichletCondition[T[t, x, y] == 0, ts]}; ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0}; tmax = 1; Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = NDSolveValue[Flatten[{op == {0, 0, 0, HS}, 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]

It takes about 450s on my Silver Pentium. Visualization

ContourPlot[temperature[1, x, y], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, ContourStyle -> None, Contours -> 20]

Show[DensityPlot[ Norm[{xVel[1, x, y], yVel[1, x, y]}], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, PlotPoints -> 100], StreamPlot[{xVel[1, x, y], yVel[1, x, y]}, {x, y} \[Element] mesh, ColorFunction -> None, StreamStyle -> LightGray]]

Update 1. To simulate the source as a solid-block heater we can use as a first step next model

Needs["NDSolve`FEM`"]; (*geometry*)L = 300; sizes = {length -> 300, height -> 280, ductW -> 10, offsetx1 -> 5, offsetx2 -> 130, offsety1 -> 30, offsety2 -> 1, offsety3 -> 4, slab -> 1, r -> 8/10}; sizes[[All, 2]] = sizes[[All, 2]]/L; sizes; \[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes; \[CapitalOmega]d1 = Rectangle[{ductW, height - offsety1 - slab}, {ductW + 3*slab + offsetx1, height - offsety1}] /. sizes; \[CapitalOmega]d2 = Rectangle[{ductW, height/2}, {ductW + slab, height - offsety1}] /. sizes; \[CapitalOmega]d = RegionUnion[\[CapitalOmega]d1, \[CapitalOmega]d2]; (*duct surfaces*) ds = First@{(ductW <= x <= ductW + slab && y == height/2) || (ductW <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1) || (height/2 <= y <= height - offsety1 && x == ductW) || (height/2 <= y <= height - offsety1 - slab && x == ductW + slab) || (height - offsety1 - slab <= y <= height - offsety1 && x == ductW + 3*slab + offsetx1) || (ductW + slab <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1 - slab)} /. sizes; tubes = Table[{Disk[{ductW + slab, height - offsety2 - n*offsety3}, r], Disk[{ductW + slab + offsetx1, height - offsety2 - n*offsety3}, r]}, {n, 0, 7}] /. sizes; (*tube surfaces*) t1 = Table[(x - (ductW + slab))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; t2 = Table[(x - (ductW + slab + offsetx1))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; (*put into a form appropriate for bcs*) ts = t1[[1]] || t1[[2]] || t1[[3]] || t1[[4]] || t1[[5]] || t1[[6]] || t1[[7]] || t1[[8]] || t2[[1]] || t2[[2]] || t2[[3]] || t2[[4]] || t2[[5]] || t2[[6]] || t2[[7]] || t2[[8]]; HSsurf = First[{(130/L <= x <= 170/L && y == 170/L) || (0 <= y <= 170/L && x == 130/L) || (0 <= y <= 170/L && x == 170/L)}]; HS = ImplicitRegion[130/L <= x <= 170/L && 0 <= y <= 170/L, {x, y}]; \[CapitalOmega]3 = RegionDifference[\[CapitalOmega]1, RegionUnion[Flatten@tubes]]; \[CapitalOmega]4 = RegionDifference[\[CapitalOmega]3, HS]; \[CapitalOmega] = RegionDifference[\[CapitalOmega]4, \[CapitalOmega]d]; (*mesh*) mesh = ToElementMesh[\[CapitalOmega], MaxCellMeasure -> .01, "MaxBoundaryCellMeasure" -> .01(*1*), AccuracyGoal -> 5]; mesh["Wireframe"]

parameters = {Pr -> 0.7, 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; walls = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, ts || x == 0 || x == length || y == 0 || y == height || ds || HSsurf]} /. sizes; reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]; bcsflow = {walls, reference}; bcstemperatures = {DirichletCondition[T[t, x, y] == 0, ts],DirichletCondition[T[t, x, y] == 1, HSsurf]}; ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0}; 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]

It takes about 222s. Visualization

ContourPlot[temperature[1, x, y], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, ContourStyle -> None, Contours -> 20]

Show[DensityPlot[ Norm[{xVel[1, x, y], yVel[1, x, y]}], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, PlotPoints -> 100], StreamPlot[{xVel[1, x, y], yVel[1, x, y]}, {x, y} \[Element] mesh, ColorFunction -> None, StreamStyle -> LightGray]]

Update 2. We can resolve small disks even better using MeshRefinementFunction to generate mesh as follows

(*geometry*)L = 300; sizes = {length -> 300, height -> 280, ductW -> 10, offsetx1 -> 5, offsetx2 -> 130, offsety1 -> 30, offsety2 -> 1, offsety3 -> 4, slab -> 1, r -> 8/10}; sizes[[All, 2]] = sizes[[All, 2]]/L; \[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes; \[CapitalOmega]d1 = Rectangle[{ductW, height - offsety1 - slab}, {ductW + 3*slab + offsetx1, height - offsety1}] /. sizes; \[CapitalOmega]d2 = Rectangle[{ductW, height/2}, {ductW + slab, height - offsety1}] /. sizes; \[CapitalOmega]d = RegionUnion[\[CapitalOmega]d1, \[CapitalOmega]d2]; (*duct surfaces*) ds = First@{(ductW <= x <= ductW + slab && y == height/2) || (ductW <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1) || (height/2 <= y <= height - offsety1 && x == ductW) || (height/2 <= y <= height - offsety1 - slab && x == ductW + slab) || (height - offsety1 - slab <= y <= height - offsety1 && x == ductW + 3*slab + offsetx1) || (ductW + slab <= x <= ductW + 3*slab + offsetx1 && y == height - offsety1 - slab)} /. sizes; tubes = Table[{Disk[{ductW + slab, height - offsety2 - n*offsety3}, r], Disk[{ductW + slab + offsetx1, height - offsety2 - n*offsety3}, r]}, {n, 0, 7}] /. sizes; (*tube surfaces*) t1 = Table[(x - (ductW + slab))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; t2 = Table[(x - (ductW + slab + offsetx1))^2 + (y - (height - offsety2 - n*offsety3))^2 <= r^2, {n, 0, 7}] /. sizes; (*put into a form appropriate for bcs*) ts = t1[[1]] || t1[[2]] || t1[[3]] || t1[[4]] || t1[[5]] || t1[[6]] || t1[[7]] || t1[[8]] || t2[[1]] || t2[[2]] || t2[[3]] || t2[[4]] || t2[[5]] || t2[[6]] || t2[[7]] || t2[[8]]; HSsurf = First[{(130/L <= x <= 170/L && y == 170/L) || (0 <= y <= 170/L && x == 130/L) || (0 <= y <= 170/L && x == 170/L)}]; HS = ImplicitRegion[130/L <= x <= 170/L && 0 <= y <= 170/L, {x, y}]; \[CapitalOmega]3 = RegionDifference[\[CapitalOmega]1, RegionUnion[Flatten@tubes]]; \[CapitalOmega]4 = RegionDifference[\[CapitalOmega]3, HS]; \[CapitalOmega] = RegionDifference[\[CapitalOmega]4, \[CapitalOmega]d]; (*mesh*)f = Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices]; If[0 <= x <= ductW + 3*slab + offsetx1 && y >= 0.4, area > .00001, area > 0.01]]]; mesh = ToElementMesh[\[CapitalOmega], AccuracyGoal -> 5, MeshRefinementFunction -> f /. sizes] mesh["Wireframe"]

parameters = {Pr -> 0.7, 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; walls = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, ts || x == 0 || x == length || y == 0 || y == height || ds || HSsurf]} /. sizes; reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]; bcsflow = {walls, reference}; bcstemperatures = {DirichletCondition[T[t, x, y] == 0, ts],DirichletCondition[T[t, x, y] == 1, HSsurf]}; ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0}; 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]

It takes 226s on my laptop. Visualization

{ContourPlot[temperature[1, x, y], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, ContourStyle -> None, Contours -> 20], Show[DensityPlot[ Norm[{xVel[1, x, y], yVel[1, x, y]}], {x, y} \[Element] mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotRange -> All, PlotPoints -> 100], StreamPlot[{xVel[1, x, y], yVel[1, x, y]}, {x, y} \[Element] mesh, ColorFunction -> None, StreamStyle -> LightGray]]}

This is not much different from the previous case, but this one is easier to manage mesh.

in the first case, a constant heat flux bc was used, which is what I want. But fluid flows freely through the heat source, which is not expected. In the second case, the heat source is a block and fluid cannot flow freely through it, which is a type of heat source that I need. But the bc is a constant temperature, which is supposed to be constant flux instead, e.g. q=10. Is it possible to simulate the heat source as a solid block with a constant heat flux? Thank you. — user55777
– user55777, Commented Jul 24, 2024 at 10:20
In the wall bc of the 1st case, walls={DirichletCondition[{u[t,x,y]==0, v[t,x,y]==0}, ... ||HSsurf]}/.sizes, the surface of the heat source HSsurf is included. However, it doesn't seem to be working. — user55777
– user55777, Commented Jul 24, 2024 at 12:18
@user55777 Please, note, that first case is actually your code in normalized cavity geometry. As I understand from your HS definition there is no constant flux, but raiser HS is a volumetric heat source. This source can easily be converted to a source with a constant surface temperature. To do this, we must calculate the Nusselt number Nu on the surface. Integrate and equate to the volumetric source $\int_S( -\lambda \nabla T).\vec{dS}=q V$. This way we will determine the surface temperature $T_1$. — Alex Trounev
– Alex Trounev, Commented Jul 24, 2024 at 15:38
But if you think that the temperature is distributed unevenly, then you need to set the physical parameters of the heater. — Alex Trounev
– Alex Trounev, Commented Jul 24, 2024 at 15:40
in your answer, either we have HS = With[{q = 10}, If[130/L <= x <= 170/L && 0 <= y <= 170/L, q, 0]] or DirichletCondition[T[t, x, y] == 1, HSsurf] for the heat source, on the other hand we have T[0, x, y] == 0 as an ic, which are inconsistent with each other at t=0. Well, the solver can give plausible solutions. Do you think the solutions have some potential issues, as mentioned by user21? Thank you! — user55777
– user55777, Commented Jul 31, 2024 at 1:54

user21 · Accepted Answer · 2024-07-22 07:30:40Z

The main problem is that your initial condition for T is not consistent with the PDE. You define:

ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0};

and at the same time have:

HS = If[130 <= x <= 170 && 0 <= y <= 170, 10, 0]

This means that in some part of the domain, you have a non-zero T at t=0 with is inconsistent and that's why the solver is making smaller and smaller step sized and is not able to resolve this. You can notice this when you remove the HS from the equation, the time integration will be quick.

The main question you have to answer is if the fluid can flow through the heater area or not. You have two choices: 1) use the approach from Alex's answer where you remove the heater block and model it as a boundary condition (fluid can not flow there). The second (which I will show) is to model the heater block appropriately (fluid can flow in the heater area).

To model the heater block you will need two things. You need to ramp up the heat to make it physical and you will benefit from a geometry that has the heater as a sub-region part. This second requirement is a bit tricky to get and I hope that a future version of Mathematica will make this easier.

The first part is exactly like your, generate the mesh. Now, we extract the region holes from the mesh:

regionHoles = mesh["RegionHoles"] (* {{11., 279.}, {10.5848, 225.308}, {16., 251.}, {11., 251.}, {16., 255.}, {11., 255.}, {16., 259.}, {11., 259.}, {16., 263.}, {11., 263.}, {16., 267.}, {11., 267.}, {16., 271.}, {11., 271.}, {16., 275.}, {11., 275.}, {16., 279.}} *)

Next,

bmesh = ToBoundaryMesh[mesh] ru = ToBoundaryMesh[ RegionUnion[MeshRegion[bmesh], Line[{{130, 0}, {130, 170}, {170, 170}, {170, 0}}]]] (newMesh = ToElementMesh[ru, "RegionHoles" -> regionHoles, MaxCellMeasure -> 100])["Wireframe"]

Note, how the region now has an internal boundary where the heater is.

Here is a ramp function that ramps up d from 0 to 1 in the range from 0 to transition.

ramp[d_] := Evaluate[ With[{transition = 10^-3}, If[d <= transition, d/transition, 1]]]

With this we can define HS2:

HS2 = With[{r = ramp[t]}, If[130 <= x <= 170 && 0 <= y <= 170, r*10, 0]]

And now the time integration will start. This will be computationally more expensive to compute then to remove the heater and model it through a boundary condition.

tmax = 1; Monitor[AbsoluteTiming[{xVel, yVel, pressure, temperature} = NDSolveValue[ Flatten@{op == {0, 0, 0, HS2}, bcsflow, bcstemperatures, ic}, {u, v, p, T}, {x, y} \[Element] newMesh, {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]

One last thing you may want to consider: Double check that the boundary condition of the heater edge is what you want in terms of the heater domain.

Concerning you questions:

Your first goal should be to get a correct result and then you can think about speeding things up
I think it's unlikely that you'll see much or any quality improvement when using a quad dominant mesh.

thank you for the suggestion and comments! If fluid can flow through the heater, should we remove the heat source surface HSsurf from the wall bc, that is, walls = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, ts||x==0||x==length||y==0||y==height||ds(*||HSsurf*)]}/.sizes? — user55777
– user55777, Commented Jul 31, 2024 at 5:19
your method looks resolved the inconsistent bc and ic for the first case of Alex, in which the fluid can flow through the heater area freely. Well, for the 2nd case of Alex, where the the heater was removed and fluid cannot flow there, could this method also resolve the bc&ic inconsistency in this case? Thank you for the help! — user55777
– user55777, Commented Jul 31, 2024 at 7:07
Yes, you most likely can remove the wall DirichletCondition for the heater wall. For the second question, most likely you'd need to ramp up the system from 0 everywhere. — user21
– user21, Commented Aug 1, 2024 at 5:36

Stack Exchange Network

Speed up finite element method solver with small features

2 Answers 2

Linked

Hot Network Questions

Speed up finite element method solver with small features

2 Answers 2

Linked

Related

Hot Network Questions