Stabilized subgrid multiscale finite element formulation for advection-diffusion-reaction equation with variable coefficients coupled with Stokes-Darcy equation

TL;DR

This paper develops a stabilized finite element method for coupled advection-diffusion-reaction and Stokes-Darcy flow equations with variable coefficients, including error analysis and algebraic stabilization parameter approximation.

Contribution

It introduces a novel stabilized finite element formulation for coupled ADR and Stokes-Darcy equations with algebraic stabilization parameter approximation.

Findings

Effective stabilization of coupled equations demonstrated

Error estimates provided for the proposed method

Applicable to variable coefficient problems

Abstract

In this paper subgrid multiscale stabilized finite element method for Advection-Diffusion-Reaction (ADR) equation coupled with Stokes-Darcy flow problem has been studied. Here the advection velocity involved in ADR equation obeys Stokes-Darcy flow equation. In this study the approach of algebraic approximation of stabilization parameter has been considered. Further apriori error estimation has been elaborately carried out.

Tables

Table 1. Table 1: Comparison between Galerkin method and SGS method for small diffusion

	Galerkin Method	Galerkin Method	SGS Method	SGS Method
Mesh size	Error	Order of convergence	Error	Order of convergence
10	0.00049606		0.000265901
20	0.000152662	1.70017	6.36986 $e^{-5}$	2.06155
40	3.94719 $e^{-5}$	1.95144	1.68468 $e^{-5}$	1.91879
80	1.43178 $e^{-5}$	1.46302	4.73744 $e^{-6}$	1.83029
160	4.79682 $e^{-6}$	1.57766	5.4902 $e^{-7}$	1.91276
320	1.38697 $e^{-6}$	1.79015	4.45223 $e^{-8}$	1.89875

Table 2. Table 2: Comparison between Galerkin method and SGS method for diffusion dominated case

	Galerkin Method	Galerkin Method	SGS Method	SGS Method
Mesh size	Error	Order of convergence	Error	Order of convergence
10	0.000571444		0.000574645
20	0.000139961	2.02959	0.000141085	2.0261
40	3.20151 $e^{-5}$	2.1282	3.22881 $e^{-5}$	2.12749
80	7.74389 $e^{-6}$	2.04762	7.82478 $e^{-6}$	2.04488
160	2.02045 $e^{-6}$	1.93838	2.03639 $e^{-6}$	1.94204
320	4.69837 $e^{-7}$	2.10444	4.75272 $e^{-7}$	2.09919

Equations

$$\begin{array}{rl}-\mu \Delta \text{u}+\sigma \text{u}+▽p& =\text{f}in\Omega \\ ▽\cdot \text{u}& =0in\Omega \\ \text{u}& =\text{0}on\partial \Omega \end{array}$$

$$\begin{array}{rl}-▽\cdot \widetilde{▽}c+\text{u}\cdot ▽c+\alpha c& =gin\Omega \\ c& =0on\Omega \end{array}$$

$$\begin{array}{rl}{a}_S(\text{u},\text{v})-b(\text{v},p)+b(\text{u},q)& =l_S(\text{v})\forall \text{v}\in V\times V,q\in Q\\ a_T(c,d)& =l_T(d)\forall d\in V\end{array}$$

$$\begin{array}{rl}{a}_S({\text{u}}_h,{\text{v}}_h)-b({\text{v}}_h,p_h)+b({\text{u}}_h,q_h)& =l_S({\text{v}}_h)\forall {\text{v}}_h\in V_h\times V_h,q_h\in Q_h\\ a_T(c_h,d_h)& =l_T(d_h)\forall d_h\in V_h\end{array}$$

$$\begin{array}{rl}{a}_S({\text{u}}_h,{\text{v}}_h)-b({\text{v}}_h,p_h)+b({\text{u}}_h,q_h)+{\int }_{{\Omega }^{\prime }}(-{\mathcal{L}}_1^{\ast }{\text{V}}_h)[\tau ]{\mathcal{L}}_1{\text{U}}_h& =l_S({\text{v}}_h)+{\int }_{{\Omega }^{\prime }}(-{\mathcal{L}}_1^{\ast }{\text{V}}_h)[\tau ]\text{F}\\ \forall {\text{V}}_h& =({\text{v}}_h,q_h)\in V_h\times V_h\times Q_h\\ a_T(c_h,d_h)+{\int }_{{\Omega }^{\prime }}(-{\mathcal{L}}_2^{\ast }{d}_h){\tau }_3{\mathcal{L}}_2{c}_h& =l_T(d_h)+{\int }_{{\Omega }^{\prime }}(-{\mathcal{L}}_2^{\ast }{d}_h){\tau }_{3}g\\ & \forall d_h\in V_h\end{array}$$

$${\mathcal{L}}_1\text{U}=\left[\begin{array}{c}-\mu (c)\Delta \text{u}+\sigma \text{u}+▽p\\ ▽\cdot \text{u}\end{array}\right],itsadjoint{\mathcal{L}}_1^{\ast }\text{U}=\left[\begin{array}{c}-\mu (c)\Delta \text{u}+\sigma \text{u}-▽p\\ -▽\cdot \text{u}\end{array}\right],\text{F}=\left[\begin{array}{c}\text{f}\\ 0\end{array}\right]$$

$$[\tau ]=diag({\tau }_1,{\tau }_1,{\tau }_2)=\left[\begin{array}{ccc}(c_1\frac{\mu }{h^2}+\sigma {)}^{-1}& 0& 0\\ 0& (c_1\frac{\mu }{h^2}+\sigma {)}^{-1}& 0\\ 0& 0& c_2\mu \end{array}\right]and{\tau }_3=(\frac{9D}{4h^2}+\frac{3U}{2h}+\alpha {)}^{-1}$$

$$\Vert u-P_u^{h}u{\Vert }_l=\Vert E_u^I{\Vert }_l\le C(m,\Omega )h^{m+1-l}\Vert v{\Vert }_{m+1}$$

$$\begin{array}{rl}\Vert E_{u1}^A{\Vert }_V^2+\Vert E_{u2}^A{\Vert }_V^2+\Vert E_p^A{\Vert }_Q^2& \le C_1(\text{u},p)h^2\\ \Vert E_c^A{\Vert }_V^2& \le C_2(c)h^2\end{array}$$

$$\begin{array}{rl}{a}_S(\text{u},{\text{v}}_h)-b({\text{v}}_h,p)+b(\text{u},q_h)& =l_S({\text{v}}_h)\forall {\text{v}}_h\in V_h\times V_h,q_h\in Q_h\\ a_T(c,d_h)& =l_T(d_h)\forall d_h\in V_h\end{array}$$

$$a_S(\text{u}-{\text{u}}_h,{\text{v}}_h)-b({\text{v}}_h,p-p_h)+b(\text{u}-{\text{u}}_h,q_h)={\int }_{{\Omega }^{\prime }}({\mathcal{L}}_1^{\ast }{\text{V}}_h)[\tau ](\text{F}-{\mathcal{L}}_1{\text{U}}_h)$$

$$\begin{gathered} \mu {\int }_{\Omega }▽E_{\text{u}}^I:▽{\text{v}}_h+\mu {\int }_{\Omega }▽E_{\text{u}}^A:▽{\text{v}}_h+\sigma {\int }_{\Omega }{E}_{\text{u}}^I\cdot {\text{v}}_h+\sigma {\int }_{\Omega }{E}_{\text{u}}^A\cdot {\text{v}}_h-{\int }_{\Omega }(▽\cdot {\text{v}}_h)E_p^I- \\ {\int }_{\Omega }(▽\cdot {\text{v}}_h)E_p^A+{\int }_{\Omega }(▽\cdot E_{\text{u}}^I)q_h+{\int }_{\Omega }(▽\cdot E_{\text{u}}^A)q_h \\ ={\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta v_{1h}-\sigma v_{1h}+\frac{\partial q_h}{\partial x})(-\mu \Delta (E_{u1}^I+E_{u1}^A)+\sigma (E_{u1}^I+E_{u1}^A)+\frac{\partial E_p^I}{\partial x}+\frac{\partial E_p^A}{\partial x})+ \\ {\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta v_{2h}-\sigma v_{2h}+\frac{\partial q_h}{\partial y})(-\mu \Delta (E_{u2}^I+E_{u2}^A)+\sigma (E_{u2}^I+E_{u2}^A)+\frac{\partial E_p^I}{\partial y}+\frac{\partial E_p^A}{\partial y})+ \\ {\int }_{\Omega }{\tau }_2(▽\cdot {\text{v}}_h)(▽\cdot E_{\text{u}}^I)+{\int }_{\Omega }{\tau }_2(▽\cdot {\text{v}}_h)(▽\cdot E_{\text{u}}^A) \end{gathered}$$

$$\begin{gathered} \mu {\int }_{\Omega }▽E_{\text{u}}^A:▽E_{\text{u}}^A+\sigma {\int }_{\Omega }{E}_{\text{u}}^A\cdot E_{\text{u}}^A={\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta E_{u1}^A-\sigma E_{u1}^A+\frac{\partial E_p^A}{\partial x})(-\mu \Delta E_{u1}^I+\sigma E_{u1}^I+ \\ \frac{\partial E_p^I}{\partial x})+{\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta E_{u1}^A-\sigma E_{u1}^A+\frac{\partial E_p^A}{\partial x})(-\mu \Delta E_{u1}^A+\sigma E_{u1}^A+\frac{\partial E_p^A}{\partial x})+{\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta E_{u2}^A-\sigma E_{u2}^A \\ +\frac{\partial E_p^A}{\partial y})(-\mu \Delta E_{u2}^I+\sigma E_{u2}^I+\frac{\partial E_p^I}{\partial y})+{\int }_{{\Omega }^{\prime }}{\tau }_1(\mu \Delta E_{u2}^A-\sigma E_{u2}^A+\frac{\partial E_p^A}{\partial y})(-\mu \Delta E_{u2}^A+\sigma E_{u2}^A+\frac{\partial E_p^A}{\partial y}) \\ +{\int }_{{\Omega }^{\prime }}{\tau }_2(▽\cdot E_{\text{u}}^A)(▽\cdot E_{\text{u}}^I)+{\int }_{{\Omega }^{\prime }}{\tau }_2(▽\cdot E_{\text{u}}^A)(▽\cdot E_{\text{u}}^A)-\mu {\int }_{\Omega }▽E_{\text{u}}^I:▽E_{\text{u}}^A-{\int }_{\Omega }(▽\cdot E_{\text{u}}^I)E_p^A \\ =S_1+S_2+S_3+S_4+S_5+S_6+S_7+S_8(say) \end{gathered}$$

$$\begin{array}{rl}{S}_1& \le \mid {\tau }_1\mid \{(\sum _{k=1}^{n_{el}}({\mu }^2{M}_{2k}+\sigma \mid \mu \mid M_{1k}+\mid \mu \mid M_{3k}))\Vert E_{u1}^I{\Vert }_2+(\sum _{k=1}^{n_{el}}(\sigma \mid \mu \mid M_{2k}+\\ & {\sigma }^2{M}_{1k}+\mid \sigma \mid M_{3k}))\Vert E_{u1}^I\Vert +(\sum _{k=1}^{n_{el}}(\mid \mu \mid M_{2k}+\sigma M_{1k}+M_{3k}))\Vert E_p^I{\Vert }_1\}\\ & \le C\mid {\tau }_1\mid \{(\sum _{k=1}^{n_{el}}({\mu }^2{M}_{2k}+\sigma \mid \mu \mid M_{1k}+\mid \mu \mid M_{3k}))\Vert u_1{\Vert }_2+h^2(\sum _{k=1}^{n_{el}}(\sigma \mid \mu \mid M_{2k}\\ & +{\sigma }^2{M}_{1k}+\mid \sigma \mid M_{3k}))\Vert u_1{\Vert }_2+(\sum _{k=1}^{n_{el}}(\mid \mu \mid M_{2k}+\sigma M_{1k}+M_{3k}))\Vert p{\Vert }_1\}\end{array}$$

$$\begin{array}{rl}{S}_2& \le \mid {\tau }_1\mid \sum _{k=1}^{n_{el}}({\mu }^2{M}_{2k}^2+2\sigma \mid \mu \mid M_{1k}{M}_{2k}+{\sigma }^2{M}_{1k}^2+M_{3k}^2)\\ S_3& \le C\mid {\tau }_1\mid \{(\sum _{k=1}^{n_{el}}({\mu }^2{M}_{2k}^{\prime }+\sigma \mid \mu \mid M_{1k}^{\prime }+\mid \mu \mid M_{3k}^{\prime }))\Vert u_2{\Vert }_2+h^2(\sum _{k=1}^{n_{el}}(\sigma \mid \mu \mid M_{2k}^{\prime }\\ & +{\sigma }^2{M}_{1k}^{\prime }+\mid \sigma \mid M_{3k}^{\prime }))\Vert u_2{\Vert }_2+(\sum _{k=1}^{n_{el}}(\mid \mu \mid M_{2k}^{\prime }+\sigma M_{1k}^{\prime }+M_{3k}^{\prime }))\Vert p{\Vert }_1\}\\ S_4& \le \mid {\tau }_1\mid \sum _{k=1}^{n_{el}}({\mu }^2{M}_{2k}^{\prime 2}+2\sigma \mid \mu \mid M_{1k}^{\prime }{M}_{2k}^{\prime }+{\sigma }^2{M}_{1k}^{\prime 2}+M_{3k}^{\prime 2})\\ S_5& \le \mid {\tau }_2\mid Ch\sum _{k=1}^{n_{el}}(M_{4k}\Vert u_1{\Vert }_2+M_{5k}^{\prime }\Vert u_2{\Vert }_2)\\ S_6& \le \mid {\tau }_2\mid \sum _{k=1}^{n_{el}}(M_{4k}^2+M_{5k}^{\prime 2})\end{array}$$

$$\begin{array}{rl}{S}_7& =-\mu {\int }_{\Omega }(\frac{\partial E_{u1}^I}{\partial x}\frac{\partial E_{u1}^A}{\partial x}+\frac{\partial E_{u1}^I}{\partial y}\frac{\partial E_{u1}^A}{\partial y}+\frac{\partial E_{u2}^I}{\partial x}\frac{\partial E_{u2}^A}{\partial x}+\frac{\partial E_{u2}^I}{\partial y}\frac{\partial E_{u2}^A}{\partial y})\end{array}$$

$$\begin{array}{rl}& \le \mid \mu \mid \{\frac{1}{2{ϵ}_1}(\Vert \frac{\partial E_{u1}^I}{\partial x}{\Vert }^2+\Vert \frac{\partial E_{u1}^I}{\partial y}{\Vert }^2+\Vert \frac{\partial E_{u2}^I}{\partial x}{\Vert }^2+\Vert \frac{\partial E_{u2}^I}{\partial y}{\Vert }^2)+\frac{{ϵ}_1}{2}(\Vert \frac{\partial E_{u1}^A}{\partial x}{\Vert }^2+\\ & \Vert \frac{\partial E_{u1}^A}{\partial y}{\Vert }^2+\Vert \frac{\partial E_{u2}^A}{\partial x}{\Vert }^2+\Vert \frac{\partial E_{u2}^A}{\partial y}{\Vert }^2)\}\\ & \le \mid \mu \mid \{C\frac{h^2}{2{ϵ}_1}(\Vert u_1{\Vert }_2^2+\Vert u_2{\Vert }_2^2)+\frac{{ϵ}_1}{2}(\Vert \frac{\partial E_{u1}^A}{\partial x}{\Vert }^2+\Vert \frac{\partial E_{u1}^A}{\partial y}{\Vert }^2+\Vert \frac{\partial E_{u2}^A}{\partial x}{\Vert }^2+\Vert \frac{\partial E_{u2}^A}{\partial y}{\Vert }^2)\}\end{array}$$

$$\begin{array}{rl}{S}_8={\int }_{\Omega }(\frac{\partial E_{u1}^I}{\partial x}+\frac{\partial E_{u2}^I}{\partial y})E_p^A& \le (\Vert \frac{\partial E_{u1}^I}{\partial x}\Vert +\Vert \frac{\partial E_{u2}^I}{\partial y}\Vert )\Vert E_p^A\Vert \\ & \le C\frac{h^2}{{ϵ}_2}(\Vert u_1{\Vert }_2^2+\Vert u_2{\Vert }_2^2)+\frac{{ϵ}_2}{2}\Vert E_p^A{\Vert }^2\end{array}$$

$$\Vert E_{u1}^A{\Vert }_V^2+\Vert E_{u2}^A{\Vert }_V^2+\Vert E_p^A{\Vert }_Q^2\le C(\text{u},p)h^2$$

$$a_T(c-c_h,d_h)={\int }_{{\Omega }^{\prime }}({\mathcal{L}}_2^{\ast }{d}_h){\tau }_3(g-{\mathcal{L}}_2{c}_h)={\int }_{{\Omega }^{\prime }}({\mathcal{L}}_2^{\ast }{d}_h){\tau }_3{\mathcal{L}}_2(c-c_h)$$

$$\begin{gathered} {\int }_{\Omega }\widetilde{▽}{E}_c^I\cdot ▽d_h+{\int }_{\Omega }\widetilde{▽}{E}_c^A\cdot ▽d_h+{\int }_{\Omega }{d}_h\text{u}\cdot ▽E_c^I+{\int }_{\Omega }{d}_h\text{u}\cdot ▽E_c^A+\alpha {\int }_{\Omega }{E}_c^I{d}_h+ \\ \alpha {\int }_{\Omega }{E}_c^A{d}_h={\int }_{{\Omega }^{\prime }}(-▽\cdot \widetilde{▽}{d}_h-\overline{u}\cdot ▽d_h+\alpha d_h){\tau }_3(-▽\cdot \widetilde{▽}{E}_c^I+\overline{u}\cdot ▽E_c^I+\alpha E_c^I) \\ {\int }_{{\Omega }^{\prime }}(-▽\cdot \widetilde{▽}{d}_h-\overline{u}\cdot ▽d_h+\alpha d_h){\tau }_3(-▽\cdot \widetilde{▽}{E}_c^A+\overline{u}\cdot ▽E_c^A+\alpha E_c^A) \end{gathered}$$

$$\begin{gathered} {\int }_{\Omega }\widetilde{▽}{E}_c^A\cdot ▽E_c^A+\alpha {\int }_{\Omega }{E}_c^A{E}_c^A={\int }_{{\Omega }^{\prime }}(-▽\cdot \widetilde{▽}{E}_c^A-\overline{u}\cdot ▽E_c^A+\alpha E_c^A){\tau }_3(-▽\cdot \widetilde{▽}{E}_c^I+ \\ \overline{u}\cdot ▽E_c^I+\alpha E_c^I)+{\int }_{{\Omega }^{\prime }}(-▽\cdot \widetilde{▽}{E}_c^A-\overline{u}\cdot ▽E_c^A+\alpha E_c^A){\tau }_3(-▽\cdot \widetilde{▽}{E}_c^A+\overline{u}\cdot ▽E_c^A+\alpha E_c^A) \\ -{\int }_{\Omega }\widetilde{▽}{E}_c^I\cdot ▽E_c^A-{\int }_{\Omega }{E}_c^A\text{u}\cdot ▽E_c^I-{\int }_{\Omega }{E}_c^A\text{u}\cdot ▽E_c^A \\ =T_1+T_2+T_3+T_4+T_5(say) \end{gathered}$$

$$\begin{array}{rl}{T}_1& \le \mid {\tau }_3\mid (\sum _{k=1}^{n_{el}}(D_{1l}{M}_{8k}+D_{2l}{M}_{8k}^{\prime }+D_{u1}{M}_{7k}+D_{u2}{M}_{7k}^{\prime }+\alpha M_{6k}))\{(D_{1l}+D_{2l})\\ & \Vert E_c^I{\Vert }_2+({\bar{D}}_{u1}+{\bar{D}}_{u2})\Vert E_c^I{\Vert }_1+\alpha \Vert E_c^I\Vert \}\\ & \le \mid {\tau }_3\mid C(\sum _{k=1}^{n_{el}}(D_{1l}{M}_{8k}+D_{2l}{M}_{8k}^{\prime }+D_{u1}{M}_{7k}+D_{u2}{M}_{7k}^{\prime }+\alpha M_{6k}))\{(D_{1l}+D_{2l})\\ & +h({\bar{D}}_{u1}+{\bar{D}}_{u2})+h^2\alpha \}\Vert c{\Vert }_2\end{array}$$

$$\begin{array}{rl}{T}_2& \le \mid {\tau }_3\mid \sum _{k=1}^{n_{el}}\{(D_{1l}{M}_{8k}+D_{2l}{M}_{8k}^{\prime }+D_{u1}{M}_{7k}+D_{u2}{M}_{7k}^{\prime }+\alpha M_{6k})(D_{1l}{M}_{8k}+\\ & D_{2l}{M}_{8k}^{\prime }+{\bar{D}}_{u1}{M}_{7k}+{\bar{D}}_{u2}{M}_{7k}^{\prime }+\alpha M_{6k})\}\end{array}$$

$$\begin{array}{rl}{T}_3& =-{\int }_{\Omega }(D_1\frac{\partial E_c^I}{\partial x}\frac{\partial E_c^A}{\partial x}+D_2\frac{\partial E_c^I}{\partial y}\frac{\partial E_c^A}{\partial y})\\ & \le \frac{D_{1l}}{2{ϵ}_3}\Vert E_c^I{\Vert }_1^2+\frac{D_{2l}{ϵ}_3}{2}(\Vert \frac{\partial E_c^I}{\partial x}{\Vert }^2+\Vert \frac{\partial E_c^I}{\partial y}{\Vert }^2)\\ & \le \frac{D_{1l}}{2{ϵ}_3}{h}^2\Vert c{\Vert }_2^2+\frac{D_{2l}{ϵ}_3}{2}(\Vert \frac{\partial E_c^I}{\partial x}{\Vert }^2+\Vert \frac{\partial E_c^I}{\partial y}{\Vert }^2)\end{array}$$

$$\begin{array}{rl}{T}_4& =-{\int }_{\Omega }(u_1\frac{\partial E_c^I}{\partial x}+u_2\frac{\partial E_c^I}{\partial y})E_c^A\\ & \le \frac{(u_{1l}+u_{2l})}{2{ϵ}_4}\Vert E_c^I{\Vert }_1^2+\frac{(u_{1l}+u_{2l})}{2}{ϵ}_4\Vert E_c^A{\Vert }^2\\ & \le \frac{(u_{1l}+u_{2l})}{2{ϵ}_4}{h}^2\Vert c{\Vert }_2^2+\frac{(u_{1l}+u_{2l})}{2}{ϵ}_4\Vert E_c^A{\Vert }^2\end{array}$$

$$\begin{array}{rl}{T}_5& =-{\int }_{\Omega }(u_1\frac{\partial E_c^A}{\partial x}+u_2\frac{\partial E_c^A}{\partial y})E_c^A\\ & \le \frac{(u_{1l}+u_{2l})}{2{ϵ}_4}(\Vert \frac{\partial E_c^A}{\partial x}{\Vert }^2+\Vert \frac{\partial E_c^A}{\partial y}{\Vert }^2)+\frac{(u_{1l}+u_{2l})}{2}{ϵ}_4\Vert E_c^A{\Vert }^2\end{array}$$

$$\Vert \frac{\partial E_c^A}{\partial x}{\Vert }^2+\Vert \frac{\partial E_c^A}{\partial y}{\Vert }^2+\Vert E_c^A{\Vert }^2\le C_2(c)h^2$$

$$\begin{array}{rl}\Vert u_1-u_{1h}{\Vert }_V^2+\Vert u_2-u_{2h}{\Vert }_V^2+\Vert p-p_h{\Vert }_Q^2+\Vert c-c_h{\Vert }_V^2& \le C^{\prime }(\text{u},p,c)h^2\end{array}$$

$$\begin{gathered} \Vert u_1-u_{1h}{\Vert }_V^2+\Vert u_2-u_{2h}{\Vert }_V^2+\Vert p-p_h{\Vert }_Q^2+\Vert c-c_h{\Vert }_V^2 \\ =\Vert E_{u1}^I+E_{u1}^A{\Vert }_V^2+\Vert E_{u2}^I+E_{u2}^A{\Vert }_V^2+\Vert E_c^I+E_c^A{\Vert }_V^2+\Vert E_c^I+E_c^A{\Vert }_V^2 \\ \le \bar{C}(\Vert E_{u1}^I{\Vert }_V^2+\Vert E_{u2}^I{\Vert }_V^2+\Vert E_p^I{\Vert }_V^2+\Vert E_c^I{\Vert }_V^2+\Vert E_{u1}^A{\Vert }_V^2+\Vert E_{u2}^A{\Vert }_V^2+\Vert E_p^A{\Vert }_V^2+\Vert E_c^A{\Vert }_V^2) \\ \le C^{\prime }(\text{u},p,c)h^2 \end{gathered}$$

Full text

Keywords Stokes-Darcy-Brinkman equation $\cdot$ Advection-diffusion-reaction equation $\cdot$ Subgrid scale method $\cdot$ A priori error estimation

1 Introduction

Transport phenomena, mathematically expressed in terms of Advection-Diffusion-Reaction(ADR) equation has always been an active area of research in the fields of Biomedical Engineering, Environmental Sciences, Chemical Engineering etc. Our previous works [4],[6] have focused on transport equation with spatially variable diffusion and advection coefficients. In this paper we have considered only diffusion coefficients as variable and the advection velocity comes from Stokes-Darcy flow problem. There are few studies [1], [2] available in literature on coupled Stokes-Darcy transport equation, but none of them has studied stabilized subgrid multiscale finite element method for both the equations. Being most general stabilization method now-a-days subgrid method has been considered more suitable finite element method to be dealt with. It involves split of the unknown true solution is chosen to be the standard Galerkin finite element solution and obtaining the unresolvable solution in terms of the known resolvable one, we will finally arrive at the subgrid formulation. [4] has derived an expression of the stabilization parameter for subgrid formulation of ADR equation with variable coefficients and another study [5] has found out the same for Stokes-Darcy equation. Both these studies have followed the approach of algebraic approximation of the parameter. Combining those results we will have the required stabilization parameter for this coupled problem. Here further we have derived apriori error estimation the coupled problem.

The paper is organised as: section 2 presents model problem along with weak formulation and subgrid formulation for the weakly coupled problem. In the next section we have carried out apriori error estimation for this stabilized method with the help of interpolation estimates for introducing projection operators.

2 Statement of the problem

Let $\Omega$ be an open bounded domain in $R^{d}$, d=2,3. Here we have worked with two dimensional model for the sake of simplicity in further derivation, but it can be extended to three dimensional model straightforward.

For an incompressible solvent fluid the Stokes-Darcy (or Brinkman) equation representing its flow is to find $u:\Omega\rightarrow R^{2}$ and $p:\Omega\rightarrow R$ such that

[TABLE]

where u= ($u_{1},u_{2}$),p and $\mu$ are velocity, pressure and viscosity of the fluid, $\sigma$ is the inverse of permeability and f is the body force.

And the ADR equation with spatially variable coefficients representing the concentration c of transporting solute in the same domain $\Omega$ along with homogeneous Dirichlet boundary condition is to find c: $\Omega\rightarrow R$ such that,

[TABLE]

where the notation, $\tilde{\bigtriangledown}:=(D_{1}\frac{\partial}{\partial x},D_{2}\frac{\partial}{\partial y})$

$D_{1},D_{2}$ are spatially variable diffusion coefficients along x-axis and y-axis respectively, $\alpha$ is the reaction coefficient and g denotes the source of solute mass. Here we make some assumptions on the coefficients as follows:

(i) $\mu$, $\sigma$ and $\alpha$ are positive constants.

(ii) The diffusion coefficients $D_{1},D_{2}$ are continuous functions on bounded domain $\Omega$and hence bounded. Let $D_{1l},D_{2l}$ be upper bounds of them respectively.

(iii) The body force $\textbf{f}\in(L^{2}(\Omega))^{2}$ and source $g\in L^{2}(\Omega)$

2.1 Weak formulation

Let us consider the standard space $V=H^{1}_{0}(\Omega)$ as the admissible space for both velocity fields and concentration and the space $Q=L^{2}(\Omega)$ for pressure.

The weak formulation is to find $(\textbf{u},p)\in V\times V\times Q$ and $c\in V$ such that

[TABLE]

where $a_{S}(\textbf{u},\textbf{v})=\int_{\Omega}\mu\bigtriangledown\textbf{u}:\bigtriangledown\textbf{v}+\sigma\int_{\Omega}\textbf{u}\cdot\textbf{v}$ and   $b(\textbf{v},q)=\int_{\Omega}(\bigtriangledown\cdot\textbf{v})q$

$a_{T}(c,d)=\int_{\Omega}\tilde{\bigtriangledown}c\cdot\bigtriangledown d+\int_{\Omega}d\textbf{u}\cdot\bigtriangledown c+\alpha\int_{\Omega}cd$

$l_{S}(\textbf{v})=\int_{\Omega}\textbf{f}\cdot\textbf{v}$ and $l_{T}(d)=\int_{\Omega}gd$

The stability of the continuous problem (3) has been shown in [9] and [6].

2.2 Subgrid formulation

Let $\Omega$ be discretized into total $n_{el}$ number of sub-domains $\Omega_{k}$ for k=1,2,…,$n_{el}$ and $h_{k}$ be the diameter of each sub-domain respectively. Let h= $\underset{k=1,2,...,n_{el}}{max}h_{k}$ and $\tilde{\Omega}=\bigcup_{k=1}^{n_{e}l}\Omega_{k}$ be the union of interior elements. Let $V_{h}$ and $Q_{h}$ be the suitable finite dimensional subspaces of $V$ and $Q$ respectively where

$V_{h}=\{v\in V:v(\Omega_{k})=\mathcal{P}^{2}(\Omega_{k})\}$ and $Q_{h}=\{q\in Q_{s}:q(\Omega_{k})=\mathcal{P}^{1}(\Omega_{k})\}$

where $\mathcal{P}^{1}(\Omega_{k})$ and $\mathcal{P}^{2}(\Omega_{k})$ denote complete polynomial of order 1 and 2 respectively over each $\Omega_{k}$ for k=1,2,…,$n_{el}$.

Now the standard Galerkin finite element formulation of (3) is to find $(\textbf{u}_{h},p_{h})\in V_{h}\times V_{h}\times Q_{h}$ and $c\in V_{h}$ such that

[TABLE]

As we earlier mention that for subgrid formulation the finite element solution will be taken as resolvable scale and we will obtain the subgrid formulation by expressing the unresolvable scale in terms of known solution. Hence the subgrid formulation for (3) is to find $(\textbf{u}_{h},p_{h})\in V_{h}\times V_{h}\times Q_{h}$ and $c\in V_{h}$ such that

[TABLE]

where

[TABLE]

$\mathcal{L}_{2}c=-\bigtriangledown\cdot\widetilde{\bigtriangledown}c+\overline{u}\cdot\bigtriangledown c+\alpha c$ and its $adjoint$ $\mathcal{L}_{2}^{*}c=-\bigtriangledown\cdot\widetilde{\bigtriangledown}c-\overline{u}\cdot\bigtriangledown c+\alpha c$

the stabilization parameters[],[] are obtained as

[TABLE]

3 Error estimation

In this section we are going to derive apriori estimate in $V$ norm that is in standard $H^{1}$ norm. Before that we introduce projection operators and carry out error splitting as follows:

Let $\textbf{e}=(e_{\textbf{u}},e_{p},e_{c})$ denote the error where the components are $e_{\textbf{u}}=(e_{u1},e_{u2})=(u_{1}-u_{1h},u_{2}-u_{2h}),e_{p}=(p-p_{h})$ and $e_{c}=(c-c_{h})$. Let us introduce the projection operator corresponding to each component as the following,

(i)For any $\textbf{u}\in V\times V$ let there exist an interpolation $P^{h}_{\textbf{u}}:V\times V\longrightarrow V_{h}\times V_{h}$ satisfying $b(\textbf{u}-P^{h}_{\textbf{u}}\textbf{u},q_{h})=0$    $\forall q_{h}\in Q_{h}$ and component wise satisfy $L^{2}$-orthogonality condition i.e. for any $u_{i}\in V$   $(u_{i}-P^{h}_{u_{i}}u_{i},v_{h})=0$   $\forall v_{h}\in V_{s}^{h}$ for i=1,2.

(ii) Let $P^{h}_{p}:Q\longrightarrow Q_{h}$ be the $L^{2}$ orthogonal projection given by

$\int_{\Omega}(p-P^{h}_{p}p)q_{h}=0$   $\forall p\in Q,\hskip 2.84526pt\forall q_{h}\in Q_{h}$

(iii) Similarly let $P^{h}_{c}:V\longrightarrow V_{h}$ satisfy $\int_{\Omega}(c-P^{h}_{c}c)v_{h}=0$   $\forall c\in V,\hskip 2.84526pt\forall v_{h}\in V_{h}$

Now each components of the error can be split into two parts interpolation part, $E^{I}$ and auxiliary part, $E^{A}$ as follows:

$e_{u1}=(u_{1}-u_{1h})=(u_{1}-P^{h}_{u1}u_{1})+(P^{h}_{u1}u_{1}-u_{1h})=E^{I}_{u1}+E^{A}_{u1}$

Similarly, $e_{u2}=E^{I}_{u2}+E^{A}_{u2}$, $e_{p}=E^{I}_{p}+E^{A}_{p}$, and $e_{c}=E^{I}_{c}+E^{A}_{c}$

Interpolation estimates [7, 8] : for any true solution $u$ with regularity upto (m+1)

[TABLE]

where l is a positive integer and C is a constant depending on m and the domain. For l=0 and 1 it implies standard $L^{2}(\Omega)$ and $H^{1}(\Omega)$ norms respectively. We will use $\|\cdot\|$ instead of $\|\cdot\|_{0}$ to denote $L^{2}(\Omega)$ norm.

Now we carry out apriori error estimate in two parts. In first part we will bound auxiliary error and in the later part using the result obtained in first part we will find final form of apriori estimate.

Theorem 1**.**

Auxiliary error estimate: For velocity $\textbf{u}_{h}=(u_{1h},u_{2h})$, pressure $p_{h}$ and concentration $c_{h}$ belonging to $V_{h}\times V_{h},Q_{h},V_{h}$ respectively satisfying (5), assume sufficient regularity of exact solution in equations (1)-(2). Then there exists constants $C_{1}$ and $C_{2}$, depending upon u,p and c respectively, such that

[TABLE]

Proof.

We divide the proof into two parts: first part contains estimation of auxiliary error for velocity and pressure components separately and in the later part we will find bound for auxiliary error corresponding to concentration. Now for true solution (5) becomes

[TABLE]

First part: Subtracting first equation of (8) from that of (5) we will have, $\forall\textbf{v}_{h}\in V_{h}\times V_{h},q_{h}\in Q_{h}$

[TABLE]

After explicitly writing the bilinear forms and using the error splitting we have the above equation as follows

[TABLE]

Since this holds for all $\textbf{v}_{h}\in V_{h}\times V_{h}$ and $q_{h}\in Q_{h}$, replacing them by $E^{A}_{\textbf{u}}$ and $E^{A}_{p}$ respectively and using properties of projection operators, we finally have

[TABLE]

Now we estimate each of the above terms separately. Before that here we like to make an important observation: by the virtue of chosen finite element spaces it can be clearly said that every element belonging to $V_{h}$ and $Q_{h}$ and their derivatives all are bounded functions over each sub-domain $\Omega_{k}$. Therefore let us consider the positive numbers $M_{1k},M_{2k},M_{3k},M_{4k},M_{5k},M_{1k}^{\prime},M_{2k}^{\prime},M_{3k}^{\prime},M_{4k}^{\prime},M_{5k}^{\prime}$ as bounds for $E^{A}_{u1},\Delta E^{A}_{u1},\frac{\partial E^{A}_{p}}{\partial x},\frac{\partial E^{A}_{u1}}{\partial x},\frac{\partial E^{A}_{u1}}{\partial y},E^{A}_{u2},\Delta E^{A}_{u2},\frac{\partial E^{A}_{p}}{\partial y},\frac{\partial E^{A}_{u2}}{\partial x},\frac{\partial E^{A}_{u2}}{\partial y}$ respectively.

To bound $T_{1}$ simply multiplying the terms, then applying Cauchy-Schwarz inequality over each of them and finally using bounds for $E^{A}s$ we will have

[TABLE]

Applying similar arguments we have the following bounds

[TABLE]

Applying Cauchy-Schwarz and Young’s inequality we get

[TABLE]

[TABLE]

Applying similar argument

[TABLE]

This completes finding bounds for each of the terms in the right hand side of (11). Now we put all the results obtained above into the equation (11) and take out those terms, which are common with the terms of $LHS$, in the $LHS$ from the $RHS$ of the equation. Consequently in the $RHS$ we remain with the terms multiplied by $h^{2}$ as the stabilization parameters are of order $h^{2}$. Now we consider $\epsilon_{1},\epsilon_{2}$ in such a manner that all the coefficients in $LHS$ can be made positive. After that taking minimum of all those coefficients and dividing both the sides by that we finally arrive at the completion of the proof as follows:

[TABLE]

Second part: Subtracting second equation of (8) from that of (5) we will have, $\forall d_{h}\in V_{h}$

[TABLE]

After explicitly writing the bilinear forms and using the error splitting we have the above equation as follows

[TABLE]

Since this holds for all $d_{h}\in V_{h}$, replacing it by $E^{A}_{c}$ in the above equation and using properties of projection operators, we finally have

[TABLE]

We will follow the same procedure as the previous part. Under the same observation let us consider the positive numbers $M_{6k},M_{7k},M_{8k},M_{7k}^{\prime},M_{8k}^{\prime}$ as bounds for $E^{A}_{c},\frac{\partial E^{A}_{c}}{\partial x},\frac{\partial^{2}E^{A}_{c}}{\partial x^{2}},\frac{\partial E^{A}_{c}}{\partial x},\frac{\partial E^{A}_{c}}{\partial y}$ respectively over each sub-domain.

First simply multiplying the terms of $T_{1}$, then applying Cauchy-Schwarz inequality over each of them and finally using bounds for $E^{A}s$ we can find estimate of $T_{1}$ as follows:

[TABLE]

Applying similar arguments we have the following bounds

[TABLE]

For estimating next terms we are going to use Cauchy-Schwarz and Young’s inequalities as follows:

[TABLE]

[TABLE]

Similarly

[TABLE]

Now combining all the results obtained above into the equation (20) and proceeding as explained at the end of the first part we finally get

[TABLE]

This completes the proof. ∎

Theorem 2**.**

Apriori error estimate: Assuming the same condition as in the previous theorem

[TABLE]

where $C^{\prime}$ is a constant depending upon $\textbf{u},p,c$.

Proof.

We prove this result by first applying triangle inequality and then interpolation estimates and the results obtained in the previous theorem as the following

[TABLE]

This completes the proof. ∎

Remark 1**.**

This error estimation is not computable as it depends upon the true solution, but it implies convergence of the stabilized method and the order of convergence for this method is 2.

4 Numerical Experiment

In this section through a test case based on hydrological importance [10], we have shown subgrid multiscale stabilization method is performing better than Galerkin method for small diffusion for this weakly coupled Stokes-Darcy/ transport equation. We have considered a source of pollutant dispersed into an incompressible fluid in a simple bounded square domain $\Omega=(0,1)\times(0,1)$. The advection velocity comes from the solution of Stokes-Darcy flow problem and the diffusion coefficients are taken as $D_{1}=10^{-7}(1+0.02x)^{2}$, $D_{2}=10^{-8}(1+0.02y)^{2}$ and the reaction term $\alpha=10$. Table 1 represents the comparison between Galerkin method and SGS method for small diffusion.

Again for higher diffusion coefficients compared to the reaction term both the Galerkin method and SGS method perform well. Let us consider the values of diffusion coefficients as $D_{1}=(1+0.02x)^{2}$, $D_{2}=0.1(1+0.02y)^{2}$ and reaction term $\alpha=0.001$. Table 2 presents the comparison of both Galerkin method and SGS method for diffusion dominated case.

Remark 2**.**

The numerical experiment has proved theoretically established result and shown that the order of convergence under SGS method for both the cases is 2 whereas for small diffusion Galerkin method oscillates and for diffusion dominated flow it behaves well as stabilized method.

5 Conclusion

In this paper we have studied subgrid multiscale stabilization formulation of weakly coupled Stokes-Darcy and transport equation. As well as we have carried out apriori error estimation, which implies the convergence of the method and finally the numerical experiment verifies the theoretically established result.