Robust globally divergence-free weak Galerkin finite element methods for natural convection problems

TL;DR

This paper introduces a new class of weak Galerkin finite element methods for stationary natural convection problems, ensuring divergence-free velocity solutions and demonstrating robustness and optimal error estimates through theoretical analysis and numerical experiments.

Contribution

The paper develops divergence-free weak Galerkin methods with proven stability, optimal error estimates, and an unconditionally convergent algorithm for natural convection problems.

Findings

Methods produce divergence-free velocity fields.

Numerical results confirm robustness across Rayleigh numbers.

Optimal a priori error estimates are established.

Abstract

This paper proposes and analyzes a class of weak Galerkin (WG) finite element methods for stationary natural convection problems in two and three dimensions. We use piecewise polynomials of degrees k, k-1, and k(k>=1) for the velocity, pressure, and temperature approximations in the interior of elements, respectively, and piecewise polynomials of degrees l, k, l(l = k-1,k) for the numerical traces of velocity, pressure and temperature on the interfaces of elements. The methods yield globally divergence-free velocity solutions. Well-posedness of the discrete scheme is established, optimal a priori error estimates are derived, and an unconditionally convergent iteration algorithm is presented. Numerical experiments confirm the theoretical results and show the robustness of the methods with respect to Rayleigh number.

Tables

Table 1. Table 1: Results for different methods with k = 1 𝑘 1 k=1

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	5.9412E-01		1.6959E-01		4.4819E-01		2.4656E-01		2.7341E-02		3.9988E-16
$16\times 8$	3.1494E-01	0.92	4.7778E-02	1.83	2.3637E-01	0.92	1.2464E-01	0.98	6.8747E-03	1.99	1.9062E-15
$32\times 16$	1.5988E-01	0.98	1.2396E-02	1.95	1.1983E-01	0.98	6.2498E-02	0.99	1.7191E-03	2.00	3.0715E-15
$64\times 32$	8.0247E-02	0.99	3.1249E-03	1.99	6.0122E-02	0.99	3.1272E-02	1.00	4.2894E-04	2.00	3.1834E-14
$128\times 64$	4.0162E-02	1.00	7.8018E-04	2.00	3.0087E-02	1.00	1.5639E-02	1.00	1.0704E-04	2.00	4.6475E-14

Table 2. (a) Method:WG-I

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	5.9412E-01		1.6959E-01		4.4819E-01		2.4656E-01		2.7341E-02		3.9988E-16
$16\times 8$	3.1494E-01	0.92	4.7778E-02	1.83	2.3637E-01	0.92	1.2464E-01	0.98	6.8747E-03	1.99	1.9062E-15
$32\times 16$	1.5988E-01	0.98	1.2396E-02	1.95	1.1983E-01	0.98	6.2498E-02	0.99	1.7191E-03	2.00	3.0715E-15
$64\times 32$	8.0247E-02	0.99	3.1249E-03	1.99	6.0122E-02	0.99	3.1272E-02	1.00	4.2894E-04	2.00	3.1834E-14
$128\times 64$	4.0162E-02	1.00	7.8018E-04	2.00	3.0087E-02	1.00	1.5639E-02	1.00	1.0704E-04	2.00	4.6475E-14

Table 3. (b) Method:WG-II

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	7.0486E-01		7.8104E-01		4.7353E-01		2.6104E-01		1.3922E-01		1.5492E-15
$16\times 8$	3.2996E-01	1.10	1.8899E-01	2.05	2.3962E-01	0.98	1.2868E-01	1.02	3.5017E-02	1.99	5.5321E-16
$32\times 16$	1.6192E-01	1.03	4.8031E-02	1.98	1.2025E-01	0.99	6.4066E-02	1.01	8.7749E-03	2.00	6.8348E-15
$64\times 32$	8.0518E-02	1.01	1.2196E-02	1.98	6.0178E-02	1.00	3.1996E-02	1.00	2.1989E-03	2.00	9.5579E-15
$128\times 64$	4.0158E-02	1.00	3.0774E-03	1.99	3.0095E-02	1.00	1.5993E-02	1.00	5.5203E-04	2.00	2.0390E-14

Table 4. (c) Method:WG-III

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	7.4774E-01		8.3792E-01		4.7910E-01		3.1663E-01		1.6162E-01		1.0304E-16
$16\times 8$	3.3583E-01	1.02	1.9985E-01	2.07	2.4031E-01	0.99	1.5503E-01	1.03	4.0626E-02	1.99	2.0466E-16
$32\times 16$	1.6272E-01	1.01	5.0551E-02	1.98	1.2033E-01	1.00	7.7080E-02	1.01	1.0178E-02	2.00	1.7369E-15
$64\times 32$	8.0623E-02	1.00	1.2810E-02	1.98	6.0183E-02	1.00	3.8485E-02	1.00	2.5494E-03	2.00	2.3551E-15
$128\times 64$	4.0212E-02	1.00	3.2282E-03	1.99	3.0093E-02	1.00	1.9235E-02	1.00	6.3946E-04	2.00	6.5944E-15

Table 5. Table 2: Results for different methods with k = 2 𝑘 2 k=2

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	1.6192E-01		2.8177E-02		6.6611E-02		2.3814E-02		1.5210E-03		1.5852E-15
$16\times 8$	4.2800E-02	1.92	3.6801E-03	2.94	1.7476E-02	1.92	5.9899E-03	1.98	1.9029E-04	2.99	1.3501E-14
$32\times 16$	1.0767E-02	1.99	4.6124E-04	2.99	4.4430E-03	1.98	1.4995E-03	1.99	2.3790E-05	3.00	6.8867E-14
$64\times 32$	2.6808E-03	2.01	5.7386E-05	3.01	1.1115E-03	1.99	3.7495E-04	2.00	2.9736E-06	3.00	3.8064E-14
$128\times 64$	6.7021E-04	2.00	7.1513E-06	3.00	2.7795E-04	2.00	9.3738E-05	2.00	3.7173E-07	3.00	7.2047E-14

Table 6. (a) Method:WG-I

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	1.6192E-01		2.8177E-02		6.6611E-02		2.3814E-02		1.5210E-03		1.5852E-15
$16\times 8$	4.2800E-02	1.92	3.6801E-03	2.94	1.7476E-02	1.92	5.9899E-03	1.98	1.9029E-04	2.99	1.3501E-14
$32\times 16$	1.0767E-02	1.99	4.6124E-04	2.99	4.4430E-03	1.98	1.4995E-03	1.99	2.3790E-05	3.00	6.8867E-14
$64\times 32$	2.6808E-03	2.01	5.7386E-05	3.01	1.1115E-03	1.99	3.7495E-04	2.00	2.9736E-06	3.00	3.8064E-14
$128\times 64$	6.7021E-04	2.00	7.1513E-06	3.00	2.7795E-04	2.00	9.3738E-05	2.00	3.7173E-07	3.00	7.2047E-14

Table 7. (b) Method:WG-II

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	2.5023E-01		5.9209E-02		6.6212E-02		4.1197E-02		4.9610E-03		6.5550E-16
$16\times 8$	6.3163E-02	1.98	7.4474E-03	2.99	1.7485E-02	1.92	1.0276E-02	1.99	6.1111E-04	3.02	7.4872E-15
$32\times 16$	1.5659E-02	2.01	9.3395E-04	2.99	4.4432E-03	1.98	2.5691E-03	2.00	7.5883E-05	3.01	5.6488E-15
$64\times 32$	3.8820E-03	2.01	1.1720E-04	3.00	1.1117E-03	2.00	6.4257E-04	2.00	9.4569E-06	3.00	2.4648E-14
$128\times 64$	9.6547E-04	2.00	1.4685E-05	3.00	2.7815E-04	2.00	1.6070E-04	2.00	1.1804E-06	3.00	2.1412E-13

Table 8. (c) Method:WG-III

mesh	$\frac{{\parallel \nabla 𝒖-{\nabla }_h{𝒖}_{h0}\parallel }_0}{{\parallel \nabla 𝒖\parallel }_0}$		$\frac{{\parallel 𝒖-{𝒖}_{h0}\parallel }_0}{{\parallel 𝒖\parallel }_0}$		$\frac{{\parallel p-p_{h0}\parallel }_0}{{\parallel p\parallel }_0}$		$\frac{{\parallel \nabla T-{\nabla }_h{T}_{h0}\parallel }_0}{{\parallel \nabla T\parallel }_0}$		$\frac{{\parallel T-T_{h0}\parallel }_0}{{\parallel T\parallel }_0}$		$di{v}_h{U}_h$
	error	order	error	order	error	order	error	order	error	order
$8\times 4$	1.3075E-01		6.2217E-02		6.6237E-02		2.1605E-02		5.3332E-03		1.6739E-15
$16\times 8$	3.4979E-02	1.90	7.6750E-03	3.02	1.7492E-02	1.92	5.4667E-03	1.98	6.6033E-04	3.02	3.5685E-15
$32\times 16$	8.9627E-03	1.96	9.4948E-04	3.01	4.4333E-03	1.98	1.3734E-03	1.99	8.2232E-05	3.01	3.2603E-14
$64\times 32$	2.2617E-03	1.99	1.1834E-04	3.00	1.1121E-03	2.00	3.4409E-04	2.00	1.0263E-05	3.00	7.3185E-14
$128\times 64$	5.6761E-04	2.00	1.4791E-05	3.00	2.7826E-04	2.00	8.6112E-05	2.00	1.2821E-06	3.00	2.6642E-13

$u{1}_{max}$	the maximum horizontal velocity on the vertical mid-plane of the cavity
$u{2}_{max}$	the maximum vertical velocity on the horizontal mid-plane of the cavity
$\overline{Nu}$	the average Nusselt number throughout the cavity
$N{u}_{max}$	the maximum value of the local Nusselt number on the boundary at x=0
$N{u}_{min}$	the minimum value of the local Nusselt number on the boundary at x=0

Table 10. Table 3: Natural convection in a square cavity: comparison with the benchmark solutions

Ra		WG-I,k=1	WG-I,k=1	WG-I,k=2	WG-I,k=2	Ref.[12]	Ref.[38]	Ref.[21]	Ref.[22]	Ref.[29]
		$40\times 40$	$70\times 70$	$40\times 40$	$50\times 50$		$64\times 64$	$70\times 70$	$70\times 70$	$100\times 100$
${10}^3$	$u{1}_{max}$	3.653	3.654	3.640	3.646	3.649	-	3.68	-	3.489
	$u{2}_{max}$	3.711	3.698	3.697	3.697	3.697	-	3.73	3.686	3.69
	$\overline{Nu}$	1.118	1.118	1.118	1.118	1.118	-	1.074	1.117	1.117
	$N{u}_{max}$	1.506	1.506	1.560	1.506	1.505	-	1.47	-	1.501
	$N{u}_{min}$	0.691	0.691	0.691	0.691	0.692	-	0.623	-	0.691
${10}^4$	$u{1}_{max}$	16.227	16.188	16.183	16.180	16.178	16.19	16.10	-	16.122
	$u{2}_{max}$	19.744	19.611	19.600	16.628	19.617	19.63	19.90	19.63	19.79
	$\overline{Nu}$	2.243	2.244	2.245	2.245	2.243	-	2.084	2.243	2.254
	$N{u}_{max}$	3.528	3.530	3.531	3.531	3.528	-	3.47	-	3.579
	$N{u}_{min}$	0.585	0.585	0.585	0.585	0.586	-	0.4968	-	0.577
${10}^5$	$u{1}_{max}$	34.829	34.771	34.715	34.702	34.81	34.74	34.00	-	34.00
	$u{2}_{max}$	69.049	68.736	67.875	68.290	68.22	68.48	70.00	68.85	70.63
	$\overline{Nu}$	4.515	4.519	4.522	4.522	4.519	-	4.30	4.521	4.598
	$N{u}_{max}$	7.701	7.713	7.716	7.720	7.717	-	7.71	-	7.945
	$N{u}_{min}$	0.726	0.727	0.728	0.728	0.729	-	0.614	-	0.698
${10}^6$	$u{1}_{max}$	64.977	64.710	64.835	64.541	64.63	64.81	65.40	-	65.40
	$u{2}_{max}$	217.307	221.534	208.237	220.609	219.36	220.46	228	221.6	227.11
	$\overline{Nu}$	8.797	8.813	8.825	8.825	8.800	-	8.743	8.806	8.976
	$N{u}_{max}$	17.676	17.511	17.462	17.536	17.925	-	17.46	-	17.86
	$N{u}_{min}$	0.970	0.976	0.980	0.979	0.989	-	0.716	-	0.913
${10}^7$	$u{1}_{max}$	154.770	148.802	148.454	148.596	145.267*	148.40	139.7	-	143.56
	$u{2}_{max}$	819.329	695.512	703.702	707.696	703.253*	694.14	698	702.3	714.48
	$\overline{Nu}$	16.564	16.484	16.522	16.521	-	-	13.99	16.40	16.656
	$N{u}_{max}$	47.155	40.374	40.935	40.329	41.025*	-	30.46	-	38.6
	$N{u}_{min}$	1.359	1.353	1.363	1.367	1.380*	-	0.787	-	1.298

Equations

$$\displaystyle\left\{\begin{array}[]{rl}-\Pr\Delta\bm{u}+\nabla\cdot(\bm{u}\otimes\bm{u})+\nabla p-\Pr Ra\bm{j}T=\bm{f}&in\quad\Omega_{f},\\ \nabla\cdot\bm{u}=0&in\quad\Omega_{f},\\ -\kappa\Delta T+\nabla\cdot(\bm{u}T)=g&in\quad\Omega,\\ \bm{u}\equiv\bm{0}&in\quad\Omega_{s}\bigcup\partial\Omega_{f},\\ T=0&on\quad\partial\Omega,\end{array}\right.$$

$$⟨u,v{⟩}_{\partial {\mathcal{T}}_h}:=$$

$$(u,v{)}_{{\mathcal{T}}_h}:=$$

$$\mathbf{W}:=\{\mathbf{v}\in [H_0^1({\Omega }_f){]}^d:\nabla \cdot \mathbf{v}=0\}$$

$$a(\mathbf{u},\mathbf{v})$$

$$d(T,\mathbf{v})$$

$$c(\mathbf{w};\mathbf{u},\mathbf{v})$$

$$c(\mathbf{u};\mathbf{u},\mathbf{v})$$

$$\overline{c}(\mathbf{u};T,s)$$

$$\displaystyle\left\{\begin{array}[]{rl}A(\bm{u};\bm{u},\bm{v})+b(\bm{v},p)-b(\bm{u},q)-d(T,\bm{v})&=(\bm{f},\bm{v}),\forall(\bm{v},q)\in\bm{W}\times L_{0}^{2}(\Omega_{f}),\\ \overline{A}(\bm{u};T,s)&=(g,s),\forall s\in H_{0}^{1}(\Omega),\end{array}\right.$$

$$A(\mathbf{u};\mathbf{u},\mathbf{v})$$

$$\overline{A}(\mathbf{u};T,s)$$

$$({\mathrm{Pr}}^{-1}\mathcal{N}Ra{\kappa }^{-1}+Ra\mathcal{M}{\kappa }^{-2})\Vert g{\Vert }_{-1}+\mathcal{N}{\mathrm{Pr}}^{-2}\Vert \mathbf{f}{\Vert }_{-1}<1,$$

$$\mathcal{N}:=\underset{0\ne \mathbf{w},\mathbf{u},\mathbf{v}\in \mathbf{W}}{\sup }\frac{c(\mathbf{w};\mathbf{u},\mathbf{v})}{|\mathbf{w}{|}_1|\mathbf{u}{|}_1|\mathbf{v}{|}_1},\mathcal{M}:=\underset{\begin{array}{c}0\ne \mathbf{w}\in \mathbf{W},\\ 0\ne T,s\in H_0^1(\Omega )\end{array}}{\sup }\frac{\overline{c}(\mathbf{w};T,s)}{|\mathbf{w}{|}_1|T{|}_1|s{|}_1}.$$

$$({\mathrm{Pr}}^{-1}\mathcal{N}Ra{\kappa }^{-1}+Ra\mathcal{M}{\kappa }^{-2})\Vert g{\Vert }_{-1}+\mathcal{N}{\mathrm{Pr}}^{-2}\Vert \mathbf{f}{\Vert }_{-1}<1-\delta .$$

$$({\nabla }_{w,r,K}v,\mathbf{\tau }{)}_K=-(v_0,\nabla \cdot \mathbf{\tau }{)}_K+⟨v_b,\mathbf{\tau }\cdot {\mathbf{n}}_K{⟩}_{\partial K}\forall \mathbf{\tau }\in [P_r(K){]}^d.$$

$${\nabla }_{w,r}\mathbf{v}:=({\nabla }_{w,r}{v}_1,\cdots ,{\nabla }_{w,r}{v}_d{)}^T.$$

$$({\nabla }_{w,r,K}\cdot \mathbf{v},\tau {)}_K=-({\mathbf{v}}_0,\nabla \tau {)}_K+⟨{\mathbf{v}}_b\cdot {\mathbf{n}}_K,\tau {⟩}_{\partial K}\forall \tau \in P_r(K).$$

$${\nabla }_{w,r}\cdot {|}_K={\nabla }_{w,r,K}\cdot ,\forall K\in {\mathcal{T}}_h.$$

$${\nabla }_{w,r}\cdot \mathbf{w}=({\nabla }_{w,r}\cdot {\mathbf{w}}_1,\cdots ,{\nabla }_{w,r}\cdot {\mathbf{w}}_d{)}^T.$$

$${\mathbf{V}}_h$$

$${\mathbf{V}}_h^0$$

$$Q_h$$

$$Q_h^0$$

$$S_h$$

$$S_h^0$$

$$a_h({\mathbf{u}}_h,{\mathbf{v}}_h)$$

$$b_h({\mathbf{v}}_h,q_h)$$

$$d_h(T_h,{\mathbf{v}}_h)$$

$${\overline{a}}_h(T_h,s_h)$$

$$c_h({\mathbf{w}}_h;{\mathbf{u}}_h,{\mathbf{v}}_h)$$

$${\overline{c}}_h({\mathbf{u}}_h;T_h,s_h)$$

$$c_h({\mathbf{w}}_h;{\mathbf{v}}_h,{\mathbf{v}}_h)=0,{\overline{c}}_h({\mathbf{w}}_h;s_h,s_h)=0.$$

$$\displaystyle\left\{\begin{array}[]{rl}A_{h}(\bm{u}_{h};\bm{u}_{h},\bm{v}_{h})+b_{h}(\bm{v}_{h},p_{h})-b_{h}(\bm{u}_{h},q_{h})-d_{h}(T_{h},\bm{v}_{h})&=(\bm{f},\bm{v}_{h0}),\forall(\bm{v}_{h},q_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0},\\ \overline{A}_{h}(\bm{u}_{h};T_{h},s_{h})&=(g,s_{h0}),\forall s_{h}\in S_{h}^{0},\end{array}\right.$$

$$A_h({\mathbf{w}}_h;{\mathbf{u}}_h,{\mathbf{v}}_h)$$

$${\overline{A}}_h({\mathbf{u}}_h;T_h,s_h)$$

$$r_{hb}|_{e}=\left\{\begin{array}[]{l}-(\bm{u}_{h0}\cdot\bm{n}_{e})|_{K_{1}\bigcap e}-(\bm{u}_{h0}\cdot\bm{n}_{e})|_{K_{2}\bigcap e},\quad\forall e\in\varepsilon_{h}^{f}/\partial\Omega_{f},\\ 0,\quad\forall e\in\partial\Omega_{f},\end{array}\right.$$

$$\Vert {\nabla }_h\cdot {\mathbf{u}}_{h0}{\Vert }_0^2+\sum _{e\in {\epsilon }_h^f/\partial {\Omega }_f}\Vert ({\mathbf{u}}_{h0}\cdot {\mathbf{n}}_e){|}_{K_1}+({\mathbf{u}}_{h0}\cdot {\mathbf{n}}_e){|}_{K_2}{\Vert }_{0,e}^2=0.$$

Full text

\ams

52B10, 65D18, 68U05, 68U07

1 Introduction

Let $\mathbb{R}^{d}(d=2,3)$ be a polygonal or polyhedral domain with a polygonal or polyhedral subdomain $\Omega_{f}\subset\Omega$ and $\Omega_{s}:=\Omega\setminus\Omega_{f}$, we consider the following stationary natural convection (or conduction-convection) problem: seek the velocity $\bm{u}=(u_{1},u_{2},\cdots,u_{d})^{T}$, the pressure $p$, and the temperature $T$ such that

[TABLE]

where $\otimes$ is defined by $\bm{u}\otimes\bm{v}=(u_{i}v_{j})_{d\times d}$ for $\bm{v}=(v_{1},v_{2},\cdots,v_{d})^{T}$, $\bm{j}$ is the vector of gravitational acceleration with $\bm{j}=(0,1)^{T}$ when $d=2$ and $\bm{j}=(0,0,1)^{T}$ when $d=3$, $\bm{f}\in[L^{2}(\Omega_{f})]^{d}$, $g\in L^{2}(\Omega)$ are the forcing functions, and $\Pr$, $Ra$ denote the Prandtl and Rayleigh numbers, respectively,.

The model problem (1.6), arising both in nature and in engineering applications, is a coupled system of fluid flow, governed by the incompressible Navier-Stokes equations, and heat transfer, governed by the energy equation. Due to its practical significance, the development of efficient numerical methods for natural convection has attracted a great many of research efforts; see, e.g. [1],[2],[3],[29],[23],[12],[15],[18],[20],[19],[21],[22],[25],[26],[28],[38],[39]. In [2, 3], error estimates for some finite element methods were derived in approximating stationary and non-stationary natural convection problems. [20, 19] applied Petrov-Galerkin least squares mixed finite element methods to discretize the problems. [25, 26] developed a nonconforming mixed element method and a Petrov-Galerkin least squares nonconforming mixed element method for the stationary problems. In [37], three kinds of decoupled two level finite element methods were presented. [38, 39] applied the variational multiscale method to solve the stationary and non-stationary problems.

In this paper, we consider a weak Galerkin (WG) finite element discretization of the model problem (1.6). The WG method was first proposed and analyzed to solve second-order elliptic problems [30, 31]. It is designed by using a weakly defined gradient operator over functions with discontinuity, and then allows the use of totally discontinuous functions in the finite element procedure. Similar to the hybridized discontinuous Galerkin (HDG) method [11], the WG method is of the property of local elimination of unknowns defined in the interior of elements. We note that in some special cases the WG method and the HDG method are equivalent (cf. [6, 7, 8]). In [6], a class of robust globally divergence-free weak Galerkin methods for Stokes equations were developed, and then were extended in [40] to solve incompressible quasi-Newtonian Stokes equations.We also refer to [9, 16, 13, 17, 24, 33, 32, 41, 10, 35, 34, 36] for some other developments and applications of the WG method.

This paper aims to propose a class of WG methods for the natural convection problems. The methods include as unknowns the velocity, pressure, and temperature variables both in the interior of elements and on the interfaces of elements. In the interior of elements, we use piecewise polynomials of degrees $k,k-1,$ and $k$ $(k\geq 1)$ for the velocity, pressure, and temperature approximations, respectively. On the interfaces of elements, we use piecewise polynomials of degrees $l,k,l$ $(l=k-1,k)$ for the numerical traces of velocity, pressure and temperature. The methods are shown to yield globally divergence-free velocity approximations.

The rest of the paper is organized as follows. Section 2 introduces the WG finite element scheme. Section 3 shows the existence and uniqueness of the discrete solution. Section 4 derives a priori error estimates. Section 5 discusses the local elimination property and the convergence of an iteration method for the WG scheme. Finally, Section 6 provides numerical examples to verify the theoretical results.

Throughout this paper, we use $a\apprle b$ $(a\apprge b)$ to denote $a\leq Cb$ $(a\geq Cb)$, where the constant C is positive independent of mesh size $h,h_{K},h_{e}$ and the $\Pr$, $\kappa$ and Rayleigh number.

2 WG finite element scheme

2.1 Notation

For any bounded domain $D\in\mathbb{R}^{s}(s=d,d-1)$, let $H^{m}(D)$ and $H_{0}^{m}(D)$ denote the usual $m^{th}$-order Sobolev spaces on D, and $\lVert\cdot\rVert_{m,D},\lvert\cdot\rvert_{m,D}$ denote the norm and semi-norm on these spaces. We use $(\cdot,\cdot)_{m,D}$ to denote the inner product of $H^{m}(D)$, with $(\cdot,\cdot)_{D}:=(\cdot,\cdot)_{0,D}$. When $D=\Omega$, we set $\lVert\cdot\rVert_{m}:=\lVert\cdot\rVert_{m,\Omega},\lvert\cdot\rvert_{m}:=\lvert\cdot\rvert_{m,\Omega},$ and $(\cdot,\cdot):=(\cdot,\cdot)_{\Omega}$. In particular, when $D\subset R^{d-1}$, we use $\langle\cdot,\cdot\rangle_{D}$ to replace $(\cdot,\cdot)_{D}$. For integer $k\geqslant 0$, $P_{k}(D)$ denotes the set of all polynomials on D with degree no more than $k$. We also need the following spaces:

$L_{0}^{2}(\Omega):=\{q\in L^{2}(\Omega):(q,1)=0\}$,

$\bm{H}(div,D):=\{\bm{v}\in\bm{L}^{2}(D):\nabla\cdot\bm{v}\in L^{2}(D)\}.$

Let $\mathcal{T}_{h}^{s}$ and $\mathcal{T}_{h}^{f}$ be shape-regular simplicial decompositions of the subdomains $\Omega_{s}$ and $\Omega_{f}$, respectively. Then $\mathcal{T}_{h}:=\mathcal{T}_{h}^{s}\cup\mathcal{T}_{h}^{f}=\bigcup\{K\}$ is a shape-regular simplicial decomposition of $\Omega$. Let $\varepsilon_{h}^{s}$ and $\varepsilon_{h}^{f}$ be the sets of all edges (faces) of all elements in $\mathcal{T}_{h}^{s}$ and $\mathcal{T}_{h}^{f}$, respectively, and set $\varepsilon_{h}:=\varepsilon_{h}^{s}\cup\varepsilon_{h}^{f}=\bigcup\{e\}$. For any $K\in\mathcal{T}_{h}$, $e\in\varepsilon_{h}$, we denote by $h_{K}$ and $h_{e}$ the diameters of $K$ and $e$, respectively, and set $h:=\max\limits_{K\in\mathcal{T}_{h}}h_{K}$. Let $\bm{n}_{K}$ and $\bm{n}_{e}$ be the outward unit normal vectors along the boundary $\partial K$ and $e$. We denote by $\nabla_{h}$ and $\nabla_{h}\cdot$ the piecewise-defined gradient and divergence with respect to $\mathcal{T}_{h}$. We also introduce the mesh-dependent inner products and mesh-dependent norms:

[TABLE]

2.2 Weak problem

We first introduce the space

[TABLE]

and the following bilinear and trilinear forms: for any $\bm{u},\bm{v},\bm{w}\in H_{0}^{1}(\Omega_{f})$, $q\in L_{0}^{2}(\Omega_{f})$, and $T,s\in H_{0}^{1}(\Omega)$,

[TABLE]

It is easy to see that, for $\bm{u}\in\bm{W}$,

[TABLE]

Then the variational problem of (1.6) reads as follows: seek $(\bm{u},p,T)\in\bm{W}\times L_{0}^{2}(\Omega_{f})\times H_{0}^{1}(\Omega)$ such that

[TABLE]

where

[TABLE]

Theorem 2.1**.**

[3*]*For $\bm{f}\in[H^{-1}(\Omega_{f})]^{d}$ and $g\in H^{-1}(\Omega)$, the weak problem (2.3) has at least one solution. In addition, it admits a unique solution $(\bm{u},p,T)\in\bm{W}\times L_{0}^{2}(\Omega_{f})\times H_{0}^{1}(\Omega)$ if

[TABLE]

where

[TABLE]

In what follows, we assume that the solution $(\bm{u},p,T)$ is unique and, more precisely, there exists a fixed constant $\delta>0$ such that

[TABLE]

2.3 Discrete weak operators

In order to design a WG finite element scheme for the problem (1.1), we introduce the discrete weak gradient operator $\nabla_{w,r}$ and the discrete weak divergence operator $\nabla_{w,r}\cdot$ as follows.

Definition 2.1**.**

For any $K\in\mathcal{T}_{h}$ and $v\in\mathcal{V}(K):=\left\{\{v_{0},v_{b}\}:v_{0}\in L^{2}(K),v_{b}\in H^{1/2}(\partial K)\right\}$, the discrete weak gradient $\nabla_{w,r,K}v\in[P_{r}(K)]^{d}$ on $K$ is determined by the equation

[TABLE]

Then we define the global discrete weak gradient operator $\nabla_{w,r}$ by

$\nabla_{w,r}|_{K}=\nabla_{w,r,K},\quad\forall K\in\mathcal{T}_{h}$.

For a vector $\bm{v}=(v_{1},\cdots,v_{d})^{T}\in[\mathcal{V}(K)]^{d}$, we define its discrete weak gradient $\nabla_{w,r}\bm{v}$ by

[TABLE]

Definition 2.2**.**

For any $K\in\mathcal{T}_{h}$ and $\bm{v}\in\mathcal{W}(K):=\left\{\{\bm{v}_{0},\bm{v}_{b}\}:\bm{v}_{0}\in[L^{2}(K)]^{d},\bm{v}_{b}\cdot\bm{n}_{K}\in H^{-1/2}(\partial K)\right\}$, the discrete weak divergence $\nabla_{w,r,K}\cdot\bm{v}\in P_{r}(K)$ is determined by the equation

[TABLE]

Then we define the global discrete weak divergence operator $\nabla_{w,r}\cdot$ by

[TABLE]

For a tensor $\bm{w}=(\bm{w}_{1},\cdots,\bm{w}_{d})^{T}\in[\mathcal{W}(K)]^{d\times d}$ with $\bm{w}_{i}\in[\mathcal{W}(K)]^{d}$ for $i=1,\cdots,d$, we define its discrete weak divergence $\nabla_{w,r}\cdot\bm{w}$ by

[TABLE]

2.4 WG finite element scheme

For any $K\in\mathcal{T}_{h},e\in\varepsilon_{h}$ and any integer $j\geq 0$, let $Q_{j}^{0}:L^{2}(K)\rightarrow P_{j}(K)$ and $Q^{b}_{j}:L^{2}(e)\rightarrow P_{j}(e)$ be the usual $L^{2}$ projection operators. We shall use $\bm{Q}_{j}^{b}$ to denote $Q^{b}_{j}$ for vector spaces.

For any integer $k\geq 1$ and $l=k-1,k$, we introduce the following finite dimensional spaces:

[TABLE]

For any $\bm{u}_{h}=\{\bm{u}_{h0},\bm{u}_{hb}\},\bm{v}_{h}=\{\bm{v}_{h0},\bm{v}_{hb}\}\in\bm{V}_{h}^{0}$, $q_{h}=\{q_{h0},q_{hb}\}\in Q_{h}^{0}$, and $T_{h}=\{T_{h0},T_{hb}\},s_{h}=\{s_{h0},s_{hb}\}\in S_{h}^{0}$, define the following bilinear and trilinear forms:

[TABLE]

It is easy to see that

[TABLE]

The WG finite element scheme for (1.6) is then given as follows: seek $\bm{u}_{h}=\{\bm{u}_{h0},\bm{u}_{hb}\}\in\bm{V}_{h}^{0}$, $p_{h}=\{p_{h0},p_{hb}\}\in Q_{h}^{0}$, and $T_{h}=\{T_{h0},T_{hb}\}\in S_{h}^{0}$ such that

[TABLE]

where

[TABLE]

$\tau|_{\partial T}=h_{T}^{-1}$, and m is an integer with $k-1\leq m\leq l$.

Remark 2.1**.**

It’s easy to show that the scheme (2.7) yields globally divergence-free velocity approximation $\bm{u}_{h0}$. In fact, let $K_{1},K_{2}\in\mathcal{T}_{h}$ be any two adjacent elements with a common face $e$, introduce a function $r_{hb}\in L^{2}(\varepsilon_{h})$ with

[TABLE]

and set $c_{0}:=\frac{1}{|\Omega_{f}|}\int_{\Omega_{f}}\nabla_{h}\cdot\bm{u}_{h0}d\bm{x}$. Then, taking $(\bm{v}_{h0},\bm{v}_{hb},q_{h0},q_{hb},s_{h0},s_{hb})=(\bm{0},\bm{0},\nabla_{h}\cdot\bm{u}_{h0}-c_{0},r_{hb}-c_{0},0,0)$ in (2.7) yields

[TABLE]

This indicates $\bm{u}_{h0}\in\bm{H}(div,\Omega_{f})$ and $\nabla_{h}\cdot\bm{u}_{h0}=\nabla\cdot\bm{u}_{h0}=0$, i.e. the velocity approximation $\bm{u}_{h0}$ is globally divergence-free in a pointwise sense.

3 Well-posedness of the discrete scheme

3.1 Some basic results

For the projections $Q_{j}^{0}$ and $Q_{j}^{b}$ with $j\geq 0$, the following stability and approximation results are standard.

Lemma 3.1**.**

*([27])

Let $s$ be an integer with $1\leq s\leq j+1$. Then we have, for any $K\in\mathcal{T}_{h}$ and $e\in\varepsilon_{h}$,*

[TABLE]

By using the trace theorem, the inverse inequality, and scaling arguments metioned in [27], we can get the following lemma.

Lemma 3.2**.**

For all $K\in\mathcal{T}_{h}$, $w\in H^{1}(K)$, and $1\leq\tilde{q}\leq\infty$, we have

[TABLE]

In particular, for all $w\in P_{k}(K)$,

[TABLE]

Lemma 3.3**.**

([6]) Let $0\leq k-1\leq m\leq l\leq k$. For all $K\in\mathcal{T}_{h}$ and $\bm{v}_{h}=\{\bm{v}_{h0},\bm{v}_{hb}\}\in[P_{k}(K)]^{d}\times[P_{l}(\partial K)]^{d}$, the following estimates hold:

[TABLE]

We introduce the following semi-norms: for any $(\bm{v}_{h},q_{h},s_{h})\in\bm{V}_{h}\times Q_{h}\times S_{h}$,

[TABLE]

Here we recall that $\tau|_{\partial K}=h_{K}^{-1}$. It is easy to see that the above three semi-norms are norms on $\bm{V}_{h}^{0}$, $Q_{h}^{0}$ and $S_{h}^{0}$, respectively (cf. [6]). In addition, from the lemma above it follows

[TABLE]

Remark 3.1**.**

We note that the estimates (3.1), (3.2), and (3.3) also hold for all $s_{h}\in S_{h}^{0}$ due to the fact that $s_{h}|_{K}=\{s_{h0},s_{hb}\}\in P_{k}(K)\times P_{l}(\partial K)$.

Lemma 3.4**.**

([14]) For all $s_{h0}\in S_{h0}=\{s_{h0}:s_{h0}|_{K}\in P_{k}(K),\forall K\in\mathcal{T}_{h}\}$, there exists an interpolation $I_{k}s_{h0}\in S_{h0}\cap H_{0}^{1}(\Omega)$ such that

[TABLE]

From this lemma it follows that, for all $\bm{v}_{h0}\in\bm{V}_{h0}=\{\bm{v}_{h0}:\bm{v}_{h0}|_{K}\in P_{k}(K),\forall K\in\mathcal{T}_{h}^{f}\}$, there exists an interpolation $\bm{I}_{k}\bm{v}_{h0}\in\bm{V}_{h0}\cap[H_{0}^{1}(\Omega)]^{d}$ such that

[TABLE]

Lemma 3.5**.**

For all $\bm{v}_{h}\in\bm{V}_{h}^{0}$ and $s_{h}\in S_{h}^{0}$, we have

[TABLE]

where $2\leq\tilde{q}<\infty$ when $d=2$, $2\leq\tilde{q}\leq 6$ when $d=3$, and $C_{\tilde{q}1}$, $C_{\tilde{q}2}$ are positive constants only depending on $\tilde{q}$.

Proof.

For all $\bm{v}_{h}\in\bm{V}_{h}^{0}$, we apply the Sobolev embedding theorem and Poincáre inequality to get

[TABLE]

From (3.5), (3.3), the definition of $\interleave\cdot\interleave$, and the projection property of $\bm{Q}_{l}^{b}$, it follows

[TABLE]

Using the Sobolev embedding theorem and the inverse inequality once again, by the properties of the projection-mean operator ([27]) $\bm{\Pi_{h}}:\bm{V}_{h0}=\{\bm{v}_{h0}:\bm{v}_{h0}|_{K}\in P_{k}(K),\forall K\in\mathcal{T}_{h}^{f}\}\rightarrow W^{1,2}(\Omega_{f})\cap W^{0,\tilde{q}}(\Omega_{f})$ and the fact that $2\leq\tilde{q}<\infty$ when $d=2$ and $2\leq\tilde{q}\leq 6$ when $d=3$, we have

[TABLE]

which, together with (3.8) and (3.9), yields the desired estimate (3.6).

Similarly, we can obtain (3.7). This finishes the proof. ∎

For any nonnegative integer $j$ and any $K\in\mathcal{T}_{h}$, we introduce the local Raviart-Thomas(RT) element space

[TABLE]

Lemmas 3.6-3.8 show some properties of the $RT$ projection which can be founded in ([5].Page 9-10).

Lemma 3.6**.**

For any $\bm{v}_{h0}\in\bm{RT}_{j}(K)$, $\nabla\cdot\bm{v}_{h0}|_{K}=0$ implies $\bm{v}_{h0}\in[P_{j}(K)]^{d}$.

Lemma 3.7**.**

For any $K\in\mathcal{T}_{h}$ and $\bm{v}\in[H^{1}(K)]^{d}$, there exists a unique $\bm{P}_{j}^{RT}\bm{v}\in\bm{RT}_{j}(K)$ such that

[TABLE]

If $j=0$, $\bm{P}_{j}^{RT}\bm{v}$ is determined only by (3.10). Moreover, the following approximation holds:

[TABLE]

Lemma 3.8**.**

The operator $\bm{P}_{j}^{RT}$ defined in Lemma 3.7 satisfies

[TABLE]

Lemma 3.9**.**

([6]) It holds the following commutativity properties:

[TABLE]

3.2 Stability conditions

Lemma 3.10**.**

For any $\bm{u}_{h},\bm{v}_{h}\in\bm{V}_{h}$, and $T_{h},s_{h}\in S_{h}$, the following inequalities hold:

[TABLE]

Proof.

From the definitions of $a_{h}(\cdot,\cdot),\overline{a}_{h}(\cdot,\cdot),c_{h}(\cdot;\cdot,\cdot),\overline{c}_{h}(\cdot;\cdot,\cdot),d_{h}(\cdot,\cdot)$, Cauchy-Schwarz inequality and Lemma 3.5, we can easily get (3.15),(3.17), and (3.21).

For all $\bm{u}_{h},\bm{v}_{h}\in\bm{V}_{h}$, by the definition of $\nabla_{w,k}\cdot$ we have

[TABLE]

In light of Hölder’s inequality and Lemma 3.5, we obtain

[TABLE]

From Hölder’s inequality, Lemma 3.2, Lemma 3.5, and the inverse inequality, it follows

[TABLE]

and

[TABLE]

Similarly, we can get

[TABLE]

As a result, the estimate (3.19) holds.

The estimate (3.20) follows similarly. ∎

By (2.4), Lemma 3.10, and the definitions of the trilinear forms $A_{h}(\cdot;\cdot,\cdot)$ and $\overline{A}_{h}(\cdot;\cdot,\cdot)$, we easily get the following continuity and coercivity results.

Lemma 3.11**.**

For any $\bm{w}_{h},\bm{u}_{h},\bm{v}_{h}\in\bm{V}_{h},T_{h},s_{h}\in S_{h}$, it holds

[TABLE]

By following the same routine as in the proof of ([6, Theorem 3.1]), we can obtain the following inf-sup inequality.

Lemma 3.12**.**

For any $(\bm{v}_{h},q_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}$, it holds

[TABLE]

3.3 Existence and uniqueness results

We define a space

[TABLE]

and introduce the following discretization problem: seek $(\bm{u_{h}},T_{h})\in\bm{W}_{h}\times S_{h}^{0}$

[TABLE]

It is easy to see that, by Lemma 3.12 and the theory of mixed finite element methods [5], the following conclusion holds.

Lemma 3.13**.**

The problems (2.7) and (3.28) are equivalent in the sense that (i) and (ii) hold:

(i) if $(\bm{u}_{h},p_{h},T_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}\times S_{h}^{0}$ is the solution to the problem (2.7), then $(\bm{u_{h}},T_{h})$ is the solution to the problem (3.28);

(ii) if $(\bm{u_{h}},T_{h})\in\bm{W}_{h}\times S_{h}^{0}$ is the solution to the problem (3.28), then $(\bm{u}_{h},p_{h},T_{h})$ is the solution to the problem (2.7), where $p_{h}\in Q_{h}^{0}$ is determined by

[TABLE]

In what follows we shall discuss the existence and uniqueness of the solution to the problem (3.28). To this end, we set

[TABLE]

From Lemma 3.10 we easily know that $\mathcal{N}_{h},\mathcal{M}_{h}$ are bounded from above by a positive constant independent of the mesh size $h$.

Theorem 3.1**.**

The problem (3.28) admits at least one solution $(\bm{u}_{h},T_{h})\in\bm{W}_{h}\times S_{h}^{0}$.

Proof.

First, by Lemma 3.11 it is easy to see that, for a given $\bm{u}_{h}\in\bm{W}_{h}$, the bilinear form $\overline{A}_{h}(\bm{u}_{h};\cdot,\cdot)$ is continuous and coercive on $S_{h}^{0}\times S_{h}^{0}$. Hence, by Lax-Milgram theorem there is a unique $T_{h}\in S_{h}^{0}$ such that the second equation of (3.28) holds.

Define a mapping $F:\bm{W}_{h}\rightarrow S_{h}^{0}$ by $F(\bm{u}_{h})=T_{h}$. Then the thing left is to show that there exists at least one $\bm{u}_{h}\in\bm{W}_{h}$ such that

[TABLE]

Take $s_{h}=T_{h}$ in the second equation of (3.28), and apply (3.32) and (3.25) to get

[TABLE]

which yields

[TABLE]

Take $\bm{v}_{h}=\bm{u}_{h}$ in (3.33), and we obtain

[TABLE]

This indicates

[TABLE]

By Lemma 3.10 and (3.31), we also have

[TABLE]

Now we introduce another mapping, $\mathcal{A}:\bm{W}_{h}\rightarrow\bm{W}_{h}$, defined by $\mathcal{A}(\bm{u}_{h})=\bm{w}$, where $\bm{w}\in\bm{W}_{h}$ is determined by

[TABLE]

Clearly, $\bm{u}_{h}$ is a solution to (3.33) if it is a solution to

[TABLE]

To show this system has a solution, from the Leray-Schauder’s principle it suffices to prove the following two assertions: (i) $\mathcal{A}$ is a continuous and compact mapping; (ii) for any $0\leq\lambda\leq 1$, the set $\bm{W}_{\lambda,h}:=\{\bm{v}_{h}\in\bm{W}_{h}:\bm{v}_{h}=\lambda\mathcal{A}\bm{v}_{h}\}$ is bounded.

Let $\bm{u}_{1h},\bm{u}_{2h}\in\bm{W}_{h}$, set $\bm{w}_{2}=\mathcal{A}(\bm{u}_{2h})$ and $\bm{w}_{1}=\mathcal{A}(\bm{u}_{1h})$, then we obtain

[TABLE]

Subtracting (3.38) from (3.37), and taking $\bm{v}_{h}=\bm{w}:=\bm{w}_{1}-\bm{w}_{2}$, we get

[TABLE]

Substitute $T_{h}=F(\bm{u}_{1h})$ and $T_{h}=F(\bm{u}_{2h})$ into the second equation of (3.28), respectively, and subtract the two resultant equations each other, then, in view of (2.9), we have

[TABLE]

Taking $s_{h}=F(\bm{u}_{1h})-F(\bm{u}_{2h})$ in this equation, together with (2.4), (3.34), and Lemma 3.10, leads to

[TABLE]

As a result, from (3.39) and (3.35) it follows

[TABLE]

which means that $\mathcal{A}$ is equicontinuous and uniformly bounded. Thus, $\mathcal{A}$ is compact by the Arzelá-Ascoli theorem[4].

It remains to show (ii). If $\lambda=0$, then $\bm{W}_{\lambda,h}=\{0\}$. For $\lambda\in(0,1]$ and $\bm{v}_{h}\in\bm{W}_{\lambda,h}$, by (3.36) and (2.4) we have

[TABLE]

which implies

[TABLE]

This completes the proof. ∎

We now give a global uniqueness criteria for the case of small data (small Rayleigh number $Ra$).

Theorem 3.2**.**

Suppose

[TABLE]

Then the problem (3.28) admits a unique solution $(\bm{u}_{h},T_{h})\in\bm{W}_{h}\times S_{h}^{0}$ with $T_{h}=F(\bm{u}_{h})$.

Proof.

By Theorem 3.1, let $\bm{u}_{1h},\bm{u}_{2h}\in\bm{W}_{h}$ be two solutions to the problem (3.33). Then it suffices to show $\bm{u}_{1h}=\bm{u}_{2h}$. In fact, we have

[TABLE]

Subtracting the above two equations each other with $\bm{v}_{h}=\bm{u}_{1h}-\bm{u}_{2h}$, and using (2.4), we obtain

[TABLE]

which, together with Lemma 3.10, (3.40) and (3.35), yields

[TABLE]

If $\bm{u}_{1h}\neq\bm{u}_{2h}$, then, by the assumption (3.41), we further have

[TABLE]

which contradicts. Therefore $\bm{u}_{1h}=\bm{u}_{2h}$. ∎

4 A priori error estimates

This section is devoted to the error estimation of the WG scheme (2.7). We set

[TABLE]

We recall that $k\geq 1$ and $l=k,k-1$.

Lemma 4.1**.**

For any $\bm{w},\bm{u}\in\bm{W}$, $T\in H_{0}^{1}(\Omega)$, $\bm{v}_{h}\in\bm{V}_{h}^{0}$ and $s_{h}\in S_{h}^{0}$, it holds

[TABLE]

where

[TABLE]

Proof.

From the definition of weak divergence and Green’s formula, we have

[TABLE]

which, together with the definition of the trilinear form $c_{h}(\cdot;\cdot,\cdot)$, yields (4.1).

Similarly, we can obtain (4.2). ∎

Lemma 4.2**.**

Let $j,r$ be nonnegative integers. For any $K\in\mathcal{T}_{h}$ and $\bm{v}\in[H^{r}(K)]^{d}$, the following estimates hold for the RT projection operator:

[TABLE]

Proof.

We only prove (4.3), since the estimates (4.4)-(4.6) follow similarly.

For $1\leq r\leq j+1$, by the triangle inequality, the inverse inequality, Lemma 3.1, and Lemma 3.7, we get

[TABLE]

i.e. (4.3) holds. ∎

Lemma 4.3**.**

For $\bm{u}\in[H^{k+1}(\Omega_{f})]^{d}$ with $\nabla\cdot\bm{u}=0$ and $T\in H^{k+1}(\Omega)$, it holds

[TABLE]

for $l=k$ when $d=2,3$, and for $l=k-1$ when $d=2$.

Proof.

From the hölder inequality, the sobolev inequality, and the projection properties, we have

[TABLE]

For $l=k$ when $d=2,3$, and for $l=k-1$ when $d=2$, we have

[TABLE]

[TABLE]

and

[TABLE]

Similarly, we can obtain

[TABLE]

[TABLE]

For $l=k$ when $d=2,3$, and for $l=k-1$ when $d=2$, we have

[TABLE]

[TABLE]

As a result, the two desired results follow from the definitions of $E_{N}(\bm{u};\bm{u},\bm{v}_{h})$, $\overline{E}_{N}(\bm{u};T,s_{h})$ given in Lemma 4.1. ∎

Lemma 4.4**.**

Let $(\bm{u},p,T)$ be the solution to the problem (1.6), then it holds

[TABLE]

where

[TABLE]

In addition, it holds

[TABLE]

Proof.

We first show (4.9). In fact, for all $K\in\mathcal{T}_{h},\varphi\in P_{k}(K)$, by Lemma 3.8 we get

[TABLE]

which indicates

[TABLE]

Thus, the result (4.9) follows from Lemma 3.6.

By the definition of $a_{0}(\cdot,\cdot),A(\cdot,\cdot)$ and $d_{h}(\cdot,\cdot)$, we obtain

[TABLE]

From the commutativity property (3.12), the definition of weak gradient, Green’s formula, the property of the projection $\bm{Q}_{l}^{b}$, and the relation $\langle\nabla\bm{u}\cdot\bm{n},\bm{v}_{hb}\rangle_{\partial\mathcal{T}_{h}}=0$, it follows

[TABLE]

By the definitions of the projections $\bm{Q}_{l}^{b}$ and $Q_{k}^{0}$, we have

[TABLE]

[TABLE]

By (4.1), we get

[TABLE]

The commutativity property (3.13) gives

[TABLE]

In view of (4.10), (3.10), and the definitions of $d_{h}(\cdot,\cdot)$ and weak gradient, we obtain

[TABLE]

Finally, the desired relation (4.7) follows from the combination of (4.11)-(4.17) and the first equation of (1.6).

Similarly, we can get the relation (4.8). This completes the proof. ∎

Lemma 4.5**.**

For $\bm{u}\in[H^{k+1}(\Omega_{f})]^{d}$ and $T\in H^{k+1}(\Omega)$, it holds

[TABLE]

Proof.

By Lemma 3.1 and the definition of $\bm{Q}_{m}^{0}$, we have

[TABLE]

i.e. (4.18) holds. The estimate (4.19) follows similarly. ∎

Theorem 4.1**.**

Let $(\bm{u},p,T)\in[H^{k+1}(\Omega_{f})]^{d}\times H^{k}(\Omega_{f})\times H^{k+1}(\Omega)$ and $(\bm{u}_{h},p_{h},T_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}\times S_{h}^{0}$ be the solutions to the problem (1.6) and the WG scheme (2.7), respectively. Then, under the assumption (3.41) with

[TABLE]

it holds the following estimates: for $l=k$ when $d=2,3$, and for $l=k-1$ when $d=2$,

[TABLE]

where $C_{1}:=\frac{1+{\Pr}^{-1}\lVert\bm{u}\rVert_{2}}{C_{0}+\mathcal{M}_{h}Ra\kappa^{-2}\lVert g\rVert_{h}}$.

Proof.

From (2.7) and Lemma 4.4 we easily get the following error equations:

[TABLE]

[TABLE]

Take $(\bm{v}_{h},q_{h},s_{h})=(\bm{I}_{h}\bm{u}-\bm{u}_{h},J_{h}p-p_{h},H_{h}T-T_{h})$ in the above two equations, then we have

[TABLE]

[TABLE]

By (2.4) and (3.35), it holds

[TABLE]

which, together with (4.26), (4.20), Lemma 4.3, and Lemma 4.5, leads to

[TABLE]

i.e. (4.21) holds.

Similarly, we can obtain

[TABLE]

i.e. (4.22) holds.

Finally, let us estimate $\lVert J_{h}p-p_{h}\rVert$. In light of Theorem 3.12, (4.24), Lemma 3.10, Lemma 4.3, Lemma 4.5, (4.21), and (4.22), we have

[TABLE]

i.e. (4.23) holds. ∎

From Theorem 4.1, Lemma 3.3, and the triangle inequality, it follows the following error estimates:

Theorem 4.2**.**

Under the same conditions of Theorem 4.1, it holds

[TABLE]

5 Local elimination property and iteration scheme

5.1 Local elimination

In this subsection, we shall show that in the WG scheme (2.7), the velocity, pressure, and temperature approximations, $(\bm{u}_{h0},p_{h0},T_{h0})$, defined in the interior of the elements, can be locally eliminated by using the numerical traces, $(\bm{u}_{hb},p_{hb},T_{hb})$, defined on the interface of the elements. Therefore, after the local elimination the resultant system only involves degrees of freedom of $(\bm{u}_{hb},p_{hb},T_{hb})$ as unknowns.

We rewrite the scheme (2.7) as the following form: seek $\bm{u}_{h}=(\bm{u}_{h0},\bm{u}_{hb})\in\bm{V}_{h}^{0}$, $p_{h}=(p_{h0},p_{hb})\in Q_{h}^{0}$ and $T_{h}=(T_{h0},T_{hb})\in S_{h}^{0}$ such that

[TABLE]

For all $K\in\mathcal{T}_{h}^{f}$, taking $\bm{v}_{h0}|_{\mathcal{T}_{h}^{f}/K}=\bm{0},\bm{v}_{hb}=\bm{0}$, $q_{h0}|_{\mathcal{T}_{h}^{f}/K}=0,q_{hb}=0$ and $T_{h0}|_{\mathcal{T}_{h}^{f}/K}=0,T_{hb}=0$ in (5.3), we can get the following local problem: seek $(\bm{u}_{h0},p_{h0},T_{h0})\in[P_{k}(K)]^{d}\times P_{k-1}(K)\times P_{k}(K)$ such that, for $\forall(\bm{v}_{h0},q_{h0})\in[P_{k}(K)]^{d}\times P_{k-1}(K),s_{h0}\in P_{k}(K)$,

[TABLE]

where

[TABLE]

For any $K\in\mathcal{T}_{h}$, we define the following semi-norms:

[TABLE]

It is easy to see that the above semi-norms are norms on the local spaces $[P_{k}(K)]^{d}$ and $P_{k}(K)$, respectively.

By following the same routine as in Section 3 for the global problem (2.7), we can obtain the following existence and uniqueness results for the local problem (5.3).

Theorem 5.1**.**

For any given $\bm{u}_{hb},p_{hb}$ and $T_{hb}$, and any $K\in\mathcal{T}_{h}$, the local problem (5.6) admits at least one solution. In addition, it admits a unique solution if

[TABLE]

where

[TABLE]

and $\bm{W}_{h,K}:=\{\bm{w}_{h0}\in[P_{k}(K)]^{d}:b_{h,K}(\bm{w}_{h0},q_{h0})=0,\forall q_{h0}\in P_{k-1}(K)\}.$

5.2 Iteration scheme

Since the WG scheme (2.7) is nonlinear, we introduce the following Oseen’s iteration scheme: given $\bm{u}_{h}^{0}$, for $n=1,2,\cdots,$ and $\forall(\bm{v}_{h},q_{h},s_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}\times S_{h}^{0}$,

[TABLE]

We have the following convergence theorem.

Theorem 5.2**.**

Let $(\bm{u}_{h},p_{h},T_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}\times S_{h}^{0}$ be the solution to the WG scheme (2.7), and assume that (4.20) holds. Then the Oseen’s iteration scheme (5.9) is convergent in the following sense:

[TABLE]

Proof.

Set $\bm{e}_{u}^{n}:=\bm{u}_{h}^{n}-\bm{u}_{h},e_{p}^{n}:=p_{h}^{n}-p_{h},e_{T}^{n}:=T_{h}^{n}-T_{h}$, then, from (2.7) and (5.9), we have, for $\forall(\bm{v}_{h},q_{h},s_{h})\in\bm{V}_{h}^{0}\times Q_{h}^{0}\times S_{h}^{0}$,

[TABLE]

Taking $\bm{v}_{h}=\bm{e}_{u}^{n},q_{h}=e_{p}^{n},s_{h}=e_{T}^{n}$ in (5.12), in view of (2.4) and Lemma 3.10, we get

[TABLE]

[TABLE]

which, together with (3.34), (3.35), and (4.20), implies

[TABLE]

Since $0<C_{0}<1$, the above inequality leads to the conclusion

[TABLE]

Thus, from (5.14) and (3.34) it follows

[TABLE]

Finally, in light of Lemma 3.12 and the first equation of (5.12), we obtain

[TABLE]

which, together with (5.15) and (5.16), yields $\lim\limits_{n\longrightarrow\infty}\lVert p_{h}^{n}-p_{h}\rVert=0$. ∎

6 Numerical experiments

In this section, we shall show some numerical results to examine the performance of the proposed WG methods for the natural convection equations. The Oseen’s iteration scheme (5.9) with initial guess $\bm{u}_{h}^{0}=0$ is used in all the numerical experiments.

We consider three cases of our WG methods with $k=1,2$:

[TABLE]

Example 6.1**.**

Take $\Omega=[-1,1]\times[0,1]$ and $\Omega_{f}=[0,1]\times[0,1]$. The exact solution to the problem (1.6) is given by

[TABLE]

with $\Pr=1,\kappa=1,Ra=10$. Regular triangular meshes are used for the computation (see Figure 1).

Tables 1 and 2 show the history of convergence for the velocity $\bm{u}_{h0}$, pressure $p_{h0}$, and temperature $T_{h0}$. Results of $div_{h}U_{h}=:\max\limits_{K\in\mathcal{T}_{h}^{f}}h_{K}^{-1}\lVert\nabla\cdot\bm{u}_{h0}\rVert_{0,K}$ are also listed. From the numerical results we have the following observations:

• •

The convergence rates of $\lVert\nabla\bm{u}-\nabla_{h}\bm{u}_{h0}\rVert_{0},$ $\lVert p-p_{h0}\rVert_{0}$, and $\lVert\nabla T-\nabla_{h}T_{h0}\rVert_{0}$ for the proposed WG methods with $k=1,2$ are of $k^{\text{th}}$ orders, as is consistent with the theoretical results. In addition, the convergence rates of $\lVert\bm{u}-\bm{u}_{h0}\rVert_{0}$ and $\lVert T-T_{h0}\rVert_{0}$ are of $(k+1)^{\text{th}}$ orders.

• •

Since $\lVert\nabla_{h}\cdot\bm{u}_{h0}\rVert_{0,\infty}\apprle\max\limits_{K\in\mathcal{T}_{h}^{f}}h_{K}^{-1}\lVert\nabla\cdot\bm{u}_{h0}\rVert_{0,K}$, the velocity approximations obtained by our methods are globally divergence-free, which are conformable to the conclusion in Remark 2.1.

Example 6.2**.**

We consider the well-known test cave for the natural convection codes which is called buoyancy-driven cavity problem. This problem describes the two-dimensional flow of a Boussinesq fluid in an upright square cavity of side $L=1$. Fig.2 shows the physical domain with the boundary conditions. The velocity is zero on all the boundaries. The horizontal walls are insulated with $\frac{\partial T}{\partial\bm{n}}=0$, and the vertical sides are at temperatures $T_{H}=1$ and $T_{C}=0$. We take $\Omega=\Omega_{f}=[0,1]\times[0,1]$, $\kappa=1,\Pr=0.71,\bm{f}=0$, and $g=0$.

For different Rayleigh numbers, i.e. $Ra=10^{3},10^{4},10^{5},10^{6},10^{7}$, we use the WG-I method with $k=1,2$ to compute the following quantities at different mesh sizes:

The results are listed in Table 3 and compared with the famous benchmark solutions of de Vahl Davis [12] and of some other authors such as Manzari [21], Massarotti et al [22], Wan et al [29], and Zhang et al [38]. Figure 3 and Figure 4 show the contour maps of the stream function and the isotherms of the flow. We have the following observations:

• •

From Table 3 we can see that the WG-I method gives good results for all the quantities for different Rayleigh numbers. In particular, the method with $k=2$ behaves very well at the coarsest mesh $40\times 40$.

• •

Figure 3 demonstrates that, as Rayleigh number $Ra$ increases, the circular vortex at the cavity center begins to deform into an ellipse and then breaks up into two vortices, and then there’s a big vortex in the center.

• •

Figure 4 shows that, when Rayleigh number $Ra$ is small, the heat transfer mainly depends on heat conduction (isotherms almost vertical), with the increasing of $Ra$, the heat transfer pattern gradually turns to heat convection and boundary layers appear around the two walls (isotherms almost horizontal at the center).

7 Conclusions

In this paper, we have developed a class of weak Galerkin finite element methods with globally divergence-free velocity approximation for the steady-state natural convection problems. Well-posedness of the discrete scheme is analyzed, and optimal error estimates for the velocity, temperature and pressure approximations are derived. The proposed Oseen’s iteration algorithm is unconditionally convergent. Numerical experiments verify the theoretical results.

Acknowledgments

This work was supported by National Natural Science Foundation of China (11771312) and Major Research Plan of National Natural Science Foundation of China (91430105).