Speeding up Boussinesq equations solving?

We use scaling as in our model on this page.See also Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.
A system of equations describing free convection can be written in nondimensional form as follows $\nabla .\vec u=0, \frac{d\vec u}{dt}=Pr\nabla ^2\vec u -RaPrT\frac {\vec g}{g}$

$\frac {dT}{dt}=\nabla^2T$

$ d/dt =\partial/\partial t +(\vec u .\nabla )$

$\vec u$ is velocity field vector, $T$ - temperature, $\vec g $ - gravity vector, $Pr$ - the Prandtl number, $Ra$ - the Rayleigh number.

Note, that the temperature scale 40-(-40)=80 is included in Ra definition. Code should be modified as follows

Needs["NDSolve`FEM`"]; ResourceFunction["FEMAddOnsInstall"][]; Needs["FEMAddOns`"] ClearAll[length, height, offsetx1, offsetx2, offsety]; sizes = {length -> 21/10, height -> 15/10, offsetx1 -> 1/10, offsetx2 -> 1/10, offsety1 -> 1/10, offsety2 -> 2/10 }; Geometry \[CapitalOmega]1 = Rectangle[{0, 0}, {length, height}] /. sizes; \[CapitalOmega]2 = Rectangle[{0, 0} + {0.2, 1.5 - 0.049}, {0.69, 0.049} + {0.2, 1.5 - 0.049}]; \[CapitalOmega]D1 = RegionDifference[\[CapitalOmega]1, \[CapitalOmega]2]; \[CapitalOmega]3 = Rectangle[{0, 0} + {1.2, 1.5 - 0.049}, {0.69, 0.049} + {1.2, 1.5 - 0.049}]; \[CapitalOmega]D2 = RegionDifference[\[CapitalOmega]D1, \[CapitalOmega]3]; Mesh \[CapitalOmega] = DiscretizeRegion[\[CapitalOmega]D2, MaxCellMeasure -> 0.001, AccuracyGoal -> 5, PrecisionGoal -> 5]; meshes = ToBoundaryMesh /@ {\[CapitalOmega]1, \[CapitalOmega]2, \ \[CapitalOmega]3}; bmesh = FEMUtils`BoundaryElementMeshJoin @@ meshes; boxCoordinate = {1.05, 0.8}; PCMCoordinate1 = {0.5, 1.46}; PCMCoordinate2 = {1.5, 1.46}; markerColors = {Blue, Orange, Green}; markerCoordinate = {{boxCoordinate}, {PCMCoordinate1}, \ {PCMCoordinate2}}; markerspec = MapIndexed[{First@#1, First@#2} &, markerCoordinate]; mesh = ToElementMesh[bmesh, "RegionMarker" -> markerspec, "MaxCellMeasure" -> 5*10^-3]; mesh["Wireframe"[]] mesh["Wireframe"[ "MeshElementStyle" -> Map[Directive[FaceForm[#], EdgeForm[]] &, markerColors]]] Show[{bmesh["Wireframe"], Graphics[MapThread[{PointSize[0.02], #1, Point /@ #2} &, {markerColors, markerCoordinate}]]}] Material Properties region = <|"air" -> 1, "pcm" -> 2|>; vars = {T[t, x, y], t, {x, y}}; HeatTransferPDEComponent[vars, <|"MassDensity" -> \[Rho], "SpecificHeatCapacity" -> Cp, "ThermalConductivity" -> k|>] airpara = {Subscript[\[Rho], air] -> QuantityMagnitude[ThermodynamicData["Air", "Density"]], Subscript[Cp, air] -> QuantityMagnitude[ ThermodynamicData["Air", "IsobaricHeatCapacity"]], Subscript[k, air] -> QuantityMagnitude[ ThermodynamicData["Air", "ThermalConductivity"]]}; HS30Npara = {Subscript[Cp, pcmS] -> 2.1, Subscript[Cp, pcmL] -> 2.7, Subscript[k, pcmS] -> 2.34, Subscript[k, pcmL] -> 0.6, Subscript[\[Rho], pcmS] -> 1430, latent -> 200, tpmin -> -34, tpmax -> -26 }; PUFpara = {\[Rho] -> 50, cp -> 1500, k -> 0.025}; pars = <||>; pars["MassDensity"] = If[ElementMarker == region["air"], Subscript[\[Rho], air], Subscript[\[Rho], pcmS]]; pars["SpecificHeatCapacity"] = If[ElementMarker == region["air"], Subscript[k, air], Subscript[k, pcmS] ]; pars["ThermalConductivity"] = If[ElementMarker == region["air"], Subscript[k, air], Subscript[k, pcmS]*IdentityMatrix[2]]; op = HeatTransferPDEComponent[vars, pars] Initial and Boundary Condition Tpcm = 0; Tair = 0; ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == Piecewise[{{Tpcm, ElementMarker == region["pcm"]}, {Tair, ElementMarker == region["air"]}}]}; Subscript[\[CapitalGamma], convectiveLeft] = HeatTransferValue[ x == 0 , {T[t, x, y], t, {x, y}}, <||>, <| "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>]; Subscript[\[CapitalGamma], convectiveRight] = HeatTransferValue[ x == length /. sizes, {T[t, x, y], t, {x, y}}, <||>, <| "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>]; Subscript[\[CapitalGamma], convectiveTop] = HeatTransferValue[ y == height /. sizes, {T[t, x, y], t, {x, y}}, <||>, <| "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>]; Subscript[\[CapitalGamma], convectiveBottom] == HeatTransferValue[ y == 0 /. sizes, {T[t, x, y], t, {x, y}}, <||>, <| "HeatTransferCoefficient" -> 0.7, "AmbientTemperature" -> 1|>]; bcstemperatures = { DirichletCondition[ T[t, x, y] == 0, {x, y} \[Element] \[CapitalOmega]2 || {x, y} \[Element] \[CapitalOmega]3]; } ClearAll[\[Nu], \[Epsilon], Pr, Ra]; g = 9.8; \[Alpha] = 0.0034; \[Nu] = 1.48*10^-5; \[Beta] = 2.17*10^-5; L = (4*2.1*1.5)/(2*(2.1 + 1.5)); Temp = 40; Subscript[T, \[Infinity]] = 0; Ra -> ((g*\[Alpha])/(\[Nu]*\[Beta]))*(Temp - Subscript[T, \[Infinity]])*L^3; 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; walls = { DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, {x, y}] }; reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0]; bcsflow = {walls, reference }; 

Solution over mesh with $Ra=10^6$ takes about 120 s

tmax = 1; Monitor[ AbsoluteTiming[ {xVel, yVel, pressure, temperature} = NDSolveValue[ Flatten@{op == {0, 0, 0, Subscript[\[CapitalGamma], convectiveLeft] + Subscript[\[CapitalGamma], convectiveRight] + Subscript[\[CapitalGamma], convectiveTop]}, 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 ] 

Visualization

{StreamDensityPlot[ {xVel[tmax, x, y], yVel[tmax, x, y]}, {x, y} \[Element] xVel["ElementMesh"], AspectRatio -> Automatic, PlotRange -> All, StreamPoints -> Fine, PlotLegends -> Automatic ], DensityPlot[ temperature[tmax, x, y], {x, y} \[Element] temperature["ElementMesh"], AspectRatio -> Automatic, PlotRange -> All, ColorFunction -> "TemperatureMap", PlotLegends -> Automatic ]} 

Solution on $\Omega$ with $Ra=10^5$ takes about 90 s