**2. Discrete Problem**

To pose a discrete problem, we can use a family {T*h*}*h*><sup>0</sup> of conforming triangulations to divide the domain Ω¯ such that Ω¯ = *<sup>T</sup>*∈T*<sup>h</sup> <sup>T</sup>*, ∀*h*, where *<sup>h</sup>* > 0 represents the mesh

**Citation:** González, M.; Varela, H. On the Adaptive Numerical Solution to the Darcy–Forchheimer Model. *Eng. Proc.* **2021**, *7*, 36. https:// doi.org/10.3390/engproc2021007036


Academic Editors: Joaquim de Moura, Marco A. González, Javier Pereira and Manuel G. Penedo

Published: 18 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

size. Here we follow [2] and choose the following conforming discrete subspaces of *X* and *M*, respectively:

$$X\_h := \left\{ \mathbf{v}\_h \in [L^2(\Omega)]^2 ; \forall T \in \mathcal{T}\_h, \mathbf{v}\_h|\_T \in [\mathbb{P}\_0(T)]^2 \right\} \subset X\_h$$

$$M\_h := Q\_h^1 \cap L\_0^2(\Omega) \subset M\_{\text{-}}$$

where *Q*<sup>1</sup> *<sup>h</sup>* := *qh* ∈ C0(Ω); <sup>∀</sup>*<sup>T</sup>* ∈ T*h*, *qh*|*<sup>T</sup>* <sup>∈</sup> <sup>P</sup>1(*T*) .

Then, the discrete problem consists in finding (**u***h*, *ph*) ∈ *Xh* × *Mh* such that

$$\begin{cases} \int\_{\Omega} \left( \frac{\mu}{\rho} K^{-1} \mathbf{u}\_{h} + \frac{\beta}{\rho} |\mathbf{u}\_{h}| \mathbf{u}\_{h} \right) \cdot \mathbf{v}\_{h} \, d\mathbf{x} + \int\_{\Omega} \nabla p\_{h} \cdot \mathbf{v}\_{h} \, d\mathbf{x} &= \int\_{\Omega} \mathbf{g} \cdot \mathbf{v}\_{h} \, d\mathbf{x}, \quad \forall \mathbf{v}\_{h} \in \mathcal{X}\_{h}, \\\ \int\_{\Omega} \nabla q\_{h} \cdot \mathbf{u}\_{h} \, d\mathbf{x} &= - \int\_{\Omega} q\_{h} f \, d\mathbf{x}, \quad \forall q\_{h} \in M\_{h}. \end{cases} \tag{2}$$

It is shown in [2] that problem (2) has a unique solution and that the sequence {(**u***h*, *ph*)}*<sup>h</sup>* converges to the exact solution of problem (1) in *X* × *M*. Furthermore, under additional regularity assumptions on the exact solution, some error estimates were derived in [2].

#### **3. Novel Error Estimator and Adaptive Algorithm**

We denote by E<sup>Ω</sup> , E*∂*<sup>Ω</sup> and E*T*, respectively, the sets of edges *e* belonging to the interior domain, the boundary and the element *T*; *he* denotes the length of a particular edge *e*; and *hT* is the diameter of a given element *T*. We denote by J*e*(*v*) the jump of *v* across the edge *e* in the direction of **n***e*, a fixed normal vector to side *e*. Finally, we use the operator <sup>A</sup>(**u***h*, *ph*) :<sup>=</sup> *<sup>μ</sup> <sup>ρ</sup> <sup>K</sup>*−1**u***<sup>h</sup>* <sup>+</sup> *<sup>β</sup> <sup>ρ</sup>* |**u***h*|**u***<sup>h</sup>* + ∇*ph* − **g**.

On every triangle *T* ∈ T*h*, we propose the following a posteriori error indicator:

$$\theta\_T = \left( h\_T^2 ||\tilde{\mathcal{A}}(\mathbf{u}\_h, p\_h)||\_{[L^2(T)]^2}^2 + ||\nabla \cdot \mathbf{u}\_h - f||\_{L^2(T)}^2 + \frac{1}{2} \sum\_{\varepsilon \in \mathcal{E}\_\Omega \cap \partial T} h\_T^{-1} ||\mathbb{J}\_\varepsilon(\mathbf{u}\_h \cdot \mathbf{n})||\_{L^2(\varepsilon)}^2 \\ \quad + \sum\_{\varepsilon \in \mathcal{E}\_{\partial T} \cap \partial T} h\_T^{-1} ||\mathbf{u}\_h \cdot \mathbf{n}||\_{L^2(\varepsilon)}^2 \right)^{1/2}$$

We also define the global a posteriori error indicator *θ* := ∑ *T*∈T*<sup>h</sup> θ*2 *T* 1/2 .

**Theorem 1.** *For the primal-mixed method* (2)*, there exists a positive constant C*1*, independent of h, and a positive constant C*2*, independent of h and T, such that*

$$||(\mathbf{u} - \mathbf{u}\_{h\prime}p - p\_h)||\_{X\times M} \le C\_1 \theta\_{\prime}$$

$$\theta\_T \le C\_2 ||(\mathbf{u} - \mathbf{u}\_{h\prime}p - p\_h)||\_{[L^3(w\tau)]^2 \times \mathbb{W}^{1,3/2}(w\tau)\prime} \quad \forall \ T \in \mathcal{T}\_h\text{-}$$

where *wT* = <sup>E</sup>*T*∩E*T* =<sup>Ø</sup> *T* .

We propose an adaptive algorithm based on the a posteriori error indicator *θ*. Given an initial mesh, we follow the iterative procedure described in Figure 1. Each new mesh is generated as suggested in [3].

**Figure 1.** Adaptive algorithm flux diagram.

#### **4. Numerical Experiment**

We performed several simulations in FreeFem++ [4], validating the theoretical results. Here we select an example on an L-shaped domain, <sup>Ω</sup> = (−1, 1)2\[0, 1] 2, and focus on the data *f* and **g** so that the exact solution is

$$p(\mathbf{x}, y) = \frac{1}{\mathbf{x} - \mathbf{1}.1}, \quad \mathbf{u}(\mathbf{x}, y) = \begin{pmatrix} \exp(\mathbf{x})\sin(y) \\ \exp(\mathbf{x})\cos(y) \end{pmatrix}. \tag{3}$$

Thus the solution has a singularity in pressure close to the line *x* = 1. Figure 2 shows the mesh refinement by the adaptive algorithm. Figure 3, bottom, represents the evolution with respect to degrees of freedom (DOF) of error and indicator; on the right, we can observe the evolution of the efficiency index with DOF.

**Figure 2.** Example 1. Initial mesh (270 DOF) on the (**top**); intermediate adapted mesh with 1512 DOF on the (**bottom**).

**Figure 3.** Example 1. (**Top**): Error and indicator evolution vs. DOF. (**Bottom**): Efficiency index vs. DOF.

#### **5. Discussion**

The adaptive algorithm was tested on an example with a singularity. From Figure 2 we can observe that the algorithm refined the mesh near the singularity, as expected. Since it is an academic example with a known solution, we could compute the exact error. The graphs in Figure 3 confirm that the error was lower for the adaptive refinement. Additionally, since the exact error and estimator followed close to parallel lines, we confirm that the indicator gives a consistent measure of the error. This could also be checked by the efficiency index, which is the ratio of indicator to exact total error.

**Funding:** The authors acknowledge the support of CITIC (FEDER Program and grant ED431G 2019/01). The research of M.G. is partially supported by Xunta de Galicia Grant GRC ED431C 2018-033. The research of H.V. is partially supported by Ministerio de Educación grant FPU18/06125.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

