∎

11institutetext: Jin Zhang 22institutetext: School of Mathematics and
Statistics, Shandong Normal University, Jinan 250014, China
22email<EMAIL_ADDRESS>33institutetext: Corresponding author:
Xiaowei Liu 44institutetext: School of Mathematics and Statistics, Qilu
University of Technology (Shandong Academy of Sciences), Jinan 250353, China
44email<EMAIL_ADDRESS>
# Convergence of a finite element method on a Bakhvalov-type mesh for a
singularly perturbed convection–diffusion equation in 2D††thanks: This
research is partially supported by National Natural Science Foundation of
China (11771257,11601251).

Jin Zhang  Xiaowei Liu

(Received: date / Accepted: date)

###### Abstract

A finite element method of any order is applied on a Bakhvalov-type mesh to
solve a singularly perturbed convection–diffusion equation in 2D, whose
solution exhibits exponential boundary layers. A uniform convergence of
(almost) optimal order is proved by means of a carefully defined interpolant.

###### Keywords:

Singular perturbationConvection–diffusion equation Bakhvalov-type mesh Finite
element method Uniform convergence.

###### MSC:

65N12 65N30 65N50

##  1 Introduction

Consider the elliptic boundary value problem

$\begin{split}-\varepsilon\Delta u-\boldsymbol{b}\cdot\nabla
u+cu=&f\quad\text{ in $\Omega=(0,1)^{2}$},\\\ u=&0\quad\text{on
$\partial\Omega$},\end{split}$ (1)

where $\varepsilon\ll 1$ is a small positive parameter and
$\boldsymbol{b}(x,y)=(b_{1}(x,y),b_{2}(x,y))^{T}$. The functions
$b_{1},b_{2},c$ and $f$ are assumed to be smooth on $\bar{\Omega}$. We also
assume that

$b_{1}(x,y)\geq\beta_{1}>0,\,b_{2}(x,y)\geq\beta_{2}>0,\,c(x,y)+\frac{1}{2}\nabla\cdot\boldsymbol{b}(x,y)\geq\gamma>0\quad\text{on
$\bar{\Omega},$}$ (2)

where $\beta_{1}$, $\beta_{2}$ and $\gamma$ are some constants. These
conditions ensure that (1) has a unique solution in $H^{1}_{0}(\Omega)\cap
H^{2}(\Omega)$ for all $f\in L^{2}(\Omega)$ (see, e.g.,
Roo1Sty2Tob3:2008-Robust ). Because $\varepsilon$ is small, the problem is in
general singularly perturbed and its solution typically has exponential
boundary layers at $x=0$ and $y=0$ and a corner layer at $(0,0)$.

Layer phenomena appears in different kinds of problems, for example singularly
perturbed problems, which is one of the important topics in scientific
computing. If a priori knowledge of layers has been obtained from asymptotic
analysis etc., different kinds of meshes could be designed for uniform
convergent numerical methods (see Roo1Sty2Tob3:2008-Robust ;
Mil1Rio2Shi3:2012-Fitted ). Here “uniform” means that the convergence is
independent of the singular perturbation parameter. Bakhvalov-type meshes are
one of the most popular layer-adapted meshes and usually have better numerical
performances than Shishkin-type meshes—another popular layer-adapted meshes (
see Linb:2010-Layer ).

However, it is far from mature for convergence theories of finite element
methods on Bakhvalov-type meshes. One of the main reasons is that the standard
Lagrange interpolant does not work for Bakhvalov-type meshes (see
Roos-2006-Error ). In Zhan1Liu2:2020-Optimal , we gave a new idea for
convergence analysis on Bakhvalov-type meshes in 1D. Here we extend the
analysis to two dimensions. The extension is not trivial, because we must pay
attention to the construction of the interpolant used in the case of 2D. This
interpolant is carefully defined according to the characteristics of layer
functions and the structures of Bakhvalov-type meshes. Besides, different from
1D case, we must pay attention to the homogenous Dirichlet boundary condition
when the idea in Zhan1Liu2:2020-Optimal  is applied to 2D case. The
interpolation errors are derived in a delicate way. Then almost uniform
convergence of optimal order is proved for finite element methods.

The rest of the paper is organized as follows. In Section 2 we describe the
assumptions on the regularity of the solution, introduce a Bakhvalov-type mesh
and define a finite element method of any order. Some preliminary results for
the subsequent analysis are also given in this section. In Section 3 we
construct and analyze an interpolant to the solution for uniform convergence
on the Bakhvalov-type mesh. In Section 4 almost uniform convergence of optimal
order is obtained by means of the interpolant and careful analysis of the
convective term in the bilinear form. In Section 5, numerical results
illustrate our theoretical bounds.

Denote by $\|\cdot\|_{\infty,D}$ the norms in the Lebesgue space
$L^{\infty}(D)$. In $L^{2}(D)$, the inner product and the $L^{2}(D)$-norm are
denoted by $(\cdot,\cdot)_{D}$ and $\|\cdot\|_{D}$, respectively. In
$H^{1}(D)$, the seminorms are denoted by $|\cdot|_{1}$. Here $D$ is any
measurable subset of $\Omega$. When $D=\Omega$ we drop the subscript $D$ from
the notation for simplicity. Throughout the article, all constants $C$ are
independent of $\varepsilon$ and the mesh parameter $N$ and may take different
values in different formulas.

##  2 Decomposition of the solution, Bakhvalov-type mesh and finite element
method

In this section we present a decomposition of the solution to (1), introduce a
Bakhvalov-type mesh and define a finite element method. Some preliminary
inequalities are also presented. In the subsequent analysis, let $k$ be a
fixed integer with $k\geq 1$.

###  2.1 Regularity of the solution

We make the following assumption about the solution $u$ to (1), which
describes the layer structure of $u$. This assumption is also used in
Fra1Lin2Roo3etc:2010-Uniform .

###### Assumption 1

The solution $u$ of (1) can be decomposed as

$u=S+E_{1}+E_{2}+E_{12},$ (3a) where $S$ is the smooth part of $u$, $E_{1}$
and $E_{2}$ are exponential layers along the sides $x=0$ and $y=0$ of $\Omega$
respectively, while $E_{12}$ is an exponential corner layer at $(0,0)$.
Moreover,there exists a constant $C$ such that for all $(x,y)\in\bar{\Omega}$
and $0\leq i+j\leq k+1$ one has
$\displaystyle\left|\frac{\partial^{i+j}S}{\partial x^{i}\partial
y^{j}}(x,y)\right|$ $\displaystyle\leq C,$ (3b)
$\displaystyle\left|\frac{\partial^{i+j}E_{1}}{\partial x^{i}\partial
y^{j}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-i}e^{-\beta_{1}x/\varepsilon},$ (3c)
$\displaystyle\left|\frac{\partial^{i+j}E_{2}}{\partial x^{i}\partial
y^{j}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-j}e^{-\beta_{2}y/\varepsilon},$
$\displaystyle\left|\frac{\partial^{i+j}E_{12}}{\partial x^{i}\partial
y^{j}}(x,y)\right|$ $\displaystyle\leq
C\varepsilon^{-(i+j)}e^{-(\beta_{1}x+\beta_{2}y)/\varepsilon}.$ (3d)

###### Remark 1

For the case $k=1$, the existence of this decomposition of $u$ with the bounds
on derivatives can be guaranteed by conditions on the data of the problem (1)
(see Linb1Styn2:2001-Asymptotic ). The arguments in Linb1Styn2:2001-Asymptotic
make this assumption with $k\geq 2$ credible if we impose sufficient
compatibility conditions on $f$ (see some explanations in (Stynes:2005-Steady,
, Sect.7)).

###  2.2 Bakhvalov-type mesh

Bakhvalov mesh first appeared in Bakhvalov:1969-Towards  and is graded in the
layer. Its applications require the solution of a nonlinear equation. To avoid
this difficulty, Bakhvalov-type meshes are proposed as approximations of
Bakhvalov mesh (see Linb:2010-Layer ).

Let $N$ be an even positive integer. We introduce a Bakhvalov-type mesh in the
$x$-direction

$0=x_{0}<x_{1}<\ldots<x_{N-1}<x_{N}=1.$

To resolve the layer along $x=0$, the mesh is graded in $[x_{0},x_{N/2}]$ and
equidistant in $[x_{N/2},1]$. The mesh points $x_{i}$ is defined by

$x_{i}=\left\\{\begin{aligned}
&\frac{\sigma\varepsilon}{\beta_{1}}\varphi(t_{i})\text{ with
$t_{i}=i/N$}\quad&&\text{for $i=0,1,\ldots,N/2$},\\\
&1-(1-x_{N/2})2(N-i)/N\quad&&\text{for
$i=N/2+1,\ldots,N$},\end{aligned}\right.$ (4)

with $\sigma\geq k+1$ and $\varphi(t):=-\ln(1-2(1-\varepsilon)t)$. The
parameter $\sigma$ determines the smallness of the layer terms in $x_{N/2}$.
This Bakhvalov-type mesh is also analyzed in Roos-2006-Error . In a similar
way we define the mesh $\\{y_{j}\\}_{j=0}^{N}$ along the $y$-direction except
that we replace $\beta_{1}$ by $\beta_{2}$ in (4). Then we obtain a tensor-
product rectangular mesh $\mathcal{T}_{N}$ with mesh points $(x_{i},y_{j})$
(see Figure 1).

###### Assumption 2

Assume that $\varepsilon\leq N^{-1}$ in our analysis, as is not a restriction
in practice.

Figure 1: Bakhvalov-type mesh $\mathcal{T}_{N}$.

Moreover, assume $N\geq 4$ is an even integer. Set $h_{i,x}:=x_{i+1}-x_{i}$
and $h_{j,y}:=y_{j+1}-y_{j}$ for all $i,j$. A mesh rectangle is often written
as $\tau_{i,j}=[x_{i},x_{i+1}]\times[y_{j},y_{j+1}]$ for a specific element
and more simply as $\tau$ for a generic mesh rectangle.

According to (Zhan1Liu2:2020-Optimal, , Lemma 3), we have the following lemma.

###### Lemma 1

For the Bakhvalov-type mesh (4), one has

$\displaystyle h_{0,x}\leq h_{1,x}\leq\ldots\leq h_{N/2-2,x},$ (5)
$\displaystyle C\varepsilon N^{-1}\leq h_{0,x}\leq C\varepsilon N^{-1},$ (6)
$\displaystyle\frac{1}{4}\sigma\varepsilon\leq
h_{N/2-2,x}\leq\sigma\varepsilon,$ (7)
$\displaystyle\frac{1}{2}\sigma\varepsilon\leq h_{N/2-1,x}\leq 2\sigma
N^{-1},$ (8) $\displaystyle N^{-1}\leq h_{i,x}\leq 2N^{-1}\quad N/2\leq i\leq
N-1,$ (9) $\displaystyle x_{N/2-1}\geq C\sigma\varepsilon\ln N,\quad
x_{N/2}\geq C\sigma\varepsilon|\ln\varepsilon|,$ (10) $\displaystyle
h_{i,x}^{\mu}e^{-\beta_{1}x_{i}/\varepsilon}\leq
C\varepsilon^{\mu}N^{-\mu}\quad\text{for $0\leq i\leq N/2-2$ and
$0\leq\mu\leq\sigma$}.$ (11)

For $h_{j,y}$, $0\leq j<N$, bounds analogous to (5)–(11) also hold.

###  2.3 Finite element method

On the above Bakhvalov-type mesh, define the finite element space by

$V^{N}:=\\{v^{N}\in C(\bar{\Omega}):\;\;v^{N}|_{\partial\Omega}=0\text{ and
$v^{N}|_{\tau}\in\mathcal{Q}_{k}(\tau)$ }\;\forall\tau\in\mathcal{T}_{N}\\},$

where $\mathcal{Q}_{k}(\tau)=\mathrm{span}\\{x^{i}y^{j}:\;0\leq i,j\leq k\\}$.

The finite element method is defined as follows: Find $u^{N}\in V^{N}$ such
that

$a(u^{N},v^{N})=(f,v^{N})\quad\forall v^{N}\in V^{N},$ (12)

with $a(u^{N},v^{N}):=\varepsilon(\nabla u^{N},\nabla
v^{N})+(-\boldsymbol{b}\cdot\nabla u^{N}+cu^{N},v^{N}).$ Condition (2) implies
the coercivity

$a(v^{N},v^{N})\geq\alpha\|v^{N}\|_{\varepsilon}^{2}\quad\text{for all
$v^{N}\in V^{N}$},$ (13)

where $\alpha=\min\\{1,\gamma\\}$ and

$\|v\|_{\varepsilon}:=\left\\{\varepsilon|v|^{2}_{1}+\|v\|^{2}\right\\}^{1/2}\quad\forall
v\in H^{1}(\Omega).$

It follows that there exists a unique solution $u^{N}$ for problem (12) from
Lax-Milgram lemma. Clearly, (1) and (12) imply the Galerkin orthogonality
property

$a(u-u^{N},v^{N})=0\quad\text{for all $v^{N}\in V^{N}$}.$ (14)

##  3 Interpolation and its errors

A new interpolation operator is introduced for uniform convergence. Set
$x_{i}^{s}:=x_{i}+(s/k)h_{i,x}$ and $y_{j}^{t}:=y_{j}+(t/k)h_{j,y}$ for
$i,j=0,1,\ldots,N-1$ and $s,t=0,\ldots,k-1$. For the consistency of notation,
set $x_{N}^{0}=x_{N}$ and $y_{N}^{0}=y_{N}$. For any $v\in
C^{0}(\bar{\Omega})$ its standard Lagrange interpolant $v^{I}\in V^{N}$ on the
Bakhvalov-type mesh can be written in the following form

$\displaystyle v^{I}(x,y)=$
$\displaystyle\sum_{i=0}^{N-1}\sum_{s=0}^{k-1}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}v(x_{i}^{s},y_{j}^{t})\theta_{i,j}^{s,t}(x,y)+v(x_{i}^{s},y_{N}^{0})\theta_{i,N}^{s,0}(x,y)\right)$
$\displaystyle+\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}v(x_{N}^{0},y_{j}^{t})\theta_{N,j}^{0,t}(x,y)+v(x_{N}^{0},y_{N}^{0})\theta_{N,N}^{0,0}(x,y),$

where $\theta_{i,j}^{s,t}(x,y)\in V^{N}$ is the piecewise $k$th–order hat
function associated with the point $(x_{i}^{s},y_{j}^{t})$. For the solution
$u$ to (1), recall (3a) in Assumption 1 and define the interpolant $\Pi u$ by

$\Pi u=S^{I}+\pi_{1}E_{1}+\pi_{2}E_{2}+\pi_{12}E_{12}.$ (15)

Here $S^{I}$ is the Lagrange interpolant to $S$ and

$\displaystyle(\pi_{i}E_{i})(x,y)=$ $\displaystyle
E^{I}_{i}-\mathcal{P}_{i}E_{i}-\mathcal{B}_{i}E_{i}\quad\text{for $i=1,2$},$
(16) $\displaystyle(\pi_{12}E_{12})(x,y)=$ $\displaystyle
E^{I}_{12}-\mathcal{P}_{12}E_{12}$

where

$\displaystyle\mathcal{P}_{1}E_{1}=\sum_{i=N/2-1}\sum_{s=0}^{k-1}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}E_{1}(x_{i}^{s},y_{j}^{t})\theta_{i,j}^{s,t}(x,y)+E_{1}(x_{i}^{s},y_{N}^{0})\theta_{i,N}^{s,0}(x,y)\right),$
(17)
$\displaystyle\mathcal{P}_{2}E_{2}=\sum_{j=N/2-1}\sum_{t=0}^{k-1}\left(\sum_{i=0}^{N-1}\sum_{s=0}^{k-1}E_{2}(x_{i}^{s},y_{j}^{t})\theta_{i,j}^{s,t}(x,y)+E_{2}(x_{N}^{0},y_{j}^{t})\theta_{N,j}^{0,t}(x,y)\right),$
(18)
$\displaystyle\mathcal{P}_{12}E_{12}=\sum_{i=N/2-1}\sum_{j=N/2-1}\left(\sum_{s=0}^{k-1}\sum_{t=0}^{k-1}E_{12}(x_{i}^{s},y_{j}^{t})\theta_{i,j}^{s,t}(x,y)\right)$
(19)

and

$\displaystyle\mathcal{B}_{1}E_{1}=\sum_{i=N/2-1}\sum_{s=0}^{k-1}\sum_{j=0,N}E_{1}(x_{i}^{s},y_{j}^{0})\theta_{i,j}^{s,0}(x,y),$
(20)
$\displaystyle\mathcal{B}_{2}E_{2}=\sum_{i=0,N}\sum_{j=N/2-1}\sum_{t=0}^{k-1}E_{2}(x_{i}^{0},y_{j}^{t})\theta_{i,j}^{0,t}(x,y).$

Clearly we have

$\displaystyle\Pi u\in V^{N},\quad\Pi
u=u^{I}-\sum_{i=1,2,12}\mathcal{P}_{i}E_{i}.$ (21)

###### Remark 2

The definitions of interpolation operators $\pi_{i}$, $i=1,2,12$ arise from
layer functions $E_{i}$, $i=1,2,12$ and the construction of Bakhvalov-type
mesh (4). The idea for operators $\mathcal{P}_{1}$, $\mathcal{P}_{2}$ and
$\mathcal{P}_{12}$ is same to the operator $\mathcal{P}$ in
Zhan1Liu2:2020-Optimal , that is, new interpolations for layer functions are
zero at certain degrees of freedom. The operators $\mathcal{B}_{1}$ and
$\mathcal{B}_{2}$ are introduced in order to maintain homogenous Dirichlet
boundary conditions.

From (Apel:1999-Anisotropic, , Theorem 2.7), we have the following anisotropic
interpolation results.

###### Lemma 2

Let $\tau\in\mathcal{T}_{N}$ and $v\in H^{k+1}(\tau)$.Then there exists a
constant $C$ such that Lagrange interpolation $v^{I}$ satisfies

$\displaystyle\|v-v^{I}\|_{\tau}\leq$ $\displaystyle
C\sum_{i+j=k+1}h^{i}_{x,\tau}h^{j}_{y,\tau}\left\|\frac{\partial^{k+1}v}{\partial
x^{i}\partial y^{j}}\right\|_{\tau}$
$\displaystyle\|(v-v^{I})_{x}\|_{\tau}\leq$ $\displaystyle
C\sum_{i+j=k}h^{i}_{x,\tau}h^{j}_{y,\tau}\left\|\frac{\partial^{k+1}v}{\partial
x^{i+1}\partial y^{j}}\right\|_{\tau}$

and

$\|(v-v^{I})_{y}\|_{\tau}\leq
C\sum_{i+j=k}h^{i}_{x,\tau}h^{j}_{y,\tau}\left\|\frac{\partial^{k+1}v}{\partial
x^{i}\partial y^{j+1}}\right\|_{\tau},$

where $h_{x,\tau}$ and $h_{y,\tau}$ denote the lengths along $x-$axis and
$y-$axis of the rectangle $\tau$, respectively.

Set

$\displaystyle\Omega_{11}:=[x_{0},x_{N/2-1}]\times[y_{0},y_{N/2-1}],\quad\Omega_{12}:=[x_{N/2-1},x_{N}]\times[y_{0},y_{N/2-1}],$
$\displaystyle\Omega_{21}:=[x_{0},x_{N/2-1}]\times[y_{N/2-1},y_{N}],\quad\Omega_{22}:=[x_{N/2-1},x_{N}]\times[y_{N/2-1},y_{N}].$

###### Lemma 3

Let Assumptions 1 and 2 hold true. Let $E^{I}_{i}$, $i=1,2$, denote the
Lagrange interpolants of $E_{i}$, $i=1,2$, respectively, on the Bakhvalov-type
mesh $\mathcal{T}_{N}$. Then there exists a constant $C$ such that the
following interpolation error estimates hold true:

$\displaystyle\|E_{1}-E^{I}_{1}\|_{\Omega\setminus\big{(}(x_{N/2-1},x_{N/2})\times[0,1]\big{)}}+\|E_{2}-E^{I}_{2}\|_{\Omega\setminus\big{(}[0,1]\times(y_{N/2-1},y_{N/2})\big{)}}\leq
C\varepsilon^{1/2}N^{-(k+1/2)},$
$\displaystyle\|E_{1}-E^{I}_{1}\|+\|E_{2}-E^{I}_{2}\|\leq CN^{-(k+1)},$
$\displaystyle\|E_{1}-E^{I}_{1}\|_{\varepsilon}+\|E_{2}-E^{I}_{2}\|_{\varepsilon}\leq
CN^{-k},$
$\displaystyle\|\mathcal{P}_{1}E_{1}\|_{\varepsilon}+\|\mathcal{P}_{2}E_{2}\|_{\varepsilon}+\|\mathcal{B}_{1}E_{1}\|_{\varepsilon}+\|\mathcal{B}_{2}E_{2}\|_{\varepsilon}\leq
CN^{1/2-\sigma},$

where $\mathcal{P}_{i}E_{i}$, $i=1,2$, are defined in (17) and (18),
respectively.

###### Proof

We just consider $E_{1}$, since $E_{2}$ can be analyzed in a similar way. To
consider $\|E_{1}-E^{I}_{1}\|$, we decompose it as follows

$\displaystyle\|E_{1}-E^{I}_{1}\|^{2}=$
$\displaystyle\|E_{1}-E^{I}_{1}\|^{2}_{[x_{0},x_{N/2-1}]\times[0,1]}+\|E_{1}-E^{I}_{1}\|^{2}_{[x_{N/2-1},x_{N/2}]\times[0,1]}$
$\displaystyle+\|E_{1}-E^{I}_{1}\|^{2}_{[x_{N/2},x_{N}]\times[0,1]}$
$\displaystyle=:$ $\displaystyle A_{1}+A_{2}+A_{3}.$

Note $h_{i,x}\leq C\varepsilon$ for $i=0,\ldots,N/2-2$. Lemmas 1 and 2 yield

$\displaystyle A_{1}=$
$\displaystyle\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\|E_{1}-E^{I}_{1}\|^{2}_{\tau_{i,j}}\leq
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k+1}h^{2l}_{i,x}h^{2m}_{j,y}\left\|\frac{\partial^{k+1}E_{1}}{\partial
x^{l}\partial y^{m}}\right\|^{2}_{\tau_{i,j}}$ (22) $\displaystyle\leq$
$\displaystyle
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k+1}h^{2l}_{i,x}h^{2m}_{j,y}\;(\varepsilon^{-2l}e^{-2\beta_{1}x_{i}/\varepsilon}h_{i,x}h_{j,y})$
$\displaystyle\leq$ $\displaystyle
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k+1}(\varepsilon^{2l}N^{-2l})\;h_{j,y}^{2m+1}\varepsilon^{1-2l}$
$\displaystyle\leq$ $\displaystyle C\varepsilon N^{-(2k+1)}.$

Now we consider the term $A_{2}$. Set $D_{0}:=[x_{N/2-1},x_{N/2}]\times[0,1]$.
Recall $|E_{1}(x_{N/2-1}^{i},y_{j}^{t})|\leq CN^{-\sigma}$. Then we have

$\displaystyle\|E_{1}^{I}\|_{D_{0}}^{2}\leq$ $\displaystyle
CN^{-2\sigma}\left(\sum_{i=N/2-1}\sum_{s=0}^{k-1}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}\|\theta_{i,j}^{s,t}\|_{D_{0}}^{2}+\|\theta_{i,N}^{s,0}\|_{D_{0}}^{2}\right)\right)$
(23)
$\displaystyle+CN^{-2\sigma}\left(\sum_{i=N/2}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}\|\theta_{i,j}^{0,t}\|_{D_{0}}^{2}+\|\theta_{i,N}^{0,0}\|_{D_{0}}^{2}\right)\right)$
$\displaystyle\leq$ $\displaystyle
CN^{-2\sigma}\sum_{j=0}^{N-1}h_{N/2-1,x}h_{j,y}\leq CN^{-(2\sigma+1)},$

where we have used Lemma 1. Direct calculations yield

$\|E_{1}\|_{D_{0}}\leq C\varepsilon^{1/2}N^{-\sigma}.$ (24)

From (23) and (24) we obtain

$A_{2}\leq C(\varepsilon+N^{-1})N^{-2\sigma}.$ (25)

Note $|E_{1}(x_{i}^{s},y_{j}^{t})|\leq C\varepsilon^{\sigma}$ for $i\geq N/2$.
Then the triangle inequality and Hölder inequalities yield

$\displaystyle A_{3}\leq$ $\displaystyle
C\left(\|E_{1}\|^{2}_{[x_{N/2},x_{N}]\times[0,1]}+\|E_{1}^{I}\|^{2}_{[x_{N/2},x_{N}]\times[0,1]}\right)$
(26) $\displaystyle\leq$ $\displaystyle
C\left(\|E_{1}\|^{2}_{\infty,[x_{N/2},x_{N}]\times[0,1]}+\|E_{1}^{I}\|^{2}_{\infty,[x_{N/2},x_{N}]\times[0,1]}\right)$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{2\sigma}.$

Collecting (22), (25) and (26), we prove the first and the second bounds.

Now we consider
$|E_{1}-E_{1}^{I}|^{2}_{1}=\|(E_{1}-E_{1}^{I})_{x}\|^{2}+\|(E_{1}-E_{1}^{I})_{y}\|^{2}$.
Lemmas 1 and 2 give

$\displaystyle\|(E_{1}-E_{1}^{I})_{x}\|^{2}_{\Omega_{11}\cup\Omega_{21}}=\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\|(E_{1}-E^{I}_{1})_{x}\|^{2}_{\tau_{i,j}}$
(27) $\displaystyle\leq$ $\displaystyle
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k}h^{2l}_{i,x}h^{2m}_{j,y}\left\|\frac{\partial^{k+1}E_{1}}{\partial
x^{l+1}\partial y^{m}}\right\|^{2}_{\tau_{i,j}}$ $\displaystyle\leq$
$\displaystyle
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k}h^{2l}_{i,x}h^{2m}_{j,y}\;(\varepsilon^{-2(l+1)}e^{-2\beta_{1}x_{i}/\varepsilon}h_{i,x}h_{j,y})$
$\displaystyle\leq$ $\displaystyle
C\sum_{i=0}^{N/2-2}\sum_{j=0}^{N-1}\sum_{l+m=k}(\varepsilon^{2l+1}N^{-(2l+1)})\;h_{j,y}^{2m+1}\varepsilon^{-2(l+1)}$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{-1}N^{-2k}.$

Note $\|(E_{1})_{x}\|_{\Omega_{12}\cup\Omega_{22}}\leq
C\varepsilon^{-1/2}N^{-\sigma}$. Then from the triangle inequality one has

$\displaystyle\|(E_{1}-E_{1}^{I})_{x}\|^{2}_{\Omega_{12}\cup\Omega_{22}}\leq
2\|(E_{1})_{x}\|^{2}_{\Omega_{12}\cup\Omega_{22}}+2\|(E_{1}^{I})_{x}\|^{2}_{\Omega_{12}\cup\Omega_{22}}$
(28) $\displaystyle\leq$ $\displaystyle
2\|(E_{1})_{x}\|^{2}_{\Omega_{12}\cup\Omega_{22}}+2\sum_{j=0}^{N-1}\|(E_{1}^{I})_{x}\|^{2}_{\tau_{N/2-1,j}}+2\sum_{i=N/2}^{N-1}\sum_{j=0}^{N-1}\|(E_{1}^{I})_{x}\|^{2}_{\tau_{i,j}}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}N^{-2\sigma}+C\varepsilon^{2\sigma}N^{2},$

where inverse inequalities (Ciarlet:1978-finite, , Theorem 3.2.6) and Lemma 1
yield

$\displaystyle\sum_{j=0}^{N-1}\|(E_{1}^{I})_{x}\|^{2}_{\tau_{N/2-1,j}}+\sum_{i=N/2}^{N-1}\sum_{j=0}^{N-1}\|(E_{1}^{I})_{x}\|^{2}_{\tau_{i,j}}$
$\displaystyle\leq$ $\displaystyle
C\sum_{j=0}^{N-1}h^{-2}_{N/2-1,x}\|E_{1}^{I}\|^{2}_{\tau_{N/2-1,j}}+C\sum_{i=N/2}^{N-1}\sum_{j=0}^{N-1}h^{-2}_{i,x}\|E_{1}^{I}\|^{2}_{\tau_{i,j}}$
$\displaystyle\leq$ $\displaystyle
C\sum_{j=0}^{N-1}h^{-2}_{N/2-1,x}\|E_{1}^{I}\|^{2}_{\infty,\tau_{N/2-1,j}}h_{N/2-1,x}h_{j,y}+C\sum_{i=N/2}^{N-1}\sum_{j=0}^{N-1}h^{-2}_{i,x}\|E_{1}^{I}\|^{2}_{\infty,\tau_{i,j}}h_{i,x}h_{j,y}$
$\displaystyle\leq$ $\displaystyle
C\varepsilon^{-1}N^{-2\sigma}+C\varepsilon^{2\sigma}N^{2}.$

Similar to the derivations of (27), we have

$\|(E_{1}-E_{1}^{I})_{y}\|^{2}_{\Omega_{11}\cup\Omega_{21}}\leq C\varepsilon
N^{-2k}.$ (29)

Note $\|(E_{1})_{y}\|_{\Omega_{12}\cup\Omega_{22}}\leq
C\varepsilon^{1/2}N^{-\sigma}$ and $h_{j,y}\geq C\varepsilon N^{-1}$ for
$j=0,\ldots,N-1$. Then one has

$\displaystyle\|(E_{1}-E_{1}^{I})_{y}\|^{2}_{\Omega_{12}\cup\Omega_{22}}\leq$
$\displaystyle
2\|(E_{1})_{y}\|^{2}_{\Omega_{12}\cup\Omega_{22}}+2\|(E_{1}^{I})_{y}\|^{2}_{\Omega_{12}\cup\Omega_{22}}$
(30) $\displaystyle\leq$ $\displaystyle C\varepsilon
N^{-2\sigma}+C\varepsilon^{-1}N^{2-2\sigma},$

where

$\displaystyle\|(E_{1}^{I})_{y}\|^{2}_{\Omega_{12}\cup\Omega_{22}}=$
$\displaystyle\sum_{i=N/2-1}^{N-1}\sum_{j=0}^{N-1}\|(E_{1}^{I})_{y}\|^{2}_{\tau_{i,j}}\leq
C\sum_{i=N/2-1}^{N-1}\sum_{j=0}^{N-1}h^{-2}_{j,y}\|E_{1}^{I}\|^{2}_{\tau_{i,j}}$
$\displaystyle\leq$ $\displaystyle
C\sum_{i=N/2}^{N-1}\sum_{j=0}^{N-1}h^{-2}_{j,y}(N^{-2\sigma}h_{i,x}h_{j,y})$
$\displaystyle\leq$ $\displaystyle
CN^{1-2\sigma}\max_{j}h_{j,y}^{-1}\leq\varepsilon^{-1}N^{2-2\sigma}.$

Collecting (27)–(30) and considering $\|E_{1}-E_{1}^{I}\|\leq CN^{-(k+1)}$, we
prove the third bound.

Now we consider $\|\mathcal{P}_{1}E_{1}\|_{\varepsilon}$. From (17) and
$|E(x_{N/2-1}^{i},y_{j}^{t})|\leq CN^{-\sigma}$, we can easily obtain

$\displaystyle\|\mathcal{P}_{1}E_{1}\|_{\varepsilon}^{2}\leq$ $\displaystyle
CN^{-2\sigma}\sum_{s=0}^{k-1}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}\|\theta_{N/2-1,j}^{s,t}\|^{2}_{\varepsilon}+\|\theta_{N/2-1,N}^{s,0}\|^{2}_{\varepsilon}\right)$
(31) $\displaystyle\leq$ $\displaystyle
CN^{-2\sigma}\sum_{j=0}^{N-1}\left(\varepsilon
h^{-1}_{N/2-1,x}h_{j,y}+\varepsilon
h_{N/2-1,x}h^{-1}_{j,y}+h_{N/2-1,x}h_{j,y}\right)$ $\displaystyle\leq$
$\displaystyle CN^{1-2\sigma},$

where we have used $h_{j,y}^{-1}\leq C\varepsilon^{-1}N$ from Lemma 1. The
terms $\mathcal{P}_{2}E_{2}$, $\mathcal{B}_{1}E_{1}$ and
$\mathcal{B}_{2}E_{2}$ can be analyzed in a similar way.

###### Lemma 4

Let Assumptions 1 and 2 hold true. Let $E^{I}_{12}$ denote the Lagrange
interpolant of $E_{12}$ on the Bakhvalov-type mesh $\mathcal{T}_{N}$. Then
there exists a constant $C$ such that the following interpolation error
estimates hold true:

$\displaystyle\|E_{12}-E^{I}_{12}\|_{\Omega\setminus\tau_{N/2-1,N/2-1}}\leq
C\varepsilon N^{-k}+C\varepsilon^{1/2}N^{-(k+1)},$
$\displaystyle\|E_{12}-E^{I}_{12}\|\leq C\varepsilon
N^{-k}+C\varepsilon^{1/2}N^{-(k+1)}+CN^{-1-2\sigma},$
$\displaystyle\|E_{12}-E^{I}_{12}\|_{\varepsilon}\leq CN^{-k},$
$\displaystyle\|\mathcal{P}_{12}E_{12}\|_{\varepsilon}\leq CN^{-1/2-2\sigma},$

where $\mathcal{P}_{12}E_{12}$ is defined in (19).

###### Proof

Consider

$\displaystyle\|E_{12}-E_{12}^{I}\|^{2}=$
$\displaystyle\|E_{12}-E_{12}^{I}\|^{2}_{\Omega_{11}}+\|E_{12}-E_{12}^{I}\|^{2}_{\Omega_{12}\cup\Omega_{21}\cup(\Omega_{22}\setminus\tau_{N/2-1,N/2-1})}$
(32) $\displaystyle+\|E_{12}-E_{12}^{I}\|^{2}_{\tau_{N/2-1,N/2-1}}.$

Similar to (22), the estimation for $\|E_{12}-E_{12}^{I}\|^{2}_{\Omega_{11}}$
is as follows:

$\|E_{12}-E_{12}^{I}\|^{2}_{\Omega_{11}}\leq C\varepsilon^{2}N^{-2k}.$ (33)

To analyze the second term in (32), we frequently use the following estimation

$\|E_{12}^{I}\|^{2}_{\tau_{i,j}}\leq\|E_{12}^{I}\|^{2}_{\infty,\tau_{i,j}}h_{i,x}h_{j,y}\leq
Ce^{-2(\beta_{1}x_{i}+\beta_{2}y_{j})/\varepsilon}h_{i,x}h_{j,y}.$ (34)

From (3d), direct calculations yield

$\|E_{12}\|_{\Omega\setminus\Omega_{11}}\leq C\varepsilon N^{-\sigma}.$ (35)

From (34) and Lemma 2, one has

$\displaystyle\|E_{12}^{I}\|^{2}_{\Omega_{12}}=$
$\displaystyle\sum_{i=N/2-1}^{N-1}\sum_{j=0}^{N/2-2}\|E_{12}^{I}\|^{2}_{\tau_{i,j}}\leq
C\sum_{i=N/2-1}^{N-1}\sum_{j=0}^{N/2-2}e^{-2(\beta_{1}x_{i}+\beta_{2}y_{j})/\varepsilon}h_{i,x}h_{j,y}$
(36) $\displaystyle\leq$ $\displaystyle
C\sum_{i=N/2-1}^{N-1}e^{-2\beta_{1}x_{i}/\varepsilon}h_{i,x}\;\sum_{j=0}^{N/2-2}e^{-2\beta_{2}y_{j}/\varepsilon}h_{j,y}\leq
C\varepsilon N^{-2\sigma},$

where Lemma 1 has been used. Similarly, we have

$\|E_{12}^{I}\|^{2}_{\Omega_{21}}\leq C\varepsilon N^{-2\sigma}.$ (37)

Using (34) again, we have

$\displaystyle\|E_{12}^{I}\|^{2}_{(\Omega_{22}\setminus\tau_{N/2-1,N/2-1})}=\sum_{j=N/2}^{N-1}\|E_{12}^{I}\|^{2}_{\tau_{N/2-1,j}}+\sum_{i=N/2}^{N-1}\sum_{j=N/2-1}^{N-1}\|E_{12}^{I}\|^{2}_{\tau_{i,j}}$
(38) $\displaystyle\leq$
$\displaystyle\sum_{j=N/2}^{N-1}e^{-2(\beta_{1}x_{N/2-1}+\beta_{2}y_{j})/\varepsilon}h_{N/2-1,x}h_{j,y}+\sum_{i=N/2}^{N-1}\sum_{j=N/2-1}^{N-1}e^{-2(\beta_{1}x_{i}+\beta_{2}y_{j})/\varepsilon}h_{i,x}h_{j,y}$
$\displaystyle\leq$ $\displaystyle
CN^{-1-2\sigma}\varepsilon^{2\sigma}+CN^{-2\sigma}\varepsilon^{2\sigma}$

and

$\|E_{12}^{I}\|^{2}_{\tau_{N/2-1,N/2-1}}\leq CN^{-2-4\sigma}.$ (39)

Then the triangle inequality and (33), (35)–(39) yield the first and the
second bounds.

Now we consider $\|(E_{12}-E^{I}_{12})_{x}\|$. By similar derivations for
(27), we have

$\|(E_{12}-E^{I}_{12})_{x}\|_{\Omega_{11}}\leq CN^{-k}.$ (40)

Similar to (28), the following bound can be obtained

$\|(E_{12}-E^{I}_{12})_{x}\|_{\Omega\setminus\Omega_{11}}\leq
C\varepsilon^{-1/2}N^{1-\sigma}.$ (41)

Combing (40) and (41), we prove

$\|(E_{12}-E^{I}_{12})_{x}\|\leq C\varepsilon^{-1/2}N^{1-\sigma}$

and the same bound for $\|(E_{12}-E^{I}_{12})_{y}\|$. Thus we prove the third
bound. Similar to (31), the final bound can be proved.

###### Lemma 5

Let Assumptions 1 and 2 hold true. Let $S^{I}$ and $u^{I}$ denote the Lagrange
interpolants of $S$ and $u$ on the Bakhvalov-type mesh $\mathcal{T}_{N}$,
respectively. Let $\Pi u$ and $\pi_{i}E_{i}$, $i=1,2,12$, be defined in (15)
and (16), respectively. Then there exists a constant $C$ such that the
following interpolation error estimates hold true:

$\displaystyle\sum_{i=1,2,12}\|\pi_{i}E_{i}-E_{i}\|\leq CN^{-(k+1)},$
$\displaystyle\|\nabla(S-S^{I})\|+\|u-u^{I}\|_{\varepsilon}+\|u-\Pi
u\|_{\varepsilon}\leq CN^{-k}.$

###### Proof

Check the derivations in (31) and one finds that $\|\mathcal{P}_{1}E_{1}\|\leq
CN^{-1/2-\sigma}$. Thus

$\|\pi_{1}E_{1}-E_{1}\|\leq\|E_{1}^{I}-E_{1}\|+\|\mathcal{P}_{1}E_{1}\|\leq
CN^{-(k+1)}.$

Similarly, we can prove the bounds for $E_{2}$ and $E_{12}$.

From Lemma 2 and (3b), we prove $\|\nabla(S-S^{I})\|\leq CN^{-k}$ and
$\|S-S^{I}\|_{\varepsilon}\leq C(\varepsilon^{1/2}+N^{-1})N^{-k}$ easily. Then
(3a), the triangle inequality, Lemmas 3 and 4 yield
$\|u-u^{I}\|_{\varepsilon}\leq CN^{-k}$. Besides, from (21), Lemmas 3 and 4
one has $\|u-\Pi u\|_{\varepsilon}\leq CN^{-k}$.

##  4 Uniform convergence

Set $\chi:=\Pi u-u^{N}$. From (13), (14), (3a), (15) and integration by parts,
one has

$\begin{split}&\alpha\|\chi\|_{\varepsilon}^{2}\leq a(\chi,\chi)=a(\Pi
u-u,\chi)\\\ =&\varepsilon\int_{\Omega}\nabla(\Pi
u-u)\nabla\chi\mathrm{d}x\mathrm{d}y+\sum_{i=1,2}\int_{\Omega}(E^{I}_{i}-\mathcal{P}_{i}E_{i}-E_{i})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\\\
&+\int_{\Omega}(E^{I}_{12}-\mathcal{P}_{12}E_{12}-E_{12})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y-\int_{\Omega}\boldsymbol{b}\cdot\nabla(S^{I}-S)\;\chi\mathrm{d}x\mathrm{d}y\\\
&+\sum_{i=1,2,12}\int_{\Omega}(\nabla\cdot\boldsymbol{b})(\pi_{i}E_{i}-E_{i})\;\chi\mathrm{d}x\mathrm{d}y+\int_{\Omega}c(\Pi
u-u)\chi\mathrm{d}x\mathrm{d}y\\\
&+\sum_{i=1,2}\int_{\Omega}(-\mathcal{B}_{i}E_{i})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\\\
=:&\mathrm{I}+\mathrm{II}+\mathrm{III}+\mathrm{IV}+\mathrm{V}+\mathrm{VI}+\mathrm{VII}.\end{split}$
(42)

Now we analyze the terms on the right-hand side of (42). The Cauchy-Schwarz
inequality and Lemma 5 yield

$\displaystyle(\mathrm{I}+\mathrm{VI})+(\mathrm{IV}+\mathrm{V})$ (43)
$\displaystyle\leq$ $\displaystyle C\|\Pi
u-u\|_{\varepsilon}\|\chi\|_{\varepsilon}+C\big{(}\|\nabla(S^{I}-S)\|+\sum_{i=1,2,12}\|\pi_{i}E_{i}-E_{i}\|\big{)}\|\chi\|$
$\displaystyle\leq$ $\displaystyle CN^{-k}\|\chi\|_{\varepsilon}.$

We put the arguments for $\mathrm{II}$, $\mathrm{III}$ and $\mathrm{VII}$ in
the following three lemmas.

###### Lemma 6

Let Assumptions 1 and 2 hold true. Let $\pi_{i}E_{i}$ with $i=1,2$ be defined
in (16). Then one has

$|\mathrm{II}|=\left|\sum_{i=1,2}\int_{\Omega}(E^{I}_{i}-\mathcal{P}_{i}E_{i}-E_{i})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\right|\leq
CN^{-(k+1/2)}\|\chi\|_{\varepsilon}.$ (44)

###### Proof

For the term $\mathrm{II}$, we just consider $E_{1}$ since we can analyze
$E_{2}$ in a similar way. Set $D_{0}:=[x_{N/2-1},x_{N/2}]\times[0,1]$.
According to (16) and (17), one has

$\displaystyle\int_{\Omega}(E^{I}_{1}-\mathcal{P}_{1}E_{1})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y$
(45) $\displaystyle=$
$\displaystyle\left(E^{I}_{1}-\mathcal{P}_{1}E_{1}-E_{1},\boldsymbol{b}\cdot\nabla\chi\right)_{\Omega\setminus
D_{0}}+\left(E^{I}_{1}-\mathcal{P}_{1}E_{1}-E_{1},\boldsymbol{b}\cdot\nabla\chi\right)_{D_{0}}$
$\displaystyle=$
$\displaystyle\left(E_{1}^{I}-E_{1},\boldsymbol{b}\cdot\nabla\chi\right)_{\Omega\setminus
D_{0}}+\left(-\mathcal{F}_{1},\boldsymbol{b}\cdot\nabla\chi\right)_{\Omega\setminus
D_{0}}$
$\displaystyle+\left(\mathcal{F}_{2}-E_{1},\boldsymbol{b}\cdot\nabla\chi\right)_{D_{0}}$
$\displaystyle=:$ $\displaystyle\mathrm{T}_{1}+\mathrm{T}_{2}+\mathrm{T}_{3},$

where $(E^{I}_{1}-\mathcal{P}_{1}E_{1})|_{\Omega\setminus
D_{0}}=(E_{1}^{I}-\mathcal{F}_{1})|_{\Omega\setminus D_{0}}$,
$(E^{I}_{1}-\mathcal{P}_{1}E_{1})|_{D_{0}}=\mathcal{F}_{2}|_{D_{0}}$ and

$\displaystyle\mathcal{F}_{1}:=$
$\displaystyle\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}E_{1}(x_{N/2-1},y_{j}^{t})\theta_{N/2-1,j}^{0,t}+E_{1}(x_{N/2-1},y_{N})\theta_{N/2-1,N}^{0,0},$
$\displaystyle\mathcal{F}_{2}:=$
$\displaystyle\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}E_{1}(x_{N/2},y_{j}^{t})\theta_{N/2,j}^{0,t}+E_{1}(x_{N/2},y_{N})\theta_{N/2,N}^{0,0}.$

From Lemma 3 and the Cauchy-Schwarz inequality, one has

$\displaystyle|\mathrm{T}_{1}|\leq$ $\displaystyle
C\|E_{1}-E^{I}_{1}\|_{\Omega\setminus D_{0}}\|\nabla\chi\|_{\Omega\setminus
D_{0}}\leq C\varepsilon^{1/2}N^{-(k+1/2)}\|\nabla\chi\|\leq
CN^{-(k+1/2)}\|\chi\|_{\varepsilon}.$ (46)

Note $|E_{1}(x_{N/2-1},y_{j}^{t})|\leq CN^{-\sigma}$ for any $j,t$ and
$h_{N/2-2,x}\leq C\varepsilon$. The Cauchy-Schwarz inequality yields

$\displaystyle|\mathrm{T}_{2}|\leq$ $\displaystyle
C\|\mathcal{F}_{1}\|_{D_{1}}\|\nabla\chi\|_{D_{1}}$ (47) $\displaystyle\leq$
$\displaystyle
CN^{-\sigma}\|\nabla\chi\|_{D_{1}}\left(\sum_{j=0}^{N-1}\sum_{t=0}^{k-1}\|\theta_{N/2-1,j}^{0,t}\|_{D_{1}}+\|\theta_{N/2-1,N}^{0,0}\|_{D_{1}}\right)$
$\displaystyle\leq$ $\displaystyle
CN^{-\sigma}\|\nabla\chi\|_{D_{1}}\sum_{j=0}^{N-1}h^{1/2}_{N/2-2,x}h^{1/2}_{j,y}$
$\displaystyle\leq$ $\displaystyle CN^{1/2-\sigma}\|\chi\|_{\varepsilon},$

where $D_{1}:=[x_{N/2-2},x_{N/2-1}]\times[0,1]$ and we have made use of the
supports of hat functions $\theta^{s,t}_{i,j}$.

Now we deal with the term $\mathrm{T}_{3}$. Note
$|E_{1}(x_{N/2},y_{j}^{t})|\leq C\varepsilon^{\sigma}$ for any $j,t$ and
$h_{N/2-1,x}\leq CN^{-1}$. The Cauchy-Schwarz inequality yields

$\displaystyle|\mathrm{T}_{3}|\leq$ $\displaystyle
C(\|E_{1}\|_{D_{0}}+\|\mathcal{F}_{2}\|_{D_{0}})\|\nabla\chi\|_{D_{0}}$ (48)
$\displaystyle\leq$ $\displaystyle
C\left(\varepsilon^{1/2}N^{-\sigma}+\varepsilon^{\sigma}\sum_{j=0}^{N-1}h^{1/2}_{N/2-1,x}h^{1/2}_{j,y}\right)\|\nabla\chi\|_{D_{2}}$
$\displaystyle\leq$ $\displaystyle
C(N^{-\sigma}+\varepsilon^{\sigma-1/2})\|\chi\|_{\varepsilon}$

where direct calculations and (3c) yield $\|E_{1}\|_{D_{0}}\leq
C\varepsilon^{1/2}N^{-\sigma}$ and $\|\mathcal{F}_{2}\|_{D_{0}}$ is analyzed
in a similar way to $\|\mathcal{F}_{1}\|_{D_{1}}$ in (47).

Substituting (46)–(48) into (45), we are done.

###### Lemma 7

Let Assumptions 1 and 2 hold true. Let $\pi_{12}E_{12}$ be defined in (16).
Then one has

$|\mathrm{III}|=\left|\int_{\Omega}(\pi_{12}E_{12}-E_{12})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\right|\leq
CN^{-(k+1/2)}\|\chi\|_{\varepsilon}.$ (49)

###### Proof

Set $D_{3}=:\Omega\setminus\tau_{N/2-1,N/2-1}$. According to (16), one has

$\displaystyle\int_{\Omega}(\pi_{12}E_{12}-E_{12})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y$
(50) $\displaystyle=$
$\displaystyle\left(\pi_{12}E_{12}-E_{12},\boldsymbol{b}\cdot\nabla\chi\right)_{D_{3}}+\left(\pi_{12}E_{12}-E_{12},\boldsymbol{b}\cdot\nabla\chi\right)_{\tau_{N/2-1,N/2-1}}$
$\displaystyle=$
$\displaystyle\left(E_{12}^{I}-E_{12},\boldsymbol{b}\cdot\nabla\chi\right)_{D_{3}}-\left(\mathcal{F}_{3},\boldsymbol{b}\cdot\nabla\chi\right)_{D_{3}}+\left(\mathcal{F}_{4}-E_{12},\boldsymbol{b}\cdot\nabla\chi\right)_{\tau_{N/2-1,N/2-1}}$
$\displaystyle=:$ $\displaystyle\mathrm{Q}_{1}+\mathrm{Q}_{2}+\mathrm{Q}_{3},$

where $(\pi_{12}E_{12})|_{D_{3}}=(E_{12}^{I}-\mathcal{F}_{3})|_{D_{3}}$,
$(\pi_{12}E_{12})|_{\tau_{N/2-1,N/2-1}}=\mathcal{F}_{4}|_{\tau_{N/2-1,N/2-1}}$
and

$\displaystyle\mathcal{F}_{3}:=$
$\displaystyle\sum_{t=1}^{k-1}E_{12}(x_{N/2-1},y_{N/2-1}^{t})\theta_{N/2-1,N/2-1}^{0,t}+E_{12}(x_{N/2-1},y_{N/2-1})\theta_{N/2-1,N/2-1}^{0,0}$
$\displaystyle+\sum_{s=1}^{k-1}E_{12}(x_{N/2-1}^{s},y_{N/2-1})\theta_{N/2-1,N/2-1}^{s,0},$
$\displaystyle\mathcal{F}_{4}:=$
$\displaystyle\sum_{t=0}^{k-1}E_{12}(x_{N/2},y_{N/2-1}^{t})\theta_{N/2,N/2-1}^{0,t}+E_{12}(x_{N/2},y_{N/2})\theta_{N/2,N/2}^{0,0}$
$\displaystyle+\sum_{s=0}^{k-1}E_{12}(x_{N/2-1}^{s},y_{N/2})\theta_{N/2-1,N/2}^{s,0}.$

From Lemma 4 and the Cauchy-Schwarz inequality, one has

$\displaystyle|\mathrm{Q}_{1}|\leq
C\|E_{12}^{I}-E_{12}\|_{D_{3}}\|\nabla\chi\|_{D_{3}}\leq
C\varepsilon^{1/2}N^{-(k+1/2)}\|\nabla\chi\|\leq
CN^{-(k+1/2)}\|\chi\|_{\varepsilon}.$ (51)

Note $|E_{12}(x_{N/2-1}^{s},y_{N/2-1}^{t})|\leq CN^{-2\sigma}$ for any $0\leq
s,t\leq k-1$ and $h_{N/2-2,x},h_{N/2-2,y}\leq C\varepsilon$. Then one has

$\displaystyle\|\mathcal{F}_{3}\|_{D_{3}}\leq$ $\displaystyle
CN^{-2\sigma}\left(\sum_{t=1}^{k-1}\|\theta_{N/2-1,N/2-1}^{0,t}\|_{D_{3}}+\|\theta_{N/2-1,N/2-1}^{0,0}\|_{D_{3}}+\sum_{s=1}^{k-1}\|\theta_{N/2-1,N/2-1}^{s,0}\|_{D_{3}}\right)$
$\displaystyle\leq$ $\displaystyle
CN^{-2\sigma}\left(h_{N/2-1,x}^{1/2}h_{N/2-2,y}^{1/2}+h_{N/2-2,x}^{1/2}h_{N/2-2,y}^{1/2}+h_{N/2-2,x}^{1/2}h_{N/2-1,y}^{1/2}\right)$
$\displaystyle\leq$ $\displaystyle C\varepsilon^{1/2}N^{-1/2-2\sigma}.$

Thus the Cauchy-Schwarz inequality yields

$\displaystyle|\mathrm{Q}_{2}|\leq
C\|\mathcal{F}_{3}\|_{D_{3}}\|\nabla\chi\|_{D_{3}}\leq
CN^{-1/2-2\sigma}\|\chi\|_{\varepsilon}.$ (52)

Note
$|E_{12}(x_{N/2}^{s},y_{N/2-1}^{t})|+|E_{12}(x_{N/2-1}^{s},y_{N/2}^{t})|+|E_{12}(x_{N/2},y_{N/2})|\leq
C\varepsilon^{\sigma}N^{-\sigma}$ for any $0\leq s,t\leq k-1$. The Cauchy-
Schwarz inequality yields

$\displaystyle|\mathrm{Q}_{3}|\leq$ $\displaystyle
C(\|\mathcal{F}_{4}\|_{\tau_{N/2-1,N/2-1}}+\|E_{12}\|_{\tau_{N/2-1,N/2-1}})\|\nabla\chi\|_{\tau_{N/2-1,N/2-1}}$
(53) $\displaystyle\leq$ $\displaystyle
C\left(\varepsilon^{\sigma}N^{-\sigma}h^{1/2}_{N/2-1,x}h^{1/2}_{N/2-1,y}+\varepsilon
N^{-2\sigma}\right)\|\nabla\chi\|_{\tau_{N/2-1,N/2-1}}$ $\displaystyle\leq$
$\displaystyle C\varepsilon^{1/2}N^{-2\sigma}\|\chi\|_{\varepsilon}.$

Substituting (51)–(53) into (50), we are done.

###### Lemma 8

Let Assumptions 1 and 2 hold true. Let $\mathcal{B}_{i}E_{i}$ with $i=1,2$ be
defined in (16). Then one has

$|\mathrm{VII}|=\left|\sum_{i=1,2}\int_{\Omega}(-\mathcal{B}_{i}E_{i})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\right|\leq
CN^{-k}R(N,\varepsilon)\|\chi\|_{\varepsilon},$ (54)

where $R(N,\varepsilon)=N^{-3/2}|\ln(\varepsilon N)|^{1/2}$.

###### Proof

We just present the analysis on the term involved with $\mathcal{B}_{1}E_{1}$
and the other could be analyzed in a similar way. Hölder inequalities

$\displaystyle\left|\int_{\Omega}(-\mathcal{B}_{1}E_{1})\;\boldsymbol{b}\cdot\nabla\chi\mathrm{d}x\mathrm{d}y\right|$
$\displaystyle\leq$
$\displaystyle\sum_{j=0,N-1}\|\mathcal{B}_{1}E_{1}\|_{\infty,\tau_{N/2-1,j}}\|\boldsymbol{b}\cdot\nabla\chi\|_{1,\tau_{N/2-1,j}}$
$\displaystyle\leq$ $\displaystyle
C\sum_{j=0,N-1}\left(N^{-\sigma}h_{N/2-1,x}^{1/2}h_{j,y}^{1/2}\|\boldsymbol{b}\cdot\nabla\chi\|_{\tau_{N/2-1,j}}\right)$
$\displaystyle\leq$ $\displaystyle
C(N^{-(1+\sigma)}+N^{-(1/2+\sigma)}|\ln(\varepsilon
N)|^{1/2})\|\chi\|_{\varepsilon}$ $\displaystyle\leq$ $\displaystyle
CN^{-(1/2+\sigma)}|\ln(\varepsilon N)|^{1/2}\|\chi\|_{\varepsilon},$

where we have used $h_{0,y}\leq C\varepsilon N^{-1}$ in Lemma 1 and
$h_{N/2-1,x}=-\frac{\sigma\varepsilon}{\beta}\ln(\varepsilon N)$ from (4).

###### Remark 3

Practically $R(N,\varepsilon)$ is bounded: If we assume $N\geq 10$ and
$\varepsilon\geq 10^{-1001}$, then

$R(N,\varepsilon)\leq\sqrt{\ln 10}.$

Now we are in a position to present the main result.

###### Theorem 4.1

Let Assumptions 1 and 2 hold true. Let $u$ and $u^{N}$ be the solutions of (1)
and (12), respectively. Then one has

$\displaystyle\|u-u^{N}\|_{\varepsilon}\leq$ $\displaystyle
CN^{-k}R(N,\varepsilon),$

where $R(N,\varepsilon)=N^{-3/2}|\ln(\varepsilon N)|^{1/2}$.

###### Proof

Substituting (43), (44), (49) and (54) into (42), we obtain $\|\Pi
u-u^{N}\|_{\varepsilon}\leq CN^{-k}R(N,\varepsilon)$. From a triangle
inequality and Lemma 5, one has

$\|u-u^{N}\|_{\varepsilon}\leq\|u-\Pi u\|_{\varepsilon}+\|\Pi
u-u^{N}\|_{\varepsilon}\leq CN^{-k}R(N,\varepsilon).$

Thus we are done.

###### Remark 4

In Theorem 4.1, we show the almost optimal estimation, which is slightly
better than the bound in (Roos1Scho2:2012-Analysis, , Theorem 1).

##  5 Numerical Experiments

In this section we present numerical experiments that support our theoretical
results. All calculations were carried out using Intel Visual Fortran 11 and
the discrete problems were solved by the nonsymmetric iterative solver GMRES;
see, e.g., Ben1Gol2Lie3:2005-Numerical .

For our experiments we used the boundary value problem

$\displaystyle-\varepsilon\Delta u-(2+2x-y)u_{x}-(3-x+2y)u_{y}+u$
$\displaystyle=f(x,y)\quad$ $\displaystyle\text{in $\Omega=(0,1)^{2}$},$
$\displaystyle u$ $\displaystyle=0\quad$ on $\partial\Omega$,

where the right-hand side $f$ is chosen such that

$u(x,y)=2\sin(\pi
x)\left(1-e^{-\frac{2x}{\varepsilon}}\right)(1-y)^{2}\left(1-e^{-\frac{y}{\varepsilon}}\right)$

is the exact solution. This solution exhibits typical exponential layer
behaviour as described in Assumption 1.

Numerical results are presented in Tables 1 and 2 and Figure 2, which support
our main result. Tables 1 and 2 list errors in the energy norm, i.e.,
$\|u-u^{N}\|_{\varepsilon}$, for $\varepsilon=10^{-4},10^{-5},\ldots,10^{-8}$
and $N=8,16,32,64,128,256$, in the cases of $k=1,2$. These data show uniform
convergence with respect to the singular perturbation parameter $\varepsilon$.
In the cases of $k=3,4$, errors and convergence orders also show uniform
convergence, which are plotted in Figure 2.

Here we also compare two meshes, which are Bakhvalov-type mesh (4) (denoted by
B-mesh) and Bakhvalov-Shishkin mesh (denoted by B-S-mesh). B-S-mesh
Linb:2000-Analysis  yields an optimal convergence order for bilinear element
and has excellent performances. According to Linb:2000-Analysis , B-S-mesh is
defined by tensor product and the mesh points in $x-$direction are defined by

$\varphi(i/N)=\left\\{\begin{split}&-\frac{(k+1)\varepsilon}{\beta_{1}}\ln\left(\frac{N^{2}-2i(N-1)}{N^{2}}\right)\quad&&\text{for
$i=0,\ldots,N/2$},\\\
&\varphi(1/2)+2(1-\varphi(1/2))/(i/N-1/2)\quad&&\text{for
$i=N/2+1,\ldots,N$}.\end{split}\right.$

The mesh points along $y-$direction can be defined similarly.

To compare B-mesh and B-S-mesh, we plot errors in the energy norm on these
meshes for $\varepsilon=10^{-8}$ and $N=8,16,32,64,128$, for the cases
$k=1,2,3,4$ on log-log chart in Figure 2. We can see that these meshes have
similar performances and B-S-mesh yields slightly smaller errors.

Besides, from numerical experiments we find that the linear systems become
harder to be solved by iterative solvers when $k$ and $N$ become bigger and
$\varepsilon$ becomes smaller.

Table 1: Errors and orders in the energy norm on B-mesh for $k=1$ $\varepsilon$ $N$ | 8 | 16 | 32 | 64 | 128 | 256
---|---|---|---|---|---|---
$10^{-4}$ | 0.227E+0 | 0.109E+0 | 0.540E-1 | 0.269E-1 | 0.135E-1 | 0.673E-2
$10^{-4}$ | 1.06 | 1.02 | 1.00 | 1.00 | 1.00 | —
$10^{-5}$ | 0.228E+0 | 0.109E+0 | 0.540E-1 | 0.269E-1 | 0.135E-1 | 0.673E-2
$10^{-5}$ | 1.06 | 1.02 | 1.00 | 1.00 | 1.00 | —
$10^{-6}$ | 0.229E+0 | 0.109E+0 | 0.540E-1 | 0.269E-1 | 0.135E-1 | 0.673E-2
$10^{-6}$ | 1.07 | 1.02 | 1.00 | 1.00 | 1.00 | —
$10^{-7}$ | 0.231E+0 | 0.110E+0 | 0.541E-1 | 0.269E-1 | 0.135E-1 | 0.673E-2
$10^{-7}$ | 1.08 | 1.02 | 1.00 | 1.00 | 1.00 | —
$10^{-8}$ | 0.234E+0 | 0.110E+0 | 0.541E-1 | 0.269E-1 | 0.135E-1 | 0.673E-2
$10^{-8}$ | 1.09 | 1.02 | 1.01 | 1.00 | 1.00 | —
Table 2: Errors and orders in the energy norm on B-mesh for $k=2$ $\varepsilon$ $N$ | 8 | 16 | 32 | 64 | 128
---|---|---|---|---|---
$10^{-4}$ | 0.502E-1 | 0.105E-1 | 0.247E-2 | 0.603E-3 | 0.150E-3
$10^{-4}$ | 2.25 | 2.10 | 2.03 | 2.01 | —
$10^{-5}$ | 0.501E-1 | 0.105E-1 | 0.247E-2 | 0.604E-3 | 0.150E-3
$10^{-5}$ | 2.25 | 2.09 | 2.03 | 2.01 | —
$10^{-6}$ | 0.500E-1 | 0.105E-1 | 0.246E-2 | 0.604E-3 | 0.150E-3
$10^{-6}$ | 2.25 | 2.09 | 2.03 | 2.01 | —
$10^{-7}$ | 0.500E-1 | 0.105E-1 | 0.246E-2 | 0.604E-3 | 0.150E-3
$10^{-7}$ | 2.25 | 2.09 | 2.03 | 2.01 | —
$10^{-8}$ | 0.502E-1 | 0.105E-1 | 0.246E-2 | 0.604E-3 | 0.150E-3
$10^{-8}$ | 2.25 | 2.10 | 2.03 | 2.01 | —
Figure 2: Error: energy norm, $\varepsilon=10^{-8}$.

##  6 Bibliography

## References

  * [1] T. Apel.  Anisotropic Finite Elements: Local Estimates and Applications.  Advances in Numerical Mathematics. B. G. Teubner, Stuttgart, 1999.
  * [2] N. S. Bahvalov.  On the optimization of the methods for solving boundary value problems in the presence of a boundary layer.  Zh. Vychisl. Mat. Mat. Fiz., 9:841–859, 1969.
  * [3] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numer., 14:1–137, 2005.
  * [4] P. G. Ciarlet.  The Finite Element Method for Elliptic Problems.  Studies in Mathematics and its Applications. North-Holland, Amsterdam, 1978.
  * [5] S. Franz, T. Linß, H.-G. Roos, and S. Schiller.  Uniform superconvergence of a finite element method with edge stabilization for convection-diffusion problems.  J. Comput. Math., 28(1):32–44, 2010.
  * [6] T. Linß.  Analysis of a Galerkin finite element method on a Bakhvalov-Shishkin mesh for a linear convection-diffusion problem.  IMA J. Numer. Anal., 20(4):621–632, 2000.
  * [7] T. Linß.  Layer-adapted meshes for reaction-convection-diffusion problems, volume 1985 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin, 2010.
  * [8] T. Linß and M. Stynes.  Asymptotic analysis and Shishkin-type decomposition for an elliptic convection-diffusion problem.  J. Math. Anal. Appl., 261(2):604–632, 2001.
  * [9] J. J. H. Miller, E. O’Riordan, and G. I. Shishkin.  Fitted Numerical Methods for Singular Perturbation Problems.  World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, revised edition, 2012.  Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimensions.
  * [10] H.-G. Roos.  Error estimates for linear finite elements on Bakhvalov-type meshes.  Appl. Math., 51(1):63–72, 2006.
  * [11] H.-G. Roos and M. Schopf.  Analysis of finite element methods on Bakhvalov-type meshes for linear convection-diffusion problems in 2D.  Appl. Math., 57(2):97–108, 2012.
  * [12] H.-G. Roos, M. Stynes, and L. Tobiska.  Robust Numerical Methods for Singularly Perturbed Differential Equations, volume 24 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, second edition, 2008.
  * [13] M. Stynes.  Steady-state convection-diffusion problems.  Acta Numer., 14:445–508, 2005.
  * [14] J. Zhang and X. Liu.  Optimal order of uniform convergence for finite element method on Bakhvalov-type meshes.  J. Sci. Comput., 85(1):2, 2020.


# Impacts of the Space Technology Evolution in the V&V of Embedded Software-
Intensive Systems

Carlos L. G. Batista1, Tania Basso2, Fátima Mattiello-Francisco1, Regina
Moraes2 3  1 National Institute for Space Research (INPE), São José dos
Campos, Brazil
carlos.batista, fatima.mattiello<EMAIL_ADDRESS>2 University of Campinas
(UNICAMP), Limeira, Brazil
taniabasso, regina<EMAIL_ADDRESS>3 University of Coimbra (UC), Coimbra,
Portugal - Full / Regular Research Paper (CSCI - ISSE)

###### Abstract

CubeSat-based nanosatellites are composed of COTS components and rely on its
structure and standardized interfaces. A challenge in the nanosatellites
context is to adapt the V&V (Verification and Validation) process to answer to
the increase importance of the embedded software, to reduce the artefacts to
be delivered aiming at cutting cost and time and still complying with
international standards. This work presents an analysis of the strategy
adopted in a real nanosatellite for the development of the OBDH software
embedded in NanosatC-BR2 mission. The goal is to discuss the impact that the
standardization of the structure and interfaces of the CubeSat impose on the
V&V process of the SiS and to highlight the challenges of “New Space Age” for
the use of existing V&V techniques and methods.

###### Index Terms:

CubeSat, nanosatellites, V&V techniques

##  I Introduction

In the last decade the space area has evolved significantly, thanks for the
advent of microprocessors and the spread of software as a service that allowed
the development of micro and nanosatellites programs. The advent of the
CubeSat standard has made the development cycle of small satellite projects
(nanosatellites) much faster and cheaper, compared to larger space missions.
In addition to the use of COTS (components off-the-shelf), that reduce cost,
the standardization of the satellite structure and its interfaces helps
improving the architectural satellite solutions based on CubeSat [1].

With the technological evolution of embedded electronics, memories and
satellites processors, the functionalities implemented by software have
increased. The so-called software-intensive systems (SiS) embedded in
satellites are increasing in terms of complexity and becoming critical
elements for the mission. Moreover, the “new space age”, i.e., greater access
to space by different segments of the society, demands for new space products
and services, requiring regulation and well-established development processes.
CubeSat can be taken by an example due to its popularity within the academia
and industry. The complexity of CubeSat structure, cabling and interfaces has
been significantly simplified with the standardization, but the software has
not.

Addressing this problem, SiS suppliers are reducing both the design artifacts
to be delivered in the reviews baseline and the number of reviews in the
development cycle of many SiS embedded in nanosatellites. However, negative
side effects of this strategy in the quality of the nanosatellite V&V process
(unknown software requirements in the begin of the development cycle,
immaturity of the stakeholders and expertise of the development time in the
space area projects) have been observed regarding the lack of compliance with
international standards. The challenge is to establish effective V&V
approaches that can contribute to make the nanosatellites V&V process shorter
and cheaper but mainly compliant with the international standards.

This work presents the V&V strategy adopted in the development of NanosatC-
Br2, a scientific nanosatellite leading by the Brazilian Institute for Space
Research (INPE) with partners in Brazilian Universities and software private
supplier. Considered the main SiS element of the mission, and the focus of
this study, the On-Board Data Handling (OBDH) software is responsible to the
implementation of the interfaces with the mission payload and the ground
systems. The approach addresses the V&V activities with the respective
techniques and methodologies adopted by INPE in the development of SiS
projects embedded in satellites based on the guidelines of the European
Cooperation for Space Standardization (ECSS) [2]The V&V activities rely on
baseline reviews and Model-Based Testing (MBT) to generate, automatically,
test cases, being able to systematize the OBDH testing.

In addition to the introduction, Section II presents the background and
related work. Section III presents the V&V process adopted by INPE in the
development of OBDH SiS. Section IV presents research opportunities to improve
the V&V process of standard Cubesat satellites, ending with the conclusions
and future work in Section V.

##  II Background and related work

The CubeSat standard emerged from an effort by California Polytechnic State
University and Stanford University in 1999. It is a simple and small satellite
model, developed with low-cost (less than one million Euro) and a short
project cycle (less than two years to flight readiness) educational
objectives, both in development and in launch [1]. It is a cube-shaped
nanosatellite, whose edges measure 10 centimeters, volume of one liter and can
hold a bus and payload with mass of 1.3 kilograms. More than one unit (or 1U)
can be combined to form larger satellites (e.g. 2U, 8U) and can use ejection
systems in standardized orbits (e.g. P-POD - Poly Picosatellite Orbital
Deployer), allowing to release several satellites through the same interface
and to launch several CubeSat in a single launch vehicle or share trips with
supplies in the refueling operations from ISS (International Space Station).

The architecture of a satellite mission is composed of four parts: the space,
the launch, the ground and the user segments. The space segment is the more
important one for this work. Within it are the payload (set of equipment
dedicated to the application of the space mission) and the platform (set of
subsystems designed to support the operation of the mission in orbit, also
referred as bus). Figure 1 shows the typical architecture of a satellite
platform, with the space segment subsystems and its communication with the
ground segment. We highlight the OBDH (On-Board Data Handling) subsystem, a
typical SiS, which is responsible for monitoring the health of the other
subsystems of the satellite platform, for interacting with the payloads of the
mission and for ensuring communication with the ground in the operation of the
mission.

Figure 1: Architecture of the space and ground segments

Besides dataprocessing, OBDH plays the role of the spacecraft command &
control and telemetry data management.

The command & control functions on satellites enable remote operation of the
satellite subsystems, allowing the ground controller to change the behavior of
the spacecraft as, for example, change the operational mode and coordinate
orbital maneuvers.

The functions associated with telemetry on board include real-time monitoring
of the platform’s equipment and subsystems and also of the payload in terms of
critical parameters to the mission’s survival. The data packaging and
temporary storage of mission on board are also functions performed by the OBDH
software, responsible to transmit those data to the ground in the window of
visibility of the satellite and Ground Stations. The platform telemetry
generally encompasses two categories of data: housekeeping (parameters needed
to check the health and state of the spacecraft) and attitude parameters. The
payload telemetry consists of the application’s data set (the mission’s data).
Computing on board of satellites is also used to reduce the telemetry, such as
compressing data on board to reduce the amount of data and reduce the needs of
high downlink bit rates.

Software engineering is normally addressed in the space project standards.
Best practices are recommended by NASA and ESA in cooperation with the
industry [2][3] . In fact, standardization is important to reduce risks, costs
and improve quality and communication in the space for any spacecraft mission,
including CubeSat [1]. Scholz [4] presents a standard handbook reviewing
existing space standards and their potential application within CubeSat
mission.

Focused on the V&V process, the interests of this work are on techniques and
tools to assure the quality of the software embedded in the spacecraft. The
most frequent techniques being applied to improve and assure quality in space
software are reviews and tests. According to ECSS standard [5], reviews assure
good quality of documents and are a good opportunity for customers to check,
gradually, the comprehension of the suppliers about the software details. The
tests demonstrate that the implementation meets the requirements and
performance constraints are correctly and completely implemented.

Figure 2 presents the basic mission (column 1) and related software (rightest
column) life cycle activities recommended in ECSS, along the typical space
mission phases: Mission Analysis/ Needs Identification (O), Feasibility (A),
Preliminary Definition (Project and Product)(B), Detailed Definition (Product)
(C), Production/Ground Qualification Testing (D), Utilization (E), Disposal
(F) depicted in the central part of this figure . The darker bars represent
the period that each space mission activity is carried out under system point
of view. The lighter bars correspond to the software activities related to the
software life-cycle development under software subsystem point of view (i.e.,
phases B, C and D).

Figure 2: ECSS Space Mission Project Life Cycle

On the bottom of Figure 2, can be observed the sequence of reviews in the
space mission project, at system engineering level. They are: MDR - Mission
Definition Review, PRR - Preliminary Requirements Review, SRR - System
Requirement Review, PDR - Preliminary Design Review, DDR - Detailed Design
Review, CDR – Critical Design Review, QR - Qualification Review, AR -
Acceptance Review, ORR - Operation Requirements Review, FRR - Flight Readiness
Review. The software project (a subsystem in the hierarchy of space system
development) shall be synchronized with those milestones.

Along the life cycle a series of models must be elected to help the V&V
processes. Three main physical models of the satellite are considered: (i)
Engineering Model (EM) - the development environment based on prototyping with
non-qualified components; (ii) Qualification Model (QM) - the satellite
integrated in the flight structure to be stressed on testing; and (iii) Flight
Model (FM) - the flight structure assembled for flying. From a software
project perspective QM and FM are considered identical. Note the bar
associated with the software verification activities. It considers the
software at CDR is already embedded in the hardware, being the OBDH subsystem
(software and hardware integrated).

In CubeSats, the OBDH grows in complexity as long as the hardware gets
miniaturized. More and more services are delegated to the on board software
replacing functions once done by hardware. More pieces of software are
developed to interface, control and operate different subsystems. These new
features increase the overall complexity of the satellite and also increase
the number of vulnerabilities . The success of the mission relies on the
perfect operation of the OBDH , even more for CubeSats.

Different works have been proposing solutions to mitigate vulnerabilities
without increase the costs of a V&V process. Some focusing on specific pieces
of software [6][7], others new methodologies [8][9], others to develop models
in model-driven approaches [10][11], and, even more, some of them exploit the
standards tailoring [12][13].

##  III V&V process for OBDH developed at INPE

CubeSat establishes in its minimum structure (1U) electrical interfaces with
standardized connectors through which the subsystems of the satellite platform
connect themselves (bus). In this way, the architectural complexity of the
satellite in terms of structure and cabling becomes significantly smaller,
adding the advantage of several COTS subsystems that the industry is able to
supply.

Due to increased demand of nanosatellite mission for the purpose of new
technology demonstration in orbit, in 2016, the European Space Agency (ESA)
released a document entitled Tailored ECSS Engineering Standard for In-Orbit
Demonstration (IOD) CubeSat Projects111https://copernicus-masters.com/wp-
content/uploads/2017/03/IOD_CubeSat_
ECSS_Eng_Tailoring_Iss1_Rev3.pdf. IOD Cubesats are considered structure (1-U
or multiple) generally characterized by the following attributes: (i) Complete
stand-alone systems including platform, payload, ground segment and
operations; (ii) Higher risk acceptance profile; (iii) Low level of complexity
(relative to other space ESA project); (iv) Low cost and short schedule; (v)
Short operational life-cycle (less than one year); (vi) Acceptance of single
point of failures; (vii) Limited redundancy and fault tolerance; (viii) Robust
safe mode; (ix) Extensive use of COTS elements and extensive tests focusing on
system level (functional and environment, qualification and acceptance); (x)
Simple project organization with well integrated teams (few suppliers or
subcontractors).

ESA states that only a few in the list of classical acquisition standards may
be selected as applicable. Thus, “The IOD verification process shall be
adapted to reflect the reduced complexity of work and the absence of one or
more disciplines”. Regarding software discipline, INPE has invested in
reducing complexity of the verification process of OBDH software developed for
being used in two ongoing Cubesat-based missions, named NanosatC-Br2 and
CONASAT-0.

The V&V process is still supported by reviews at design time as recommended in
ECSS. However, in the light of the ECSS Mission Project Life Cycle, already
presented in Figure 2 and fully implemented for the traditional satellites,
CubeSat project has eliminated the revisions DDR - Detailed Design Review and
QR - Qualification Review. The justification is that DDR aims, in traditional
satellites, to verify with the team hired for the development, the consistency
of the design of the drivers and software libraries with the hardware also
under development. In the case of Cubesat nanosatellite, the use of
standardized interfaces and protocols available as third-party components
(COTS) guarantees the necessary compatibility. Similarly, since the basic
Cubesat 1-U platform is already provided as a qualified COTS, QR review that
is justified on traditional satellites to demonstrate that the satellite
structure (QM) with its integrated subsystems has been verified and can be
manufactured for flight (FM) is no more necessary. However, the AR -
Acceptance Review was maintained, as it is a system-level review. Figure 3
shows the adaptation of the ECSS Mission Software Life Cycle carried out at
INPE for nanosatellites. The standardization of hardware benefits the embedded
software design due to the availability, during its development, of COTS on-
board computer architectures where the software should be embedded.

Figure 3: Software Development Life Cycle for Nanosatellites

Thus, with the help of Model Driven Engineering (MDE) approaches, early in the
development cycle, it is possible to support the requirements definition using
techniques that allow to verify the behavior of the embedded software as
expected by means executable models (Model-in-the-loop - MIL) in the PDR
revision and effectively embedded in the target hardware (Hardware-in-the-loop
- HIL) in the CDR revision (see Figure 3). MIL prototyping makes it possible
to anticipate, in the mission development cycle, the identification of SiS
interoperability and robustness problems that might occur on the interaction
of two communicatiing SiS [14] .

###  III-A Deliverables for Mission Reviews

This section presents an example of a typical V&V process tailored by INPE for
the development of an On-Board Data Handling (OBDH) subsystem embedded into
Cubesat-based satellite missions. The process outputs include both a set of
deliverables considered input for each mission review in Figure 3 and the
testing method used for the integration of the OBDH with other subsystems of
the platform, and other SiS such as the satellite payloads.

The deliverables related to each review of the software development life cycle
for Cubesat-based missions are in the third column of Figure I. It is shown in
bold the final delivery version and the respective reviews are the ones in
which the document is considered closed. Requirement Baseline and Interface
Requirements are not closed at the end of SRR review, as in traditional
satellites, but only at the end of CDR. This is feasible in software designing
that follows prototyping and incremental approaches.

TABLE I: OBDH (SiS) Mission Review and Deliverables

Mission Reviews | Deriverables - Traditional Satellite | Deriverables - CubeSats | Deriverable Owner
---|---|---|---
SRR | Requirements Baseline | Requirements Baseline | Customer
Interface Requirements | Interface Requirements | Customer
SW Development Plan |  | Supplier
V&V Plan | V&V Plan | V&V Team
PDR | SW Development Plan | SW Design | Supplier
Technical Specification | Requirements Baseline | Supplier
V&V Plan |  | V&V Team
| Interface Requirements | Customer
| SW Test Plan | Supplier
| SW Test Report | Supplier
CDR | SW Design | SW Design | Supplier
SW Source Code | SW Source Code | Supplier
SW Test Plan | SW Test Plan | Supplier
SW Test Reports | SW Test Reports | Supplier
User Manual | User Manual | Supplier
V&V Plan | V&V Plan | V&V Team
V&V Plan Instrument | System AIT Plan | V&V Team
V&V Specification | Mission Operation Concept | V&V Team
V&V Plan Subsystem |  | V&V Team
V&V Specification Subsystem |  | V&V Team
| Requirements Baseline | Customer
| Interface Requirements | Customer
AR | SW approval | SW approval | Customer
User Manual | User Manual | Supplier
V&V Plan System | System AIT Plan | V&V Team
V&V Specification System | System AIT Report | V&V Team
V&V Report System | Mission Operation Concept | V&V Team

###  III-B MIL and HIL prototyping for SiS integration purpose

Inspired by the antecipation of issues related to OBDH SiS integration with
the nanosatC-BR2 payloads, the InRob testing process originally proposed in
[15] and its extension in InRob-UML [16] have been used to guide the
construction of communicating models (MIL). These models represent the
expected behaviour when two communicating SiS perform a given service
collaboratively. The models execution allows to both verify interoperability
and robustness requirements by simulating and deriving tests cases using MBT
tools. Figure 4 presents InRob Testing Process in three phases: A) Modeling;
B) Test case Generation; C) Test case Execution.

Figure 4: InRob Testing Process

The inputs for the service modeling (phase A) are: (i) the project artifacts
that specify among others the system functional requirements, the subsystem
user manual and time constraints; (ii) the service profile development, which
supports the prioritization of a particular service to be modeled. The nominal
behavior of the elected service is formally specified. An interoperability
TIOA model results in the nominal TIOA model. From the nominal service model
and time constraints imposed on the communicating SIS, the concept of Major
and / or Minor timing deviations are used to deal with timing hazards related
to the communication channel failures. Thus, the nominal service model is
extended to represent the robust behavior of each communicating SIS, facing
those considered deviations. The output is the extended model of the service.

The extended TIOA models outputs in Phase A are the basis for the test cases
generation methods oriented to test purpose (TP) approaches in Phase B. A TP
is a formal description of a property (or a sequence of properties)
represented in the TIOA models. On-the-fly methods for test case generation
search for a specific property during the simulation and the test case is
generated by recording a trace on traversing the model aiming at the property
verification. This test case is able to validate the behavior of the
implementation concerning that property. Robustness test case generation
considers hazard scenarios (e.g. rushed and duplicate messages) specified by
the test purpose approach.

In Phase C of the testing process, the test cases were performed in real SUT.
The testing architecture relies on a Failure Emulator Mechanism (FEM), which
is considered as part of the Testing System (TS).

The InRob approach has been used in the V&V activities of the nanosatellite
named NanosatC-Br2, a scientific 2U Cubesat-based nanosatellite under
development at INPE. NanosatC-BR2 carries on Earth orbit 5 payloads being 2
scientific instruments and 3 new technologies for validation in the space.
InRob aids to create representative models for interoperability and robustness
analysis of the OBDH expected behavior (requirements) in the interactions with
the nanosatellite payloads. Following InRob Fase A, the modeling was built
with the aid of UPPAAL tool, using timed automata formalism, enabling the
analysis of the requirements through model simulation. Ten models in total
output from Fase A, being two for each payload. Figure 5 shows the nominal
models of the interactions between OBDH and payload SLP (Scientific Langmuir
Proble payload) that occur on the channel I2C in a master-slave communication.
From those models, test cases can be automatically derived through the
verification of properties in the simulation of the models using model-
checking technique, in Fase B.

Figure 5: OBDH and Payload SLP interoperability models

To start the reading, the on-board computer will send a command to the reading
initialization SLP, which will contain bytes referring to the local date and
time (a timestamp to deal with the lack of a Real Time Clock). So, the SLP
will print at the beginning of each reading the date and time of the beginning
of the collection. The SLP must respond to the OBDH the correct receipt of the
command (acknowledge). After starting the reading, the OBDH should estimate a
period of more than 5 minutes to start asking for data transmission (avoiding
to interrupt the SLP in the collection). If an error occurs and the OBDH
requests data transmission to the SLP, the SLP must ignore the request.

The implementation of Fase C in the NanosatC-BR2 count on a Scalable
Architecture Test System (SATS) based on Arduino computer boards proposed in
Figure 6 [14] . SATS supports the V&V process on the MIL and HIL perspectives.
In MIL the codes generated from models presented in Figure 5 are embedded in
the Arduinos and the interactions between OBDH and SLP can be verified from a
behavioral point of view. That facility allows us to anticipate requirement’s
ambiguities and specification faults. Moreover, functional test cases can be
derived from those models. Those test cases will be executed in the HIL
environment when the communicating SiS are embedded in the real hardware, for
acceptance purposes.

The Failure Emulator Mechanism (FEM) in the InRob Test Architecture (see
Figure 6) supports faut injections in the communication channel that is used
for the SiSs (SUT1 and SUT2) interactions in order to validate robustness
requirements. FEM is able to intercept the exchanged messages between two
subsystems under testing and inject different faults: (i) time related faults,
i.e. delay; (ii) value related faults, i.e. bit-flip and (iii) specific
communication bus faults, i.e. verbose subsystem.

Figure 6: InRob Test Architecture instantiated by SATS [14]

Figure 7: Number of Test Cases generated by InRob

Figure 7 presents the number of Test Cases derived from the ten models
generated with the aid of UPPAAL tool in Fase A. For each interaction between
the SUTs models of interoperability, which is considered a Test Purpose (TP),
one Test Case classified as nominal was automatically generated. Regarding the
interoperability models OBDH and SLP, there are 8 interactions specified,
therefore 8 nominal TCs. Moreover, for each nominal TC, at least 3 TCs are
supplemented by FEM during test execution, associated with the three fault
models implemented. In total 32 TCs were executed. For the five payloads of
NanosatC-Br2, from 10 interoperability models specified using InRob, 232 test
cases were automatically created, 58 nominal functional and 174 fault
injection test cases (applied by FEM).

The conception of a HIL environment starts in the EM integration as a
continuous integration until the CDR. Each subsystem, in special the payloads,
needs to be validated. In terms of interoperability, behaviour and
requirements conformity, all the payloads engineering models must present a
certain level of acceptance from the OBDH point of view. The HIL environment
is a great opportunity to isolate the possible sources of failures, once each
payload and even the OBDH can be tested on 1-by-1 operation. The Test Reports
from the MIL environment (SATS) should guide this phase and point out the non-
conformities and misunderstandings of the development team (OBDH and
payloads). Due to some lack of expertise of the NanoSatC-BR2 development team,
the HIL approach was not that much exercised and the project was rushed to be
ready for complete integration. We will discuss the EM and FM integration at
the Section IV.

##  IV Learning Lessons and Reduce Complexity of Verification Work

This section discusses the effect of the V&V processes adopted by INPE for
CubeSat standard nanosatellites projects, highlighting the changes introduced
by the use of Model-Based approaches. The process is supported by ECSS
standards, but it is affected by COTS at very low procurement cost and a
simple standardized architecture. The way a space organization deals with
downsizing its traditional V&V process to carry on V&V of nanosatellites is
not obvious.

It is observed in the SiS V&V process presented in Section III that the
process applied to nanosatellites is incremental due to both the
characteristics of the requirements that are more volatile and the greater
importance of interoperability in the integration of the nanosatellites
subsystems. Due to the volatility of the requirements and the higher
acceptable risk in the nanosatellites project, the closing of the review
documents is quite delayed, reaching cases in which the release is made a
couple of phases later in the development process (see Figure I, the release
of the Requirements Baseline and Interface Requirement Document) or even not
being formally produced (in Software Development Plan and Technical
Specification). However, problems are identified when we start to deal with
the spacecraft integration, EM and FM.

Both integration phases will point to all kind of misunderstandings, non-
conformities and errors, even mechanical, electrical or logical, or all of
them, at the same time. A great part of the documents are released during this
phases as result of gaps found during the integration. Thus, many requirements
and constraints are added to the project to fit the solutions – reworking
hardware/software development is more difficult than reworking a document –
instead of the solutions to fit the requirements and constraints. This ends up
with lack of information on the documents, lack of traceability and
vulnerabilities associated with conflicting implementations.

Regarding the SiS discipline, the reduced complexity of verification of
nanosatellites has been obtained by INPE with the lessons learned to support
an automatic test generation adopted in the V&V process of traditional
satellites. Although the modelling activities are very time consuming, the
benefits of the effective test suite provided allows INPE to increase
confidence in model-based approaches. However, instead of effort in acceptance
testing, which does not make sense in the nanosatellite projects, INPE focused
on the use of modeling approaches to anticipate the detection and treatment of
embedded software-intensive communicating subsystems in the development life
cycle of SiS. The goal is to avoid software rework due to interoperability
faults usually evidenced only in the integration testing.

The EM integration is aimed at the verification of every software subsystem
interoperability, beyond the mechanical and electrical integration. Being the
OBDH the core of the satellite coordination, all the EM integration tests
focus on how the OBDH software deals with other subsystems, centralizing their
control and monitoring. The Test Cases used at MIL and HIL environments can be
used at this point as a guide for the expected behaviour of the OBDH and its
payloads. In this case, the models can be used as the single source of truth
(SSOT). However, to ensure the closing of the EM integration phase, the space-
to-ground integration must be addressed.

Not only the spacecraft must be working properly on the inside, it must be
controllable from the ground. The model-based approach can be useful here, as
the concept of spacecraft operation can be modelled in a higher level than the
way it was done before with the subsystems. However, between the space and
ground segments there is a communication channel that cannot be modelled in
the same way preventing us to directly model the validation of
interoperability between them. So, we use both the EM model and the Satellite
Control Software directly to verify the complete operation of the spacecraft
and exercise possible flight plans for the mission.

Ending the CDR, the FM integration can start. This phase does not differ much
from the EM integration as we are trying to replicate the expected behaviour
found at the EM. The main differences occur in the replacement of the
subsystem engineering models by their flight models. Some of them are more
fragile (e.g. solar panels instead of dummy Printed Circuit Boards - PCBs) or
more robust (e.g. aluminium structure instead of a mock-up). The possible
problems that can be found here are more related to the assembly and cabling,
items that it is not possible to model or measure in the initial phases. At
this stage, the chance of finding logical errors in the spacecraft software
and its operation is lower.

For the NanoSatC-BR2 mission, all these lessons were pointed out. They allowed
us to learn how to reduce the complexity of the V&V process, in qualitative
and quantitative ways. The gaps in the team expertise is difficult to by-pass
as, most of time, the missions are university-class projects.

Therefore, INPE has made progress on research that contributes with a model-
based approach in terms of (i) alternative modeling formalisms [17]; (ii) new
tools to aim deriving test cases from models; (iii) open architectures to
support the test cases execution. A Scalable Architecture Test System (SATS)
based on Arduino computer boards is proposed in [14] . It allows to
systematize test scenarios based on the concept of model and hardware in the
loop. The approach combines Model-Based Testing (MBT) and Model-Driven
Engineering (MDE) techniques to support: (i) specification of interoperability
requirements, whose behavior is modeled in state machines; (ii) automated
software code generation from those models to be embedded in Arduino computer
boards; (iii) automated generation of nominal test cases focusing on
interoperability properties; (iv) extension of models with robustness aspects;
and finally (v) automated generation of test cases focusing on robustness
properties. Moreover, a failure emulator mechanism framework (FEM) is proposed
in [18] for robustness testing of interoperable software-intensive subsystems
on-board nanosatellite. FEM acts in the communication channel being part of
the integration test workbench in two phases of nanosatellite design: (i)
robustness requirement specification using model in the loop (MIL) and (ii)
robustness validation using hardware in the loop (HIL). Both approaches have
been applied in the V&V process of NanosatC-BR2 [17].

##  V Conclusions and future work

This work presented the V&V process, used at INPE, applied in nanosatellites
projects that follow the CubeSat standard. The structure or even the
subsystems of these nanosatellites are fully acquired as COTS. Therefore, the
V&V process usually adopted in traditional satellites, as recommended by ECSS
standards, needs to be adapted to reflect the reduced complexity of work and
the absence of one or more disciplines.

We presented the rationale to conceive a V&V process for software-intensive
subsystem on-board nanosatellites developed by INPE aiming to downsize the
practices used in traditional satellite V&V process. The rationale highlights
the importance of the interoperability among the OBDH SiS and the satellite
payloads. Being key elements, the success of the mission in orbit depends on
confidence in their operation.

Research opportunities were identified to improve the use of MBT approach,
firstly applied in the process of the traditional satellites, and then to
anticipate integration testing in the development life cycle of the Cubesat-
based space missions using InRob. The development of standardized test
environments using model-based methods and tools allowed the establishment of
integration testing processes for satellite subsystems that are more effective
in time and cost. The complexity of the subsystems, increasing in
functionalities implemented by software (SiS) is supported by the
systematization of the generation and execution of interoperability and
robustness tests provided by inRob, which has shown a strong alignment with
the new philosophy of small satellite projects. So, standardisation,
systematisation and automation of the V&V process may help to work around
problems (unknown requirements in the begin of the development cycle and their
constant evolution, immaturity of the stakeholders and the lack of expertise
of the development team regarding the specificities of software engineering
for space area and low budget forcing the allocation of inexperienced
developers in the SiS context) and reach a more compliant results.
Furthermore, we believe that this V&V process can have a broader impact to
several embedded systems and infrastructures for critical machines.

In the near future, through scientific partnerships with national and
international academic institutions, INPE hopes to consolidate the V&V
process, based on a test system that supports the approach of evolutionary
integration tests for the development of SiS embedded in nanosatellites. The
goal is to use the InRob method and associated tools at scale (model
construction, automatic test generation and automated instrumentation in the
execution of tests), systematizing the verification of the SiS shipped in the
different nanosatellites that will constitute the new generation of The
Brazilian Environmental Data Collecting Systems (BEDCS), in accordance with
the CONASAT222http://www.crn.inpe.br/conasat1/nanosatt.php project.

##  VI Acknowledgment

This work has been partially supported by the proj. ADVANCE - Addressing V&V
Challenges in Future Cyber-Physical Syst. (https://www.advance-rise.eu/
H2020-MSCA-RISE-2018, n. 823788) and partially supported by the Centro de
Inform. e Sistemas da Univers. de Coimbra (CISUC).

## References

  * [1] R. Twiggs, “Origin of cubesat,” _Small Satellites: Past, Present, Future, Eds: Helvajian H., Janson SW, The Aerosp Press, California_ , 2008.
  * [2] E. C. for Space Standardization, “SECSS-E-ST-10C: Space eng. - System eng. general requirements,” ECSS, Active Standards, 2017.
  * [3] NASA, _NASA System Engineering Handbook Revision 2_. NASA Headquarters, 2016.
  * [4] A. Scholz, “Cubesat standards handbook,” _The LibreCube Initiative_ , 2017\.
  * [5] E. C. for Space Standard., “ECSS-E-40 Part 1B: Space eng. – Software – Part 1: Principles and requirements,” ECSS, Active Standards, 2003.
  * [6] J. Kiesbye, D. Messmann, M. Preisinger, G. Reina, D. Nagy, F. Schummer, M. Mostad, T. Kale, and M. Langer, “Hardware-in-the-loop and software-in-the-loop testing of the move-ii cubesat,” _Aerospace_ , vol. 6, no. 12, p. 130, 2019.
  * [7] H. Leppinen, P. Niemelä, N. Silva, H. Sanmark, H. Forstén, A. Yanes, R. Modrzewski, A. Kestilä, and J. Praks, “Developing a linux-based nanosatellite on-board computer: flight results from the aalto-1 mission,” _IEEE Aerospace and Elect. Syst. Magaz._ , vol. 34, no. 1, pp. 4–14, 2019\.
  * [8] A. Alanazi and J. Straub, “Engineering methodology for student-driven cubesats,” _Aerospace_ , vol. 6, no. 5, p. 54, 2019.
  * [9] C. Venturini, B. Braun, D. Hinkley, and G. Berg, “Improving mission success of cubesats,” 2018.
  * [10] J. Eickhoff, R. Hendricks, and J. Flemmig, “Model based development and verification environment,” in _54th Int. Astronautical Cong. of the Int. Astronautical Fed. (IAF)_ , vol. 3. AIAA, dec 2003, pp. 1275–1285.
  * [11] M. Omair Khan, M. Sievers, and S. Standley, “Model-based V&V of spacecraft avionics,” _AIAA Infotech at Aerospace Conf. and Exhibit 2012_ , pp. 1–11, 2012.
  * [12] B. Tiseo, V. Quaranta, G. Bruno, and G. Sisinni, “Tailoring of ECSS Standard for Space Qualification Test of CubeSat Nano-Satellite,” vol. 13, no. 4, pp. 295–302, 2019.
  * [13] R. Feldt, R. Torkar, E. Ahmad, and B. Raza, “Challenges with software V&V activities in the space industry,” _ICST 2010 - 3rd Int Conf on Sw Test, V &V_, pp. 225–234, 2010.
  * [14] C. Conceição, “Abordagem sistemática de testes de software comunicante embarcado em nanosatélites com foco em falhas de interoperabilidade,” Ph.D. dissertation, CSE/INPE, http://urlib.net/8JMKD3MGP3W34R/3TBGC2H, 5 2019, in Portuguese.
  * [15] F. Mattiello-Francisco, E. Martins, A. R. Cavalli, and E. T. Yano, “Inrob: An approach for testing interoperability and robustness of real-time embedded software,” _Journ. of Syst. and Soft._ , vol. 85, no. 1, pp. 3–15, 2012\.
  * [16] A. C. Weller, E. Martins, and F. Mattiello-Francisco, “Inrob-uml: uma abordagem para testes de interoperabilidade e robustez baseados em modelos,” _SAST 2015_ , p. 71, 2015.
  * [17] D. Almeida and F. Mattiello-Francisco, “Modeling of the interoperability between on-board computer and payloads of the nanosat-br2 with support of the uppaal tool,” in _1st IAA Latin American Symp. on Small Satellites. Colombia, Session_ , vol. 9, 2017.
  * [18] C. Batista, A. Weller, E. Martins, and M. Mattiello-Francisco, “Towards increasing nanosatellite subsystem robustness,” _Acta Astronautica_ , vol. 156, pp. 187–196, 2019.