**About the Editor**

#### **Patricia J. Y. Wong**

Dr. Patricia J. Y. Wong joined the School of Electrical and Electronic Engineering in 2000. She received all her degrees, BSc (first class honours), MSc, MSc (Financial Engineering) and PhD, from the National University of Singapore. Her research interests include differential equations, difference equations, integral equations, and numerical mathematics. She has conducted significant research and published 4 monographs and more than 200 research papers in international refereed journals. Presently, Dr Wong serves on the editorial boards of *Mathematics* (Switzerland), *Axioms* (Switzerland), *Abstract and Applied Analysis* (USA), *Advances in Dynamical Systems and Applications* (USA), and the *International Journal of Differential Equations* (USA).

## *Article* **Optimal Control Theory for a System of Partial Differential Equations Associated with Stratified Fluids**

**Edgardo Alvarez 1,\*, Hernan Cabrales 2,† and Tovias Castro 3,†**


**Abstract:** In this paper, we investigate the existence of an optimal solution of a functional restricted to non-linear partial differential equations, which ruled the dynamics of viscous and incompressible stratified fluids in R3. Additionally, we use the first derivative of the considered functional to establish the necessary condition of the optimality for the optimal solution.

**Keywords:** non-linear optimal control; stratified fluids; energy functional; optimal condition; state variable

#### **1. Introduction**

Following the results of the modern calculus of variations, in this article, we study the optimal solution of an energy functional constraint to a partial differential equations system, which models the dynamic of an exponential stratified fluid in a three-dimensional space. To do this, we investigate the existence of solutions of a non-homogenous and non-linear partial differential system, extending the result obtained in [1], where only a potential external force was considered. Being more specific, for a <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup>3, non-empty, open, connected, and bounded set, with boundary Σ = *∂*Ω × (0, *T*) that is smooth enough (at least Lipschitz continuous) and letting *ν* be the normal vector outside the boundary, we define *Q* := Ω × (0, *T*) as the domain of our model where the motion of the fluid takes place. Here, *<sup>T</sup>* > 0, (0, *<sup>T</sup>*) is the time interval and *<sup>t</sup>* <sup>∈</sup> (0; *<sup>T</sup>*) is the temporal variable.

We are interested in establishing the existence of the solution for the following nonlinear problem in a weaker sense:

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

Received: 10 September 2021 Accepted: 11 October 2021 Published: 21 October 2021

**Citation:** Alvarez, E.; Cabrales, H.; Castro, T. Optimal Control Theory for a System of Partial Differential Equations Associated with Stratified Fluids. *Mathematics* **2021**, *9*, 2672. https://doi.org/10.3390/math9212672 Academic Editor: Patricia J. Y. Wong

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

$$\begin{cases} \begin{aligned} \frac{\partial y\_1}{\partial t} - \mu \Delta y\_1 + y \cdot \nabla y\_1 + \frac{\partial p}{\partial x\_1} &= u\_1, \\ \frac{\partial y\_2}{\partial t} - \mu \Delta y\_2 + y \cdot \nabla y\_2 + \frac{\partial p}{\partial x\_2} &= u\_2, \end{aligned} \\ \begin{aligned} \frac{\partial y\_3}{\partial t} - \mu \Delta y\_3 + g\rho + y \cdot \nabla y\_3 + \frac{\partial p}{\partial x\_3} &= u\_3, \\ \frac{\partial \rho}{\partial t} - \frac{N^2}{g} y\_3 &= u\_4, \end{aligned} \end{cases} \tag{1}$$
 
$$\begin{aligned} \frac{\partial y\_1}{\partial x\_1} + \frac{\partial y\_2}{\partial x\_2} + \frac{\partial y\_3}{\partial x\_3} &= 0, \end{aligned} \tag{2}$$

where *x* = (*x*1, *x*2, *x*3) denotes the spatial variable, *y* = *y*(*x*, *t*)=(*y*1(*x*, *t*), *y*2(*x*, *t*), *y*3(*x*, *t*)) denotes the velocity field of the fluid and *u*(*x*, *t*)=(*u*1, *u*2, *u*3, *u*4) corresponds to a known function from *<sup>L</sup>*2(Ω). We also have the parameter *<sup>μ</sup>* > 0 as the kinematic viscosity, and *<sup>N</sup>*

$$\begin{aligned} \frac{\partial y\_1}{\partial t} - \mu \Delta y\_1 + \boldsymbol{y} \cdot \nabla y\_1 + \frac{\partial p}{\partial \boldsymbol{x}\_1} &= u\_1, \\ \frac{\partial y\_2}{\partial t} - \mu \Delta y\_2 + \boldsymbol{y} \cdot \nabla y\_2 + \frac{\partial p}{\partial \boldsymbol{x}\_2} &= u\_2, \\ -\mu \Delta y\_3 + \xi \rho + \boldsymbol{y} \cdot \nabla y\_3 + \frac{\partial p}{\partial \boldsymbol{x}\_3} &= u\_3, \\ \frac{\partial \rho}{\partial t} - \frac{N^2}{\mathcal{g}} y\_3 &= u\_4, \\ \frac{\partial y\_1}{\partial \boldsymbol{x}\_1} + \frac{\partial y\_2}{\partial \boldsymbol{x}\_2} + \frac{\partial y\_3}{\partial \boldsymbol{x}\_3} &= 0, \end{aligned}$$

and *g* are positive constants. The last equation for the non-linear system is because our fluid is incompressible, *p* denotes the scalar field of the dynamic pressure and *ρ* represents the dynamic density. For an ideal case, the Equation (1) can be founded in [2,3]. For a viscous compressible fluid, the system (1) is deduced, for example, in [4].

When we study optimal control problems, we start from a dynamical system that evolves temporarily in a period time [*t*0, *tf* ], described by a state equation of a specific variable *y*(*t*), called a state variable with an initial condition *y*0. This evolution of the system depends on a particular function *u*(*t*), called a control variable, and what is sought with it is to influence the evolution of *y*(*t*) such that we can optimize (maximize or minimize) a given functional, which depends on both the state and control variables, called the *energy functional*. To be more related to these terms of the theory of optimal control, we can see [5,6].

The primary motivation of this paper is to minimize an energy functional of the form *J*(*y*, *u*), which depends on a control variable *u* and the velocity field *y* subject to a state equation that corresponds in our case, to a non-linear system of partial differential equations given in (1).

The functional that we are going to minimize is defined by:

$$J = \frac{1}{2} \int\_0^T \left\|(y\_\prime \rho)^T - y\_d\right\|\_{L^2(\Omega)}^2 \, \mathrm{d}t + \frac{\lambda}{2} \int\_0^T \left\|u - u\_d\right\|\_{L^2(\Omega)}^2 \, \mathrm{d}t\_\prime$$

where *yd* <sup>∈</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> is the desired state, *ud* <sup>∈</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> is the desired control (or also called control change) and *<sup>λ</sup>* > 0 is a constant. From a mathematical perspective, most of the control systems involve a set of ordinary differential equations or linear partial differential equations in their restrictions, see for example [7].

In this case, we consider a non-linear model, which makes this proposal novel and attractive. On the other hand, there is some progress associated with the Navier–Stokes systems [5,8]. However, not much seems to be known about works that deal with non-linear exponential stratified fluids, making our results an open door to consider new parameters such as salinity, rotation, and temperature in future works.

This paper is distributed in five sections. In Section 1, we introduce and describe the problem; later, in Section 2, we show the essential background information to understand the problem. In Section 3, we introduce the weak formulation of the problem. In Section 4, we study the existence of solutions for the optimal problem, and finally, in Section 5, we establish the optimal condition.

#### **2. Previous Definitions and Notations**

Before starting the study and analysis of our optimal control problem, we introduce some previous elements and necessary notation to understand the non-linear motion in the dynamics of viscous and incompressible stratified fluids in R<sup>3</sup> that will be considered this paper.

Let <sup>Ω</sup> be a domain of the space <sup>R</sup>3, and let *<sup>p</sup>* in <sup>R</sup>, such that 1 <sup>≤</sup> *<sup>p</sup>* <sup>≤</sup> <sup>∞</sup>. A function *<sup>y</sup>* : <sup>Ω</sup> −→ <sup>R</sup> (or <sup>C</sup>), is said to belong to *<sup>L</sup>p*(Ω), if *<sup>y</sup>* is measurable and the norm

$$||y||\_{L^{p}(\Omega)} = \begin{cases} \left(\int\_{\Omega} |y(\mathfrak{x})|^{p} \, \mathrm{d}\mathfrak{x}\right)^{1/p} & \text{if } 1 \le \infty \\\\ \mathop{\mathrm{ess}\, \mathrm{sup}}\, |y(\mathfrak{x})| & \text{if } p = \infty, \end{cases}$$

is finite. The spaces *Lp*(Ω) are Banach spaces (see [9]). Furthermore, in the spaces *Lp*(Ω) the Hölder Inequality is fulfilled, which ensures that, for *<sup>y</sup>* <sup>∈</sup> *<sup>L</sup>p*(Ω) and *<sup>v</sup>* <sup>∈</sup> *<sup>L</sup>q*(Ω) with 1 *p* + 1 *<sup>q</sup>* <sup>=</sup> 1 for 1 <sup>≤</sup> *<sup>p</sup>*, *<sup>q</sup>* <sup>≤</sup> <sup>∞</sup>, it holds:

$$\int\_{\Omega} |y(\mathfrak{x})v(\mathfrak{x})| \mathrm{d}\mathfrak{x} \le \|y(\mathfrak{x})\|\_{L^{p}(\Omega)} \cdot \|v(\mathfrak{x})\|\_{L^{q}(\Omega)}.$$

In particular, when we have that *p* = 2, then *L*2(Ω) is a Hilbert space with the scalar product

$$(y,v)\_2 = \int\_{\Omega} y(\mathfrak{x}) \cdot v(\mathfrak{x}) \, \mathrm{d}x.$$

It is known that *L*2(Ω) is one of the essential Hilbert spaces in the mathematical analysis since they appear very frequently in the study of partial differential equations, and it is the space where the kinetic energy is automatically well defined. As the variational form of a mathematical physics problem appears, we cross the Sobolev's spaces denoted by *<sup>W</sup>k*,*p*(Ω), and defined as the set of all functions *<sup>y</sup>*(*x*) <sup>∈</sup> *<sup>L</sup>p*(Ω) that have all the generalized derivatives up to the order *p*, which also belongs to *Lp*(Ω). The associated norm defined in this space is given by

$$||y||\_{W^{k,p}(\Omega)} = \left(\sum\_{|\alpha| \le k} ||D^{\alpha}y||\_{L^p(\Omega)}^p\right)^{1/p}$$

where *Dαy* is the weak derivate of order *α*. We also find other types of Sobolev spaces such as *Wk*,*<sup>p</sup>* <sup>0</sup> (Ω).

Note that when *p* = 2, we can simply write *Hk*(Ω) and *H<sup>k</sup>* <sup>0</sup>(Ω) instead of *<sup>W</sup>k*,2(Ω) and *Wk*,2 <sup>0</sup> (Ω), respectively (see for example [10]). Furthermore remember that when *k* = 1 and *p* = 2, we have that the space *W*1,2(Ω) is better known as *H*1(Ω), since it is a Hilbert space, endowed with the scalar product:

$$\mathbb{P}\left(y, v\right)\_{H^1(\Omega)} = \int\_{\Omega} y(\mathbf{x}) \cdot v(\mathbf{x}) \, d\mathbf{x} + \int\_{\Omega} \nabla y(\mathbf{x}) \cdot \nabla v(\mathbf{x}) \, d\mathbf{x} \text{ for all } y, v \in H^1(\Omega),$$

where

$$
\nabla y = \left(\frac{\partial y}{\partial x\_1}, \frac{\partial y}{\partial x\_2}, \frac{\partial y}{\partial x\_3}\right),
$$

and

$$\nabla v = \left(\frac{\partial v}{\partial \mathbf{x}\_1}, \frac{\partial v}{\partial \mathbf{x}\_2}, \frac{\partial v}{\partial \mathbf{x}\_3}\right).$$

The norm induced by the previous scalar product is given by:

$$\|y\|\_{H^1(\Omega)} = \left(\|y\|\_{L^2(\Omega)}^2 + \sum\_{i=k}^3 \left\|\frac{\partial y}{\partial \mathfrak{x}\_i}\right\|\_{L^2(\Omega)}^2\right)^{1/2}$$

.

On the other hand, let us denote by <sup>D</sup>(Ω) the space of functions *<sup>ϕ</sup>* : <sup>Ω</sup> −→ <sup>R</sup><sup>3</sup> of class *<sup>C</sup>*∞(Ω) with compact support and by <sup>D</sup> (Ω) the space of distributions on Ω. Throughout this paper, we will use the standard notations for the Lebesgue and Sobolev spaces, in particular the norm in *L*2(Ω) and the scalar product in *L*2(Ω) will be represented by · and (·, ·) respectively.

Let us define

$$(u, v) := \int\_{\Omega} \sum\_{j=1}^{3} u\_j \cdot v\_j \, \mathrm{d}x, \; u = (u\_1, u\_2, u\_3), \; v = (v\_1, v\_2, v\_3) \in L^2(\Omega)^3,$$

$$((u, v)) := \int\_{\Omega} \sum\_{j=1}^{3} \nabla u\_j \cdot \nabla v\_j \, \mathrm{d}x, \; u = (u\_1, u\_2, u\_3), \; v = (v\_1, v\_2, v\_3) \in H\_0^1(\Omega)^3,$$

and the associated norms are given to from | *u* | 2:= (*u*, *<sup>u</sup>*) and *u*<sup>2</sup> := ((*u*, *<sup>u</sup>*)).

Consider the following notation for the solenoidal Banach spaces *H* and *V*, which intrinsically satisfy the condition ∇ · *y* = 0, and which we can represent as:

$$H = \{ y \in L^2(\Omega)^3 : \nabla \cdot y = 0 \text{ in } \Omega; \ \gamma\_\mathbb{H} \ y = 0 \text{ on } \partial\Omega \},$$

and

$$V = \{ y \in H\_0^1(\Omega)^3 : \nabla \cdot y = 0 \text{ in } \Omega \}. \tag{2}$$

Here, ∇ · *y* denotes the divergence of *y* and *γ<sup>n</sup>* denotes the normal component of the trace operator, where *γ<sup>n</sup>* : *y* → *n* · *y ∂*Ω = 0, here *n* denotes the external normal to the boundary. These spaces are used very frequently in equations of the dynamics of the stratified fluids and are defined as the closure of Θ in *L*2(Ω)<sup>3</sup> and of Θ in *H*<sup>1</sup> <sup>0</sup> (Ω)3, respectively, where

$$\Theta = \{ y \in \mathcal{D}(\Omega)^3 : \nabla \cdot y = 0 \text{ in } \Omega \}.$$

It is well-known that *H* and *V* are Hilbert spaces with the scalar product (·, ·) and · respectively. Furthermore,

$$V \subset H \equiv H' \subset V'\_{\prime \prime}$$

where injections are dense and continuous.

On the other hand, if *V* is a Banach space with dual space *V* , then the duality between the spaces *V* and *V* is denoted by ·, ·*V*,*<sup>V</sup>* and its associated norm in *V* is denoted by ·*V* .

We introduce the following space of functions *y* whose derivative *yt* exists as an abstract function:

$$\mathcal{W}^{\mathfrak{a}}([0,T];X) := \{ y \in L^2([0,T];X) : y\_t \in L^{\mathfrak{a}}([0,T];X') \}, \quad 1 \le \mathfrak{a} \le 2\pi$$

$$\mathcal{W}(0,T) := \mathcal{W}^2([0,T];X).$$

The spaces defined above are endowed with the following norms:

$$\begin{aligned} \|\|y\|\|\_{W^a([0,T];\mathcal{X})} &:= \left(\|\|y\|\|\_{L^2([0,T];\mathcal{X})}^2 + \|\|y\_t\|\|\_{L^a([0,T];\mathcal{X}')}^2\right)^{1/2} \\\\ \|\|y\|\|\_{\mathcal{W}} &:= \|\|y\|\|\_{\mathcal{W}^{2,\epsilon}} \end{aligned}$$

are Banach spaces. When *X* is a Hilbert space, we have that *Wα*([0, *T*]; *X*) is a Hilbert space. In particular, the space *Wα*([0, *T*]; *X*) is endowed with the following scalar product:

$$(y, v)\_{W^u([0, T]; \mathcal{X})} = \int\_0^T (y(t), v(t)) \, \mathrm{d}t + \int\_0^T (y\_t(t), v\_t(t)) \, \mathrm{d}t.$$

In this way, we have the following results for 1 ≤ *α* ≤ 2 (see [10–12]):

$$\begin{aligned} \mathcal{W}^{\mathfrak{a}}([0,T];X) &\hookrightarrow \mathcal{C}([0,T];X) \text{ is continuous},\\ \mathcal{W}^{\mathfrak{a}}([0,T];H^{1}(\Omega)^{3}) &\hookrightarrow L^{2}(\Omega)^{3} \text{ is compact},\\ \mathcal{W}^{\mathfrak{a}}([0,T];H^{1}(\Omega)^{3}) &\hookrightarrow \mathcal{C}([0,T];L^{2}(\Omega)^{3}) \text{ is compact}.\end{aligned}$$

Now, we defined our set of admissible controls, which denoted by *Uad*, and its elements are called admissible controls, which satisfy the inequality constraints of our non-linear system given from:

$$\mathcal{U}\_{\text{all}} = \{ u \in L^2(Q\_T)^3 : u\_{\text{if}}(\mathbf{x}, t) \le u\_i(\mathbf{x}, t) \le u\_{\text{if}}(\mathbf{x}, t) \text{ with a.e. on } Q\_T \text{ for } i = 1, 2, 3, 4 \}. \tag{3}$$

where the control constraints *ua*, *ub* <sup>∈</sup> *<sup>L</sup>*2(*QT*) with

$$
u\_{a,i}(\mathbf{x},t) \le \mathcal{u}\_{b,i}(\mathbf{x},t).$$

**Remark 1.** *Note that our set of admissible controls defined in* (3) *is a non-empty, convex, and closed subset in L*2(*QT*)4*.*

Now, let us recall the following classic result that we will need later to show the existence of optimal controls.

**Definition 1** ([6])**.** *Let <sup>X</sup> be a Banach space and let <sup>J</sup>* : *<sup>X</sup>* −→ <sup>R</sup> *be a functional. We say that <sup>J</sup> is weakly lower semicontinuous, if for any sequence* (*xn*)*n*∈<sup>N</sup> <sup>⊂</sup> *<sup>X</sup> such that xn <sup>x</sup> when <sup>n</sup>* −→ <sup>∞</sup> *we have that:*

$$J(\mathfrak{x}) \le \liminf\_{n \to \infty} J(\mathfrak{x}\_n).$$

#### *Formulation of the Optimal Control Problem Associated with the Non-Linear Model*

In this part, we will formulate our optimal control problem associated with the partial differential equation described by (1). In order to show the existence of solutions of the non-linear system (1), we represent our model system in a simpler way using the following notation:

*y*(*x*, *t*) = ⎛ ⎜⎜⎝ *y*1 *y*2 *y*3 *ρ* ⎞ ⎟⎟⎠ = *y ρ* ; *∂ ∂t* ⎛ ⎜⎜⎝ *y*1 *y*2 *y*3 *ρ* ⎞ ⎟⎟⎠ <sup>=</sup> *<sup>∂</sup><sup>y</sup> <sup>∂</sup><sup>t</sup>* ; (*y* · ∇)*<sup>y</sup>* <sup>=</sup> ⎛ ⎜⎜⎝ *y* · ∇*y*<sup>1</sup> *y* · ∇*y*<sup>2</sup> *y* · ∇*y*<sup>3</sup> 0 ⎞ ⎟⎟⎠ ; *μ*Δ*y* = ⎛ ⎜⎜⎝ *μ*Δ*y*<sup>1</sup> *μ*Δ*y*<sup>2</sup> *μ*Δ*y*<sup>3</sup> 0 ⎞ ⎟⎟⎠,

and

$$M\mathfrak{y} = \begin{pmatrix} 0\\ 0\\ \mathcal{g}\rho\\ -\frac{N^2}{\mathcal{S}}\mathfrak{y}\_3 \end{pmatrix}; \quad \nabla p = \begin{pmatrix} \frac{\partial p}{\partial \mathcal{X}\_1}\\ \frac{\partial p}{\partial \mathcal{X}\_2}\\ \frac{\partial p}{\partial \mathcal{X}\_3}\\ \frac{\partial \mathcal{X}\_3}{\partial \mathcal{X}\_4} \end{pmatrix}; \quad \mathfrak{u} = \begin{pmatrix} u\_1\\ u\_2\\ u\_3\\ u\_4 \end{pmatrix}.$$

Then, we can rewrite (1) in more compact form as,

$$\begin{cases} \frac{\partial y}{\partial t} + (y' \cdot \nabla)y - \mu \Delta y + My + \nabla p = u, \quad y \in \Omega \times \mathbb{R}, \\ \text{div } (y') = 0, \ x \in \Omega \text{ and } t \ge 0, \\ y(t, \cdot) = 0 \text{ on the boundary of } \Sigma = \partial \Omega \times (0, T), \\ y(0, \cdot) = y\_0 \text{ in } \Omega. \end{cases} \tag{4}$$

In this way, we can introduce our energy functional that we want to minimize, which depends on the state and the control (*y*, *u*) and that we define by:

$$J(y, u) = \frac{1}{2} \int\_0^T \left\| y - y\_d \right\|\_{L^2(\Omega)}^2 \, \mathrm{d}t + \frac{\lambda}{2} \int\_0^T \left\| u - u\_d \right\|\_{L^2(\Omega)}^2 \, \mathrm{d}t, \,\,\, \mathrm{d}t$$

where *yd* <sup>∈</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> is the desired state, *ud* <sup>∈</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> is the desired control (or also called control change) and *<sup>λ</sup>* > 0 is a constant.

Now, we can introduce the functional space given by:

$$\mathbb{V} = \{ \boldsymbol{\varrho}(\mathbf{x}) = (\boldsymbol{\varrho}\_1, \boldsymbol{\varrho}\_2, \boldsymbol{\varrho}\_3, \boldsymbol{\varrho}\_4) : \boldsymbol{\varrho} \in H\_0^1(\Omega)^4 : \nabla \boldsymbol{\varrho}' = 0 \}, \tag{5}$$

where ∇ · *<sup>ϕ</sup>* <sup>=</sup> <sup>0</sup> denotes the divergence of *<sup>ϕ</sup>* = (*ϕ*1, *<sup>ϕ</sup>*2, *<sup>ϕ</sup>*3). The space <sup>V</sup> endowed with the inner product and the usual norm of space *H*<sup>1</sup> <sup>0</sup> (Ω)4, the space of all functions *<sup>y</sup>* <sup>∈</sup> *<sup>H</sup>*1(Ω) with null trace:

$$(y\_\prime v)\_{H^1\_0(\Omega)} := \sum\_{|\alpha|=1} (D^\alpha\_x y\_\prime D^\alpha\_x v)\_{2\alpha}$$

and

$$\|\|y\|\|\_{H^1\_0(\Omega)} = \sqrt{(y,y)\_{H^1\_0(\Omega)}} \text{ for all } u,v \in \mathbb{V}.$$

The space given by (5) will be of great importance to us, since through it we can find the functions *<sup>y</sup>* : [0, *<sup>T</sup>*] −→ <sup>V</sup>, which are weak solutions of our non-linear problem given by (4). On the other hand, our space <sup>V</sup> is clearly a Banach space with the norm ·*H*<sup>1</sup> <sup>0</sup> (Ω). In

this way, it is a Hilbert space. It is also reflexive since <sup>V</sup> <sup>⊆</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) is separable.

In summary, we can establish our optimal control problem:

$$\int \text{Minimize } J(y, u) = \frac{1}{2} \int\_0^T \|y - y\_d\|\_{L^2(\Omega)}^2 \, \text{d}t + \frac{\lambda}{2} \int\_0^T \|u - u\_d\|\_{L^2(\Omega)}^2 \, \text{d}t,\tag{6}$$

subject to the state equations that establish the dependency between the state variable *y* and the control variable *u*:

$$\begin{cases} \frac{\partial y}{\partial t} + (y' \cdot \nabla)y - \mu \Delta y + My + \nabla p = u \text{ in } Q\_{T'}\\ \text{div } (y') = 0, \text{ in } Q\_{\prime} \\ y(t\_{\prime} \cdot) = 0 \text{ on the boundary of } \Sigma = \partial \Omega \times (0, T), \\ y(0\_{\prime} \cdot) = y\_0 \text{ in } \Omega, \\ u \in U\_{ad}. \end{cases} \tag{7}$$

Here, *<sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2(*QT*) is the control, it is an external force that affects the fluid, (for example gravity); *<sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup> is a divergence-free vector field in <sup>R</sup>3, the kinematic viscosity *<sup>μ</sup>* <sup>&</sup>gt; 0 and *Uad* represents our set of constraints defined as in (3).

The aim of our control problem is to find a solution *<sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2(*QT*), where *<sup>y</sup>* is the solution of (7) associated with *u* such that it minimizes our energy functional given by (6).

In this paper, we will show the existence of solutions and establish the use of the first derivative of the energy functional to derive the conditions that the optimal solutions have to satisfy Equations (6) and (7).

#### **3. Weak Formulation for the Non-Linear Problem**

In this section, whenever we refer to space V, we will work with the functional space defined by (5), we will also identify <sup>V</sup><sup>∗</sup> as its dual space and (·, ·) and · will denote the scalar product and the usual norm in *L*2(Ω), respectively.

We are interested in establishing theorems of existence and uniqueness of the solution for our non-linear problem given by (4), for which we first study the existence of a solution in a weaker sense.

First of all, suppose there are functions *<sup>y</sup>* <sup>∈</sup> *<sup>C</sup>*2,1(<sup>Ω</sup> <sup>×</sup> (0, *<sup>T</sup>*)) and <sup>∇</sup>*<sup>p</sup>* <sup>∈</sup> *<sup>C</sup>*(*QT*) classical solutions for our non-linear system of partial differential equations given by (4). Let us show the *weak formulation* for our non-linear problem:

Suppose that *<sup>y</sup>* is a solution of (4). Then, multiplying the first equation of (4) by *<sup>v</sup>* <sup>∈</sup> <sup>V</sup>, we obtain the following:

$$(y\_t)\cdot \upsilon + (y'\cdot \nabla)y \cdot \upsilon - (\mu \Delta y) \cdot \upsilon + (My)\cdot \upsilon + (\nabla p)\cdot \upsilon = \mu \cdot \upsilon.$$

Then, integrating over Ω, we have

$$
\int\_{\Omega} (y\_t) \cdot v \, \mathrm{d}x + \int\_{\Omega} (y' \cdot \nabla) y \cdot v \, \mathrm{d}x - \int\_{\Omega} (\mu \Delta y) \cdot v \, \mathrm{d}x + 
$$

$$
\int\_{\Omega} (My) \cdot v \, \mathrm{d}x + \int\_{\Omega} (\nabla p) \cdot v \, \mathrm{d}x = \int\_{\Omega} u \cdot v \, \mathrm{d}x,
$$

therefore

$$\begin{aligned} \int\_{\Omega} (y\_t) \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x} + \int\_{\Omega} (y' \cdot \nabla) \boldsymbol{y} \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x} - \mu \int\_{\Omega} (\Delta \boldsymbol{y}) \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x} + \\\\ \int\_{\Omega} (M\boldsymbol{y}) \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x} + \int\_{\Omega} (\nabla p) \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x} = \int\_{\Omega} \boldsymbol{u} \cdot \boldsymbol{v} \, \mathrm{d}\mathbf{x}, \end{aligned} \tag{8}$$

then, applying Green's Theorem in the third term of the previous equation, we obtain

$$\begin{aligned} \int\_{\Omega} (y\_t) \cdot v \, \mathrm{d}x + \int\_{\Omega} (y' \cdot \nabla) y \cdot v \, \mathrm{d}x - \mu \left[ \int\_{\partial \Omega} \frac{\partial y}{\partial \eta} v \, \mathrm{d}\theta - \int\_{\Omega} \nabla y \cdot \nabla v \, \mathrm{d}x \right] + \\\\ + \int\_{\Omega} (My) \cdot v \, \mathrm{d}x + \int\_{\Omega} (\nabla p) \cdot v \, \mathrm{d}x = \int\_{\Omega} u \cdot v \, \mathrm{d}x, \end{aligned}$$

thus, we have

$$\begin{aligned} \int\_{\Omega} (y\_t) \cdot v \, \mathrm{d}x + \int\_{\Omega} (y' \cdot \nabla) y \cdot v \, \mathrm{d}x - \mu \int\_{\partial \Omega} \frac{\partial y}{\partial \eta} v \, \mathrm{d}\theta + \mu \int\_{\Omega} \nabla y \cdot \nabla v \, \mathrm{d}x + \mu \\\\ + \int\_{\Omega} (My) \cdot v \, \mathrm{d}x + \int\_{\Omega} (\nabla p) \cdot v \, \mathrm{d}x = \int\_{\Omega} u \cdot v \, \mathrm{d}x, \end{aligned}$$

Keeping in mind the boundary condition *y ∂*Ω = 0 and div (*y* ) = 0, we obtain the following:

$$\begin{aligned} \int\_{\Omega} (y\_t) \cdot v \, \mathrm{d}x + \int\_{\Omega} (y' \cdot \nabla) y \cdot v \, \mathrm{d}x - \mu \int\_{\partial \Omega} \frac{\partial y}{\partial \eta} v \, \mathrm{d}\theta + \\\\ + \mu \int\_{\Omega} \nabla y \cdot \nabla v \, \mathrm{d}x + \int\_{\Omega} (My) \cdot v \, \mathrm{d}x = \int\_{\Omega} u \cdot v \, \mathrm{d}x, \end{aligned}$$

then

$$\int\_{\Omega} (y\_t) \cdot v \, \mathrm{d}x + \int\_{\Omega} (y' \cdot \nabla) y \cdot v \, \mathrm{d}x + \mu \int\_{\Omega} \nabla y \cdot \nabla v \, \mathrm{d}x + \int\_{\Omega} (My) \cdot v \, \mathrm{d}x = \int\_{\Omega} u \cdot v \, \mathrm{d}x. \tag{9}$$

Now, (9) can be rewritten as

$$(y\_t, v)\_2 + ((y' \cdot \nabla)y, v)\_2 + \mu((\nabla y, \nabla v))\_2 + (My, v)\_2 = (u, v)\_2. \tag{10}$$

Next, we introduce the following bilinear and trilinear form for our weak formulation of the (10), where *<sup>a</sup>*(·, ·) : *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω)<sup>3</sup> <sup>×</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω)<sup>3</sup> −→ <sup>R</sup> and *<sup>b</sup>*(·, ·, ·) : <sup>V</sup> <sup>×</sup> <sup>V</sup> <sup>×</sup> <sup>V</sup> −→ <sup>R</sup> are defined:

$$a((y, v)) = \int\_{\Omega} \nabla y \cdot \nabla v \tag{11}$$

and

$$b(y, v, w) := (y' \cdot \nabla v, w)\_{\mathcal{Q}} = \sum\_{i, j = 1}^{3} \int\_{\Omega} y\_i' \frac{\partial v\_j}{\partial \mathbf{x}\_i} w\_j \, \mathrm{d}x \text{ for all } y, v, w \in \mathbb{V}. \tag{12}$$

Then, replacing identities (11) and (12) in (10), we obtain the following:

$$(y\_{t\prime}v)\_2 + ((y^{\prime} \cdot \nabla)y\_{\prime\prime}v)\_2 + \mu((\nabla y\_{\prime\prime}\nabla v))\_2 + (My\_{\prime\prime}v)\_2 = (u,v)\_2$$

$$a(y\_t, v)\_2 + b(y, y, v) + \mu a((y, v)) + (My, v)\_2 = (u, v)\_2.$$

In summary, we have

$$a(y\_t, v)\_2 + b(y, y, v) + \mu a((y, v)) + (My, v)\_2 = (u, v)\_2. \tag{13}$$

Equation (13) suggests the following weak formulation for our non-linear system given in (4), and which we express thus:

For every *<sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) and *<sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup>, we can find a solution *<sup>y</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*], <sup>V</sup>) with *yt* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) such that:

$$\begin{cases} \frac{d}{dt}(y,v)\_2 + \mu a((y,v))\_2 + b(y,y,v) + (My,v)\_2 = (u,v)\_2 \text{ for all } v \in \mathbb{V} \text{ a.e. } t \in (0,T) \\\\ y(0) = y\_0 \end{cases} \tag{14}$$

where the term

$$\begin{split} b(y, v, w) &= (y' \cdot \nabla v, w)\_2 \\ &= \int\_{\Omega} (y' \cdot \nabla) v \cdot w \, \mathrm{d}x \\ &= \sum\_{i, j=1}^{3} \int\_{\Omega} y\_i' \frac{\partial v\_j}{\partial \boldsymbol{x}\_i} w\_j \, \mathrm{d}x \text{ for all } y, v, w \in \mathbb{V}. \end{split} \tag{15}$$

corresponds to the non-linear term of our system (4).

We call the expression (14) the variational formulation (or weak) for our non-linear system given in (4).

On the other hand, let us see some properties with respect to the non-linear operator defined in (15), which can be found in [8]: For every (*y*, *<sup>v</sup>*, *<sup>w</sup>*) <sup>∈</sup> <sup>V</sup>, we have that:

$$\begin{array}{ll} 1. & b(y,v,w) = -b(y,w,v) \text{ if } y \cdot\_{\ldots} n = 0 \text{ on } \Gamma. \\ \end{array}$$

$$\text{2.}\quad b(y,v,v) = 0 \text{ if } y \cdot n = 0 \text{ on } \Gamma \quad \text{for all} \quad u,v \in \mathbb{V}.$$

As a consequence, we have the following lemma.

**Lemma 1** ([8,13])**.** *If n* = 3*, then*

$$|\boldsymbol{b}| \, |\boldsymbol{b}(\boldsymbol{y}, \boldsymbol{v}, \boldsymbol{w}) \, | \leq \begin{cases} \mathbb{C} \, |\, \boldsymbol{y} \mid \, ^{1/4} \cdot ||\boldsymbol{y}||^{3/4} \cdot ||\boldsymbol{v}|| \cdot |\boldsymbol{w} \mid ^{1/4} \cdot ||\boldsymbol{w}||^{3/4} \\\\ \mathbb{C} \, |\boldsymbol{y}| \, |\boldsymbol{v}|| \cdot ||\boldsymbol{w}|| \, |\boldsymbol{w}|| \, \text{ for all } \boldsymbol{y}, \boldsymbol{v}, \boldsymbol{w} \in \mathbb{V}, \end{cases} \tag{16}$$

*and in particular,*

$$b(y, v, w) \le \mathbb{C} \mid y \mid^{1/2} \cdot ||y||^{3/2} \cdot ||v|| \text{ for all } y, v \in \mathbb{V}.$$

In this way, we can introduce our definition of a weak solution given from (14) as we will see below.

#### *Existence and Uniqueness of Weak Solutions*

Let us introduce the definition of a weak solution.

**Definition 2.** *Given <sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) *and <sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup>*, say that <sup>y</sup> is a weak solution to the problem* (4) *on the interval* (0, *T*) *if:*

$$\begin{cases} y \in \mathcal{C}([0,T]; \mathbb{V}), \frac{dy}{dt} \in L^2([0,T], \mathbb{V}), \\\\ (y\_t(t), v)\_2 + \mu a((y(t), v))\_2 + (My(t), v)\_2 + b(y(t), y(t), v) = (u(t), v)\_2 \\\\ \text{for all } v \in \mathbb{V} \text{ a.e. } t \in (0, T), \\\\ y(0) = y\_0. \end{cases} \tag{17}$$

For our purposes, we want to give an equivalent formulation as an equation in functional spaces. For this, we can introduce a linear and continuous operator *A* : *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>) −→ *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) such that for all *<sup>y</sup>*, *<sup>v</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>) we have the following:

$$\begin{aligned} \langle Ay, v \rangle\_{L^2([0,T]; \mathcal{V}'), L^2([0,T]; \mathcal{V})} &= \int\_0^T \langle Ay(t), v(t) \rangle\_{\mathcal{V}', \mathcal{V}} \, \mathrm{d}t \\ &= \int\_0^T (y(t), v(t)) \, \mathrm{d}t \\ &= \int\_0^T \nabla y(t) \cdot \nabla v(t) \, \mathrm{d}t, \end{aligned} \tag{18}$$

and we define a non-linear operator *<sup>B</sup>* : *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) −→ *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) such that for all *<sup>y</sup>* <sup>∈</sup> *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>), *<sup>w</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>) we have

$$\begin{aligned} \langle B(y), w \rangle\_{L^2([0,T]; \mathcal{V}'), L^2([0,T]; \mathcal{V})} &= \int\_0^T \langle B(y)(t), w(t) \rangle\_{\mathcal{V}', \mathcal{V}} \, dt \\ &= \int\_0^T b(y(t), y(t), w(t)) \, dt. \end{aligned} \tag{19}$$

The operator *B* is a bounded mapping from *Wα*([0, *T*]; V) to *L*2([0, *T*]; V ), that is, it holds *B*(*y*)*L*2([0,*T*];V) <sup>≤</sup> *<sup>c</sup>y*<sup>2</sup> *<sup>W</sup>α*([0,*T*];V) for every *<sup>y</sup>* <sup>∈</sup> *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>).

Now, with the above notations, we can establish an equivalent formulation for our definition (17) in terms of the following functional differential equation.

**Definition 3.** *Let <sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) *and <sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup> *be given. A function <sup>y</sup>* <sup>∈</sup> *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) *is called a weak solution to the problem* (4) *on the interval* (0, *T*) *if it fulfills:*

$$\begin{cases} \mathcal{Y}\_t + \mu A \mathcal{Y} + B(\mathcal{y}) + M(\mathcal{y}) = \mathcal{u} \in L^2([0, T]; \mathbb{V}'), \\\\ \mathcal{Y}(0) = \mathcal{y}\_0 \in \mathbb{V}. \end{cases}$$

Denote by *<sup>A</sup>* the stokes operator, with domain *<sup>D</sup>*(*A*)=(*H*2(*D*))<sup>3</sup> <sup>∩</sup> <sup>V</sup>, defined by *Av* <sup>=</sup> <sup>−</sup>P(Δ*v*), for all *<sup>v</sup>* <sup>∈</sup> *<sup>D</sup>*(*A*), where <sup>P</sup> is the orthogonal projection that from *<sup>L</sup>*2(*D*)<sup>3</sup> to *H* and the set

$$a(y, v) = \int\_{\Omega} \nabla y \cdot \nabla v \, \mathrm{d}x \text{ for all } y, v \in \mathbb{V}\_{\prime}$$

or equivalently

$$(Ay, v) = a(y, v) \text{ for all } y, v \in \mathbb{V}.$$

The Stokes operator *<sup>A</sup>* is self-attached on *<sup>H</sup>*, with *<sup>A</sup>* <sup>∈</sup> *<sup>L</sup>*(V, <sup>V</sup> ), where V is the dual of V with the norm denoted by ·V and

$$(Ay, y) = \left\| y \right\|^2 \text{ for all } y \in \mathbb{V}.$$

Now, we can also equivalently formulate our control problem given in (6) and (7) using the operators defined in (18) and (19):

$$\left\{ \text{Minimize } J(y, u) \right\} = \frac{1}{2} \int\_0^T \left\| y - y\_d \right\|\_{L^2(\Omega)}^2 \, \text{d}t + \frac{\lambda}{2} \int\_0^T \left\| u - u\_d \right\|\_{L^2(\Omega)}^2 \, \text{d}t,\tag{20}$$

subject to the state equations

$$\begin{cases} y\_t + \mu Ay + B(y) + M(y) = u \in L^2([0, T]; \mathbb{V}'), \\ y(0) = y\_0 \in \mathbb{V}. \end{cases} \tag{21}$$

and the control of restrictions

$$
\mu \in \mathcal{U}\_{ad}.\tag{22}
$$

In this part, we will see how we can reduce our given energy functional (6). Sometimes it is convenient to work with the reduced functional since it allows us to better establish the study on the existence of the optimal values for our optimal control problem given in (6) and (7).

We can rewrite our optimal control problem as an optimization problem only in terms of *u*, as we will see below: we define the control to state mapping denoted by Υ, which associates an element *<sup>u</sup>* <sup>∈</sup> *Uad* <sup>⊂</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)4) with an element *<sup>y</sup>* <sup>∈</sup> *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) and which is the solution of (4).

The control to state mapping for the optimal control problem (6) and (7) is given as follows:

$$\begin{aligned} \mathcal{Y}: L^2([0, T]; L^2(\Omega)^4) &\longrightarrow \mathcal{W}^\mathfrak{a}([0, T]; \mathcal{V})\\ \mathfrak{u} &\longmapsto \mathcal{Y}(\mathfrak{u}) := \mathcal{y}\_{\mathfrak{u}\prime} \end{aligned} \tag{23}$$

where *yu* is the unique solution of (17).

#### **Remark 2.**

*1. Note that if we replace y* = Υ(*u*) *in our energy functional given in* (6)*, then our functional J would be expressed in terms of the control variable u, which we will denote by* Ψ*:*

$$J(u) := J(\mathcal{Y}(u), u) = \mathcal{Y}(u),\tag{24}$$

*where we have that*

$$\Psi(u) = \frac{1}{2} \int\_0^T \|\mathbf{Y}(u) - y\_d\|\_{L^2(\Omega)}^2 \, \mathrm{d}t + \frac{\lambda}{2} \int\_0^T \|u - u\_d\|\_{L^2(\Omega)}^2 \, \mathrm{d}t,$$

*and furthermore,* Ψ *is minimized on the set*

$$\mathcal{U}L\_{\text{ad}} = \{ u \in L^2(\Omega\_{\mathbb{T}}) : u\_{\mathfrak{e},i}(\mathbf{x},t) \le u\_i(\mathbf{x},t) \le u\_{\mathfrak{b},i}(\mathbf{x},t) \text{ a.e. on } \Omega\_{\mathbb{T}} \text{ for } i = 1,2,3,4\},\tag{25}$$

*where the term* Ψ(*u*) *will be called the reduced energy functional. In our context,* Υ *is the non-linear solution mapping associated with* (21)*.*

*2. The minimization of J subject to the state Equation* (21) *is equivalent to minimizing* Ψ *over all admissible controls.*

In order to understand the proof of the following theorem, which is one of the main results of this work, we give the proof in several stages. The main idea of the proof is the following. First, we use the Faedo–Galerkin method to find out approximate solutions of (4). Then, using some auxiliary estimations, we show the convergence of these approximations to the solution of the model (4).

**Theorem 1.** *For any <sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup>*, <sup>T</sup>* <sup>&</sup>gt; <sup>0</sup> *and <sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) *given, the problem* (4) *has a unique weak solution y on the interval* (0, *T*)*.*

**Proof. Stage 1:** Existence of the approximate solution.

Note that V defined by (5) is a reflexive and separable Hilbert space. Then, by a classical result of functional analysis, there is an orthonormal and dense subset (*zi*)*i*∈<sup>N</sup> of V.

Let us considerer the finite dimensional subspace

$$Z\_m = \text{span}\{z\_1, z\_2, z\_3, \dots, z\_m\} \subseteq \mathbb{V}, \text{ } \forall m \in \mathbb{N}.$$

Restricting to the space *Zm*, we solve the system of equations given from (17):

$$\begin{cases} \text{Let us find } y\_m(x,t) = \sum\_{i=1}^m a\_{im}(t)z\_i(x) \text{ such that} \\ (y\_m'(t), z\_j) + \mu((y\_m(t), z\_j)) + b(y\_m(t), y\_m(t), z\_j) + (My\_m(t), z\_j) = \langle u(t), z\_j \rangle, \\ t \in [0, T] : j = 1, 2, 3, \dots, m, \\ y\_m(0) = y\_{0m}. \end{cases} \tag{26}$$

Here, *<sup>y</sup>*0*<sup>m</sup>* is the orthogonal projection of the initial data *<sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup> on the subspace *Vm* = span{*z*1, *z*2, *z*3,..., *zm*}.

We can observe that (26) is in fact a Cauchy initial value problem for a non-linear system of ordinary differential equations. Indeed, in the unknowns *αim*(*t*), the system can be written as

$$\sum\_{i=1}^{m} (z\_i, z\_j) a'\_{im} + \mu \sum\_{i=1}^{m} ((z\_i, z\_j)) a\_{im} + \sum\_{i,k=1}^{m} b(z\_i, z\_k, z\_j) a\_{im} a\_{km} + \sum\_{i=1}^{m} (Mz\_i, z\_j) a\_{im} = \langle u(t), z\_j \rangle.$$

Now, due to the smoothness of the coefficients, we can use a classical result from the theory of ordinary differential equations, and ensure that there exists a unique classical solution *ym* defined on a maximal interval [0, *tm*] with 0 <sup>&</sup>lt; *tm* <sup>≤</sup> *<sup>T</sup>*. For the convergence, we need to show that *tm* = *T* for all *m*. In that way, the interval of existence of solutions will not change when *<sup>m</sup>* goes to infinity. If *tm* <sup>&</sup>lt; *<sup>T</sup>*, then

$$\lim\_{t \to t\_m} \sup \parallel y\_m(t) \parallel = +\infty. \tag{27}$$

In the following stage, we prove that (*ym*(*t*))*m*∈<sup>N</sup> is bounded on [0, *<sup>T</sup>*] by a constant independent of *t* and *m*. Then, the solution is defined in [0, *T*] for all *m*.

**Stage 2:** Estimates for the approximate solution.

Next, we show that (*ym*)*m*∈<sup>N</sup> is uniformly bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>). Indeed, multiplying equation (26) by *αjm*(*t*) and summing from *j* = 1 to *m*, we obtain

$$
\langle \boldsymbol{y}\_m'(t), \boldsymbol{y}\_m(t) \rangle + \mu (\langle \boldsymbol{y}\_m(t), \boldsymbol{y}\_m(t) \rangle) + b(\boldsymbol{y}\_m(t), \boldsymbol{y}\_m(t), \boldsymbol{y}\_m(t)) + (M\boldsymbol{y}\_m(t), \boldsymbol{y}\_m(t)) = \langle \boldsymbol{u}(t), \boldsymbol{y}\_m(t) \rangle.
$$

Using the Cauchy–Schwarz inequality and keeping in mind that *b*(*y*(*t*), *y*(*t*), *y*(*t*)) = 0 for almost all *t* ∈ [0, *T*], we obtain

$$\begin{split} \frac{1}{2} \frac{\text{d}}{\text{d}t} \mid y\_{\mathfrak{m}}(t) \mid^{2} + \mu \parallel y\_{\mathfrak{m}}(t) \parallel^{2} = \langle u(t), y\_{\mathfrak{m}}(t) \rangle - (My\_{\mathfrak{m}}(t), y\_{\mathfrak{m}}(t)) \\ \leq & \left| \langle u(t), y\_{\mathfrak{m}}(t) \rangle - (My\_{\mathfrak{m}}(t), y\_{\mathfrak{m}}(t)) \right| \\ \leq & \left| u(t), y\_{\mathfrak{m}}(t) \right| + \left| \left. M y\_{\mathfrak{m}}(t), y\_{\mathfrak{m}}(t) \right| \right| \\ \leq & \left| u(t) \right|\_{L^{2}([0,T], \mathbb{V}^{\mathfrak{p}})} \parallel y\_{\mathfrak{m}}(t) \parallel + \left| \left. M y\_{\mathfrak{m}}(t) \right| \right| \left| y\_{\mathfrak{m}}(t) \right| \\ \leq & \left| u(t) \right|\_{L^{2}([0,T], \mathbb{V}^{\mathfrak{p}})} \parallel y\_{\mathfrak{m}}(t) \right| + \left| \left. M \right| \left| \left. y\_{\mathfrak{m}}(t) \right| \right|. \end{split}$$

Therefore, we obtain

$$\frac{\text{d}}{\text{dt}} \mid y\_m(t) \parallel^2 + \mu \parallel y\_m(t) \parallel^2 \le \frac{||\,\mu(t)\parallel^2 + |\,\underline{y\_m(t)}\parallel^2 + ||\,\underline{M}\parallel|| \,\underline{y\_m(t)}\parallel^2. \tag{28}$$

The last inequality implies that

$$\frac{\mathbf{d}}{\mathbf{d}t} \mid y\_m(t) \mid^2 \le |\boldsymbol{u}(t)|^2 + c\_0 \mid y\_m(t) \mid^2, \quad c\_0 = \max\{1, ||\boldsymbol{M}||\}.$$

Integrating on both sides of the previous expression and using the Grönwall's inequality, we deduce that

$$\begin{split} \left| \left| \left| y\_{m}(t) \right| \right|^{2} - \left| \left| y\_{m}(0) \right| \right|^{2} \leq \left| \left| u \right| \right|^{2} t + c\_{0} \int\_{0}^{t} \left| \left| y\_{m}(t) \right| \right|^{2} \text{ d}t \\ \leq \left| \left| y\_{m}(t) \right| \right|^{2} \leq \left| \left| y\_{m}(0) \right| \right|^{2} + \left| \left| u \right| \right|^{2} t + c\_{0} \int\_{0}^{t} \left| \left| y\_{m}(t) \right| \right|^{2} \text{ d}t \\ \leq \left| \left| y\_{m}(t) \right| \right|^{2} \leq \left( \left| \left| y\_{0} \right| \right|^{2} + \left| \left| u \right| \right|^{2} t \right) e^{\int\_{0}^{t} c\_{0} \left| \text{d}t} \\ \leq \left| \left| y\_{0} \right| \right|^{2} \leq \left( \left| \left| y\_{0} \right| \right|^{2} + \left| \left| u \right| \right|^{2} t \right) e^{\varepsilon\_{0} t} . \end{split}$$

It follows that

$$\sup\_{t \in [0, t\_m]} |\left| y\_m(t) \right|^2 \le \left( \left| \left| y\_0 \right|^2 + \left| \left| u \right| \right|^2 t \right) e^{c\_0 t}. \tag{29}$$

This tell us, that *tm* <sup>=</sup> *<sup>T</sup>* for all *<sup>m</sup>*. Moreover, (*ym*)*m*∈<sup>N</sup> is uniformly bounded in *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>H</sup>*).

On other hand, integrating from 0 to *T* on both sides of Equation (28) and using (29), we obtain

$$\|\|\|y\_m(T)\|\|^2 + \mu \int\_0^T \|\|y\_m(t)\|\|^2 \, dt \le \left( \|\|y\_0\|\|^2 + \|\|\mu\|\|^2 \, T + c\_0(\|y\_0\|^2 + \|\|\mu\|\|^2 \, T)e^{c\_0T} \right),$$

therefore,

$$\mu \int\_0^T \parallel y\_m(t) \parallel^2 \, dt \le \left( \mid y\_0 \mid^2 + \mid \mu \mid^2 \, T + \varepsilon\_0 (|y\_0|^2 + \mid \mu \mid^2 T) e^{c\_0 T} \right).$$

Consequently, (*ym*)*m*∈<sup>N</sup> is uniformly bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>).

**Stage 3:** Estimates for the derivative of approximate solution.

Now, we prove that the derivative *y <sup>m</sup>*(*t*) is uniformly bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; V). Indeed, multiplying (26) by *α jm*(*t*) and summing from *j* = 1 to *m*, we obtain

$$|\left(y\_m'(t)\right)^2 + \mu((y\_m(t), y\_m'(t))) + b(y\_m(t), y\_m(t), y\_m'(t)) + ((My\_m(t), y\_m'(t)) = \langle u(t), y\_m'(t) \rangle.$$

After integrate the last equation, we obtain

$$\begin{aligned} \int\_0^t |\, \boldsymbol{y}'\_m(t) \mid \, ^2 \, \mathrm{d}t + \mu \int\_0^t ((\boldsymbol{y}\_m(t), \boldsymbol{y}'\_m(t)) \, \mathrm{d}t + \int\_0^t ((\boldsymbol{M}\boldsymbol{y}\_m(t), \boldsymbol{y}'\_m(t))) \, \mathrm{d}t + \mu \\ + \int\_0^t \boldsymbol{b}(\boldsymbol{y}\_m(t), \boldsymbol{y}\_m(t), \boldsymbol{y}'\_m(t)) \, \mathrm{d}t = \int\_0^t \langle \boldsymbol{u}(t), \boldsymbol{y}'\_m(t) \rangle \, \mathrm{d}t. \end{aligned}$$

Therefore

$$\begin{aligned} \int\_0^t |\, |y\_m'(t)| \, ^2 \, \mathrm{d}t + \mu \int\_0^t (\langle y\_m(t), y\_m'(t) \rangle \, \mathrm{d}t + \int\_0^t (\langle My\_m(t), y\_m'(t) \rangle) \, \mathrm{d}t \\ &= \int\_0^t \langle u(t), y\_m'(t) \rangle \mathrm{d}t - \int\_0^t b(y\_m(t), y\_m(t), y\_m'(t)) \, \mathrm{d}t. \end{aligned}$$

Taking into account the Cauchy's and Ladyzhenskaya's inequalities, we have that

 *t* 0 | *y <sup>m</sup>*(*t*) | <sup>2</sup> <sup>d</sup>*<sup>t</sup>* <sup>+</sup> *<sup>μ</sup>* 2 *t* 0 ((*ym*(*t*), *y <sup>m</sup>*(*t*))) d*t* + *t* 0 ((*Mym*(*t*), *y <sup>m</sup>*(*t*))) d*t* = = *t* 0 *u*(*t*), *y <sup>m</sup>*(*t*) d*t* − *t* 0 *b*(*ym*(*t*), *y <sup>m</sup>*(*t*), *ym*(*t*)) d*t* ≤ *<sup>u</sup>*(*t*) <sup>2</sup> *L*2([0,*T*];V) *y <sup>m</sup>*(*t*) *L*2([0,*T*];V) + *t* 0 | *ym*(*t*) | 2 *L*4([0,*T*];V) | ∇*y <sup>m</sup>*(*t*) | d*t* ≤ 1 *<sup>α</sup>*<sup>2</sup> *<sup>u</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) <sup>+</sup>*α*<sup>2</sup> <sup>4</sup> *<sup>y</sup> <sup>m</sup>*(*t*) <sup>2</sup> *L*2([0,*T*];V) + *c t* 0 | *ym*(*t*) | 1/2| ∇*ym*(*t*) <sup>|</sup> 3/2 ∗|∇*y <sup>m</sup>*(*t*) | d*t* ≤ 1 *<sup>α</sup>*<sup>2</sup> *<sup>u</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) <sup>+</sup>*α*<sup>2</sup> <sup>4</sup> *<sup>y</sup> <sup>m</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) +*c t* 0 | ∇*ym*(*t*) || ∇*y <sup>m</sup>*(*t*) | d*t* ≤ 1 *<sup>α</sup>*<sup>2</sup> *<sup>u</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) <sup>+</sup>*α*<sup>2</sup> <sup>4</sup> *<sup>y</sup> <sup>m</sup>*(*t*) <sup>2</sup> *L*2([0,*T*];V) +*c t* 0 | ∇*ym*(*t*) | <sup>2</sup> d*t* 1/2 *<sup>t</sup>* 0 | ∇*y <sup>m</sup>*(*t*) | <sup>2</sup> d*t* 1/2 ≤ 1 *<sup>α</sup>*<sup>2</sup> *<sup>u</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) <sup>+</sup>*α*<sup>2</sup> <sup>4</sup> *<sup>y</sup> <sup>m</sup>*(*t*) <sup>2</sup> *L*2([0,*T*];V) <sup>+</sup> *<sup>c</sup> ym*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) <sup>+</sup>*α*<sup>2</sup> <sup>4</sup> *<sup>y</sup> <sup>m</sup>*(*t*) <sup>2</sup> *<sup>L</sup>*2([0,*T*];V) .

In summary, we have the following:

$$\begin{split} 2\int\_{0}^{t} |\, |\, y\_{m}'(t) \mid \, ^{2} \, \mathrm{d}t + \mu \int\_{0}^{t} ((\!(y\_{m}(t), y\_{m}'(t))) \, \mathrm{d}t + \int\_{0}^{t} ((\!(My\_{m}(t), y\_{m}'(t))) \, \mathrm{d}t \\ \leq \frac{2}{\alpha^{2}} \left\| \, \left\| \, u(t) \right\|\right\|\_{L^{2}([0,T];\mathcal{V}')}^{2} + c \left\| \, \left\| \, y\_{m}(t) \right\|\right\|\_{L^{2}([0,T];\mathcal{V})}^{2} \text{ for all } t \in [0, T]. \end{split}$$

In this way, due to the fact that *ym*(*t*) is uniformly bounded in *L*∞([0, *T*]; V), it follows that *y <sup>m</sup>*(*t*) is uniformly bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; V).

**Stage 4:** Extraction of subsequence and convergence to the solution.

We can extract a subsequence of (*ym*)*m*∈<sup>N</sup> that converges (in an appropriate sense) to a function *y* and then go to the limit in the approximate problem given by (26) as follows:

Since (*ym*)*m*∈<sup>N</sup> is uniformly bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>) <sup>∩</sup> *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>H</sup>*), then there is a subsequence (which we will denote in the same way) (*ym*)*m*∈<sup>N</sup> such that

$$\begin{cases} y\_{\mathfrak{m}} \rightharpoonup y \text{ weak in } L^2([0,T]; \mathbb{V})\\ y\_{\mathfrak{m}} \rightharpoonup y \text{ weak:} \ast \text{ in } L^\infty([0,T]; H) \\ \frac{\partial y\_{\mathfrak{m}}}{\partial t} \rightharpoonup \frac{\partial y}{\partial t} \text{ weakly in } L^2([0,T]; \mathbb{V}). \end{cases} \tag{30}$$

Now, note that <sup>V</sup> -→ *H* with dense and continuous injections.

**Stage 5:** Existence of solutions.

Let us take a function *η* ∈ D(0, *T*). In Equation (26), we multiply by *η*, and integrate over the interval (0, *T*). Then we obtain,

$$\begin{aligned} -\int\_0^T \left( y\_m(t), \eta'(t) z\_j \right) \, \mathrm{d}t + \mu \int\_0^T \left( (y\_m(t), \eta(t) z\_j) \right) + \int\_0^T b(y\_m(t), y\_m(t), z\_j \eta(t)) \, \mathrm{d}t \\ + \int\_0^T \left( (M y\_m(t), z\_j \eta(t)) \right) = \int\_0^T \left( u(t), z\_j \eta(t) \right) \, \mathrm{d}t. \end{aligned}$$

Consequently, we have that,

$$\begin{aligned} -\int\_0^T \left( y(t), \eta'(t)z\_j \right) \, \mathrm{d}t + \mu \int\_0^T \left( (y(t), \eta(t)z\_j) \right) + \int\_0^T b(y(t), y(t), \eta(t)z\_j) \, \mathrm{d}t \\ + \int\_0^T \left( (My(t), \eta(t)z\_j) \right) = \int\_0^T \left( u(t), \eta(t)z\_j \right) \, \mathrm{d}t, \end{aligned}$$

or equivalently

$$\begin{aligned} \int\_0^T \left( y'(t), \eta(t) z\_j \right) \, \mathrm{d}t + \mu \int\_0^T \left( (y(t), \eta(t) z\_j) \right) + \int\_0^T b(y(t), y(t), \eta(t) z\_j) \, \mathrm{d}t \\ + \int\_0^T \left( (My(t), \eta(t) z\_j) \right) = \int\_0^T \left( u(t), \eta(t) z\_j \right) \, \mathrm{d}t. \end{aligned}$$

This equality is true by linearity and by density for all *<sup>v</sup>* <sup>∈</sup> <sup>V</sup>. Thus, we have that *<sup>y</sup>* verifies the equation given by (17).

Now, let us show that *<sup>y</sup>*(0) = *<sup>y</sup>*0. Since *<sup>y</sup>* is a weak solution of (4), taking *<sup>η</sup>* <sup>∈</sup> *<sup>C</sup>*∞([0, *<sup>T</sup>*]) with *<sup>η</sup>*(*T*) = 0, for all *<sup>v</sup>* <sup>∈</sup> <sup>V</sup> we have that:

$$\begin{aligned} -\int\_0^T \left( y(t), \eta'(t)v \right) \, \mathrm{d}t + \mu \int\_0^T \left( (y(t), \eta(t)v) \right) + \int\_0^T b(y(t), \eta(t), \eta(t)v) \, \mathrm{d}t \\ + \int\_0^T \left( (My(t), \eta(t)v) \right) = \left( u\_0, \eta(0)v \right) + \int\_0^T \left( u(t), \eta(t)v \right) \, \mathrm{d}t. \end{aligned}$$

On the other hand,

$$\frac{\mathbf{d}}{\mathbf{d}t} \langle y, \eta(v) \rangle = \langle y\_{t\prime} \eta(v) \rangle + \langle y\_{\prime\prime} \eta\_t(v) \rangle\_{\ast}$$

we can integrate from 0 a *T* and we obtain the following:

$$\begin{aligned} \langle \langle \boldsymbol{y}(t), \boldsymbol{\eta}(t)(\boldsymbol{v}) \rangle - \langle \boldsymbol{y}(0), \boldsymbol{\eta}(0) \rangle &= \int\_0^T \left[ \langle \boldsymbol{y}\_t, \boldsymbol{\eta}(\boldsymbol{v}) \rangle + \langle \boldsymbol{y}, \boldsymbol{\eta}\_t(\boldsymbol{v}) \rangle \right] \, \mathrm{d}t \\ &= \int\_0^T \langle \boldsymbol{y}\_t, \boldsymbol{\eta}(\boldsymbol{v}) \rangle \, \mathrm{d}t + \int\_0^T \langle \boldsymbol{y}\_t \boldsymbol{\eta}\_t(\boldsymbol{v}) \rangle \, \mathrm{d}t. \end{aligned}$$

Therefore,

$$\int\_0^T \langle y, \eta\_t(v) \rangle \, \mathrm{d}t = \langle y(t), \eta(t)(v) \rangle - \langle y(0), \eta(0) \rangle - \int\_0^T \langle y\_{t\prime} \eta(v) \rangle \, \mathrm{d}t.$$

Thus, if *η* is such that *η*(*T*) = 0, it follows:

$$\int\_0^T \langle y, \eta\_t(v) \rangle \, \mathrm{d}t = -\langle y(0), \eta(0)v \rangle - \int\_0^T \langle y\_t, \eta(v) \rangle \, \mathrm{d}t.$$

Thus, we obtain that

$$\eta(0)(y\_0 - y(0), v) = 0 \text{ for all } \eta \in \mathbb{C}^\infty([0, T]), \ v \in \mathbb{V}.$$

Taking *<sup>η</sup>*(0) = 1, obtain that *<sup>y</sup>*(0) = *<sup>y</sup>*<sup>0</sup> <sup>∈</sup> <sup>V</sup>.

**Stage 6:** Continuity and uniqueness of the weak solution.

Let us show that *<sup>y</sup>* : [0, *<sup>T</sup>*] −→ <sup>V</sup> is continuous. Indeed, note that the non-linear operator *B*(*y*) is defined by (19), now since *yt* = *u* − *μAy* − *B*(*y*), in this way, by the Lemma (1) it follows that *<sup>B</sup>*(*y*) <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ) and also *<sup>u</sup>*, *Ay*, *My* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ), therefore, *yt* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ), thus *<sup>y</sup>* <sup>∈</sup> *<sup>C</sup>*([0, *<sup>T</sup>*]; <sup>V</sup>).

Now, suppose that there exist two weak solutions *y*1(*t*), *y*2(*t*) for the equations in the dynamics of the stratified fluids given by (4), with the initials values *<sup>y</sup>*<sup>01</sup> , *<sup>y</sup>*<sup>02</sup> <sup>∈</sup> <sup>V</sup>. Let us denote by *<sup>y</sup>*(*t*) = *<sup>y</sup>*1(*t*) <sup>−</sup> *<sup>y</sup>*2(*t*), then *<sup>y</sup>*(*t*) satisfies *yt*(*t*) <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup> ), and for every test function *δ*, we have the following:

$$
\langle \underline{\boldsymbol{y}}\_{t\prime} \boldsymbol{\delta} \rangle + \mu (\nabla \underline{\boldsymbol{y}}\_{t\prime} \nabla \boldsymbol{\delta}) + (M \underline{\boldsymbol{y}}\_{t\prime} \boldsymbol{\delta}) + b (\underline{\boldsymbol{y}}\_{1\prime} \underline{\boldsymbol{y}}\_{1\prime} \boldsymbol{\delta}) - b (\underline{\boldsymbol{y}}\_{2\prime} \underline{\boldsymbol{y}}\_{2\prime} \boldsymbol{\delta}) = 0.
$$

In this way, taking *δ* = *y*, it follows that:

$$\frac{1}{2}\frac{d}{dt}\left|\left|y(t)\right|^2 + \mu \left|\left|\nabla y(t)\right|\right|^2 + \left|\left|My(t)\right|\right|^2 + b(y\_1(t), y\_1(t), y(t)) - b(y\_2(t), y\_2(t), y(t)) = 0.$$

Now, adding and subtracting *b*(*y*1(*t*), *y*2(*t*), *y*(*t*)) and having in mind that *b*(*y*1(*t*), *y*(*t*), *y*(*t*)) = 0, we obtain the following:

$$\begin{split} \frac{d}{dt} \left| |y(t)|^2 + |\operatorname{M}y(t)|^2 + 2\mu \mid \nabla y(t) \right|^2 &= 2(-b(y\_1(t), y\_1(t), y(t)) + b(y\_2(t), y\_2(t), y(t))) \\ &= -2b(y(t), y\_2(t), y(t)) \\ &= -2\sum\_{i,j=1}^3 \int\_{\Omega} y^i \frac{\partial y\_j^j}{\partial t^i} y^j \operatorname{d}t \\ &\leq 2\sum\_{i,j=1}^3 \left( \int\_{\Omega} |y^i|^4 \operatorname{d}t \right)^{1/4} \left( \int\_{\Omega} |\nabla y\_j^j|^2 \operatorname{d}t \right)^{1/2} \left( \int\_{\Omega} |y^i|^4 \operatorname{d}t \right)^{1/4} \\ &\leq 2\left|\nabla y\_2(t)\right| \|y(t)\|\_2^2 \\ &\leq 2c\left|\nabla y\_2(t)\right| y(t)^{1/2}|\nabla y(t)|^{5/2} \\ &\leq c\left(|y(t)|^{1/2}|\nabla y(t)|\right)^{3/2} \\ &\leq c(\mu) (|y(t)|^{1/2})^4 + 2\mu (|\nabla y(t)|^{3/2})^{4/3} \\ &\leq c\left(|y(t)|^{2} + 2\mu \mid \nabla y(t) \right)^2. \end{split}$$

In summary, we have that

$$\frac{\text{d}}{\text{d}t} \mid y(t) \mid^2 + \mid My(t) \mid^2 + 2\mu \mid \nabla y(t) \mid^2 \le c \mid y(t) \mid^2 + 2\mu \mid \nabla y(t) \mid^2 \dots$$

This implies that

$$\frac{\mathbf{d}}{\mathbf{d}t} \mid y(t) \mid^2 + \mid My(t) \mid^2 \le c \mid y(t) \mid^2 \dots$$

Now, since *<sup>y</sup>*2(*t*) <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>), we can apply Gronwall's Lemma. Using the fact that *y*(0) = 0, it follows

$$|\;y(t)\;|^2 \le 0 \text{ for all } t \in [0,T]\_\prime$$

thus, *y*<sup>1</sup> = *y*2, as we wanted to prove.

#### **4. Study of the Existence of Solutions for Our Optimal Control Problem**

In this section, we will show the existence of optimal controls for our non-linear system given by (7).

Let us show that our optimal control problem formulated in (6) and (7) has a solution in *Uad*. To prove this fact, we need the following result associated with the non-linear operator.

**Lemma 2** ([14])**.** *Assume that yn converges to <sup>y</sup> in <sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) *weakly for* <sup>1</sup> <sup>≤</sup> *<sup>α</sup>* <sup>≤</sup> <sup>2</sup>*. Then, for any v* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>)*,*

*B*(*yn*), *v <sup>L</sup>*2([0,*T*];V),*L*2([0,*T*];V) −→ *B*(*y*), *v <sup>L</sup>*2([0,*T*];V),*L*2([0,*T*];V) *as n* −→ ∞.

Now, with the following result, we want to show that our optimal control problem formulated by (6) and (7) has a solution in *Uad*.

**Theorem 2.** *The optimal control problem given by* (6) *and* (7) *admits an optimal solution u* ∈ *Uad with associated state <sup>y</sup>* <sup>∈</sup> *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) *for* <sup>1</sup> <sup>≤</sup> *<sup>α</sup>* <sup>≤</sup> <sup>2</sup>*.*

**Proof.** Note that the set of admissible controls defined by (3) is non-empty, convex, and closed in *<sup>L</sup>*2(Ω)4. Then, for every control *<sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*2(Ω)4, applying Theorem 1, there is a unique weak solution of the state Equations (20) and (21). Therefore, we have that *J*(*y*, *u*) ≥ 0 for every admissible pair (*y*, *u*).

Hence, there exists the infimum of *J* over all admissible controls and states that such:

$$J(y\_n, u\_n) = \int\_{\Omega} ||y\_n - y\_d||^2 \,\mathrm{d}t + \lambda \int\_{\Omega} ||u\_n - u\_d||^2 \,\mathrm{d}t\_n$$

then,

$$\begin{aligned} \|\boldsymbol{u}\_{n}\|^{2} &= \|(\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}) + \boldsymbol{u}\_{d}\|^{2} \\ &\leq \left(\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\| + \|\boldsymbol{u}\_{d}\|\right)^{2} \\ &= \left\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\right\|^{2} + 2\left\|\boldsymbol{u}\_{d}\right\| \cdot \left\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\right\| + \left\|\boldsymbol{u}\_{d}\right\|^{2} \\ &\leq \left\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\right\|^{2} + \left\|\boldsymbol{u}\_{d}\right\|^{2} + \left\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\right\|^{2} + \left\|\boldsymbol{u}\_{d}\right\|^{2} \\ &= 2\left\|\boldsymbol{u}\_{n} - \boldsymbol{u}\_{d}\right\|^{2} + 2\left\|\boldsymbol{u}\_{d}\right\|^{2} .\end{aligned}$$

now it follows that

$$\begin{aligned} \|u\_n\|\_{L^2([0,T];V)} < \infty &\implies \int\_0^T \|u\_n(\mathbf{x},t)\|^2 \, \mathrm{d}t \le 2\int\_0^T \|u\_n - u\_d\|^2 \, \mathrm{d}t + 2\|u\_d\|^2 T\\ &\le f(u\_n, u\_n) + 2\|u\_d\|^2 T\\ &< \infty. \end{aligned}$$

In summary, we have that:

$$0 \le \overline{\lambda} = \inf\_{(y,\mathfrak{u})\text{ admissible}} J(y,\mathfrak{u})$$

$$ < \infty.$$

On the other hand, there is a sequence (*yn*, *un*)*n*∈<sup>N</sup> of admissible pairs such that *<sup>J</sup>*(*yn*, *un*) −→ *J* as *n* −→ +∞.

First, we will show that (*un*)*n*∈<sup>N</sup> and (*yn*)*n*∈<sup>N</sup> are bounded sequences in *<sup>L</sup>*2(Ω)<sup>4</sup> and *Wα*([0, *T*]; V), respectively.

From the convergence, we see that the set (*J*(*yn*, *un*))*n*∈<sup>N</sup> is bounded, this implies that the set (*un*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*2(Ω)4.

Now, we need to show that (*yn*)*n*∈<sup>N</sup> and (*ynt*)*n*∈<sup>N</sup> are bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>). Indeed,

$$\begin{cases} y\_{\text{nt}}(t) + \mu A y\_{\text{n}}(t) + M y\_{\text{n}}(t) + B(y\_{\text{n}}(t)) = \mu\_{\text{n}}(t) \text{ in } L^2([0, T]; \mathbb{V}')\\ \qquad \qquad y\_{\text{n}}(0) = y\_0 \text{ in } \mathbb{V}. \end{cases} \tag{31}$$

Since *yn*(*t*) <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>), we multiply (31) by *<sup>χ</sup>*(0, *<sup>t</sup>*)*yn*(*t*), where *<sup>χ</sup>*(0, *<sup>t</sup>*) is the characteristic function on the interval (0, *t*) and applying the identity *b*(*yn*(*t*), *yn*(*t*), *yn*(*t*)) = 0, we get

$$\int\_{0}^{t} \left( y\_{n\ell}(t), y\_{n}(t) \right) \mathrm{d}t + \mu \int\_{0}^{t} \left\| y\_{n}(t) \right\|^{2} \mathrm{d}t + \int\_{0}^{t} \left( \left( M(y\_{n}(t), y\_{n}(t)) \right) \mathrm{d}t = \int\_{0}^{t} \left\langle u\_{n}(t), y\_{n}(t) \right\rangle \mathrm{d}t,\tag{32}$$

The right-hand side can be estimated by

$$\begin{split} \left| \int\_{0}^{t} (u\_{n}(t), y\_{n}(t)) \, \mathrm{d}t \right| &\leq \int\_{0}^{t} |\, u\_{n}(t) \mid \cdot \mid y\_{n}(t) \mid \, \mathrm{d}t \\ &\leq \mathbb{C} \int\_{0}^{t} \left\| y\_{n}(t) \right\| \cdot \left| \, u\_{n}(t) \right| \, \mathrm{d}t \\ &\leq \frac{\mu}{2} \int\_{0}^{t} \left\| y\_{n}(t) \right\|^{2} \, \mathrm{d}t + \frac{\mathcal{C}^{2}}{2\mu} \int\_{0}^{t} \left| \, u\_{n}(t) \right| \, \mathrm{^{2}} \, \mathrm{d}t, \end{split} \tag{33}$$

where *C* only depends on Ω.

Then, applying the integration in parts to (32) and (33) implies that:

$$\begin{split} \left( \left| \left| y\_n(t) \right| \right|^2 + \mu \int\_0^t \left\| y\_n(t) \right\|^2 \, \mathrm{d}t + \int\_0^t \left( \left( M(y\_n(t), y\_n(t)) \right) \right) \, \mathrm{d}t \leq \left| y\_n(0) \right|^2 + \frac{\mathcal{C}^2}{\mu} \int\_0^t \left| u\_n(t) \right|^2 \, \mathrm{d}t \\ \leq \left| \left. y\_0 \right|^2 + \frac{\mathcal{C}^2}{\mu} \left\| u\_n \right\|\_{L^2(\Omega)^4}. \end{split} \right. \\ \left. \left. \left( \left| \left. y\_n \right| \right| \right) \right| + \left| \left. y\_n \right| \right|\_{L^2(\Omega)^4}. \end{split}$$

Since (*un*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*2(Ω)<sup>4</sup> and (*yn*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>H</sup>*) and *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>V</sup>*), it follows that (*yn*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>V</sup>*).

Now, multiplying (31) by *ynt* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>V</sup>*), we obtain the following:

$$\begin{split} \|y\_{nt}\|\_{L^{2}(\Omega)^{4}}^{2} + \mu \int\_{0}^{T} ((y\_{n}(t), y\_{nt}(t))) \, \mathrm{d}t + \int\_{0}^{T} (M(y\_{n}(t), y\_{nt}(t)) \, \mathrm{d}t + \\ + \int\_{0}^{T} b(y\_{n}(t), y\_{n}(t), y\_{nt}(t)) \, \mathrm{d}t = \int\_{0}^{T} (u\_{n}(t), y\_{nt}(t)) \, \mathrm{d}t. \end{split} \tag{34}$$

Thus

$$\begin{aligned} \|y\_{nt}\|\_{L^2(\Omega)^4}^2 + \mu \int\_0^T ((y\_n(t), y\_{nt}(t))) \, \mathrm{d}t + \int\_0^T (M(y\_n(t), y\_{nt}(t)) \, \mathrm{d}t \\ = \int\_0^T (u\_n(t), y\_{nt}(t)) \, \mathrm{d}t - \int\_0^T b(y\_n(t), y\_n(t), y\_{nt}(t)) \, \mathrm{d}t. \end{aligned}$$

On the side right, we have the following estimate:

$$\begin{aligned} \left| \int\_0^T (\mu\_n(t), y\_{nt}(t)) \, \mathrm{d}t \right| &\leq \int\_0^T \left( \frac{|\, \mu\_n(t) \, |^2}{4} + |\, y\_{nt}(t) \, |^2 \right) \mathrm{d}t \\ &= \frac{1}{4} \| |\mu\_n| \|\_{L^2(\Omega)^4}^2 + \| y\_{nt} \|\_{L^2(\Omega)^4}^2 \end{aligned}$$

Note that

$$\mu \int\_0^T \left( (y\_n(t), y\_{nt}(t)) \right) \, \mathrm{d}t = \frac{\mu}{2} \left( \|y\_n(T)\|^2 - \|y\_0\|^2 \right).$$

Then, since (*yn*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>H</sup>*) and *<sup>L</sup>*∞([0, *<sup>T</sup>*]; *<sup>V</sup>*), it follows that

$$\begin{aligned} \left| \int\_0^T b(y\_n(t), y\_n(t), y\_{nt}(t)) \, \mathrm{d}t \right| &\leq \mathbb{C} \int\_0^T |\, y\_n(t) \, | \, ^{1/2} \cdot \| y\_n(t) \| ^{3/2} \cdot \| y\_{nt}(t) \| \, \mathrm{d}t \\ &\leq \mathbb{C} \int\_0^T \| y\_{nt}(t) \| \, \mathrm{d}t. \end{aligned}$$

Therefore, from Equation (34), we have

$$\frac{\mu}{2} \|y\_n(T)\|^2 + \int\_0^t \left( \left(M(y\_n(t), y\_n(t))\right) \, \mathrm{d}t + \mathbb{C} \int\_0^T \|y\_{nt}(t)\| \, \mathrm{d}t \le \frac{\mu}{2} \|y\_0\|^2 + \frac{1}{4} \|u\_{tt}\|\_{L^2(\Omega)}^2 \right)$$

Since (*un*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*2(Ω)4, we can ensure that (*ynt*)*n*∈<sup>N</sup> is bounded in *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>V</sup>*). Thus, it follows that (*yn*)*n*∈<sup>N</sup> is bounded in *<sup>W</sup>α*([0, *<sup>T</sup>*]; *<sup>V</sup>*).

Then, we can extract a subsequence (*y <sup>n</sup>*, *u <sup>n</sup>*)*n*∈<sup>N</sup> converging weakly in the space *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) <sup>×</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> to some limit (*y*, *<sup>u</sup>*).

Now, let us show that (*y*, *u*) is an admissible pair, that is, it satisfies the state equations given by (21). Indeed, note that the set of admissible controls *Uad* is non-empty, convex, and closed in *<sup>L</sup>*2(Ω)3, so it is weakly closed. Therefore, *<sup>u</sup>* is admissible, that is, *<sup>u</sup>* <sup>∈</sup> *Uad*, and *y* is the state associated with *u*.

Then, let us show that the pair (*y*, *u*) satisfies the state equations given by (21), that is, for every *<sup>v</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>), we have the following convergences:

$$\begin{split} \langle \mathcal{Y}\_{\mathcal{W}'t'}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}'),L^{2}([0,T];\mathcal{V})} &\longrightarrow \langle \overline{\mathcal{Y}}\_{t'}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}),L^{2}([0,T];\mathcal{V})} \\ \langle \mathcal{A}\mathcal{y}\_{\mathcal{W}'t'}\mathcal{v} \rangle\_{L^{2}([0,T);\mathcal{V}'),L^{2}([0,T];\mathcal{V})} &\longrightarrow \langle \mathcal{A}\overline{\mathcal{y}}\_{t}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}),L^{2}([0,T];\mathcal{V})} \\ \langle \mathcal{A}\mathcal{y}\_{\mathcal{W}'t'}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}'),L^{2}([0,T];\mathcal{V})} &\longrightarrow \langle \mathcal{A}\overline{\mathcal{y}}\_{\mathcal{W}'t'}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}),L^{2}([0,T];\mathcal{V})} \\ \langle \mathcal{u}\_{\mathcal{W}'t}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}'),L^{2}([0,T];\mathcal{V})} &\longrightarrow \langle \overline{\mathcal{u}}\_{t}\mathcal{v} \rangle\_{L^{2}([0,T];\mathcal{V}'),L^{2}([0,T];\mathcal{V})} \end{split}$$

as *n* −→ ∞.

Then, according to Lemma 2, we obtain the convergence of the non-linear term from:

*B*(*yn*), *v <sup>L</sup>*2([0,*T*)];V),*L*2([0,*T*];V) −→ *B*(*y*, *v*)*<sup>L</sup>*2([0,*T*)];V),*L*2([0,*T*];V) as *n* −→ +∞.

Consequently, all the terms in the weak formulation of the state equation converge, and

$$\langle \overline{\mathcal{Y}}\_t + \mu A \overline{\mathcal{Y}} + B(\overline{\mathcal{Y}}) + M(\overline{\mathcal{Y}}) - \overline{\pi}, \upsilon \rangle\_{L^2([0,T);\mathbb{V}^\prime), L^2([0,T];\mathbb{V})} = 0,$$

is fulfilled for all *<sup>v</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; <sup>V</sup>).

Moreover, since the imbedding *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) -<sup>→</sup> *<sup>C</sup>*([0, *<sup>T</sup>*]; <sup>V</sup>) is continuous, then the mapping *<sup>w</sup>* −→ *<sup>w</sup>*(0) is linear and continuous from *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) to <sup>V</sup>, hence, we have that *yn*(0) converges weakly to *y*(0).

By the construction of the proof, we have that *y*<sup>0</sup> = *yn*(0) for all *n*, hence, it holds *y*(0) = *y*0.

Finally, it remains to show *J* = *J*(*y*, *u*) = *J*(*v*). Remember that our energy functional is given by (20), therefore we have that *J* = *J*(*v*) is a convex functional. Moreover, *J*(*v*) is continuous on *<sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>) <sup>×</sup> *<sup>L</sup>*2(Ω)4, thus by Definition 1 we have that *<sup>J</sup>* <sup>=</sup> *<sup>J</sup>*(*v*) is weakly lower semicontinuous, that is,

$$J(\mathfrak{Y}, \mathfrak{u}) \le \liminf J(y'\_{n'} \mathfrak{u}'\_n) = \mathfrak{I}.$$

Now, since (*y*, *u*) is an admissible pair, and *J* is the infimum over all admissible pairs, then it follows that *J* = *J*(*y*, *u*).

Thus, we have that (*y*, *u*) is a pair of optimal controls.

#### **5. Optimality Condition**

In this section, we will show that the optimal solution must satisfy the first-order necessary optimality condition associated with our optimal control problem given in (6).

We will study the case in which the Gâteaux derivative of the energy functional vanishes. We obtain a possible candidate solution for our optimal control, that is, if the Gâteaux derivative of our functional exists, then the optimal solution must satisfy the first-order necessary condition.

The first-order necessary condition allows conclusions to be drawn that have to do with the form and characterization of control problems.

In this part, we will establish the first-order necessary optimization condition associated with our optimal control problem given in (6). This condition will be necessary for local optimization since it is of vital importance in many aspects, that is, from the first-order necessary conditions, we can establish the candidates to be optimal controls by numerical approximations in such a way that the approximate solutions allow us to solve the first-order optimization system at a discrete level and this would be additional work that could be studied later as future research work.

Now, we can show that the optimal solution must satisfy the first-order necessary condition associated with our problem given in (20). This is performed directly using the Gâteaux derivative of our functional <sup>Ψ</sup>(*u*). In fact, for every *<sup>h</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)4) and for every *<sup>α</sup>* <sup>∈</sup> <sup>R</sup>, we have that

$$\Psi(\overline{u} + ah) \ge \Psi(\overline{u})\_\prime$$

due to the very definition of *u*. In particular, we have that

$$\forall \mathfrak{a} > 0, \frac{\Psi(\overline{\mathfrak{u}} + \mathfrak{a}h) - \Psi(\overline{\mathfrak{u}})}{\mathfrak{a}} \ge 0,$$

and

$$\forall \mathfrak{a} < 0, \; \frac{\Psi(\overline{\mathfrak{a}} + \mathfrak{a}h) - \Psi(\overline{\mathfrak{a}})}{\mathfrak{a}} \le 0,$$

which implies that the derivative at the point *<sup>α</sup>* <sup>∈</sup> <sup>R</sup> of the function *<sup>α</sup>* −→ <sup>Ψ</sup>(*<sup>u</sup>* <sup>+</sup> *<sup>α</sup>h*), is precisely the Gâteaux derivative of Ψ in the direction of *h* at the point *u* vanishes for every *<sup>h</sup>* <sup>∈</sup> *<sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)4).

Before stating our main result, let us recall the following result:

**Theorem 3** ([8])**.** *Let <sup>y</sup>*<sup>0</sup> *be in* <sup>V</sup>*; the mapping <sup>u</sup>* −→ *yu from <sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)4) *to <sup>W</sup>α*([0, *<sup>T</sup>*]; <sup>V</sup>)*, is Gâteaux differentiable* ((Υ (*u*)) · *<sup>h</sup>*1) *in every direction <sup>h</sup>*<sup>1</sup> *in <sup>L</sup>*2([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)4)*. Furthermore,* (Υ (*u*)) · *h*<sup>1</sup> = *σ*(*h*1) *is the solution of the problem given by*

$$\begin{cases} \frac{d\sigma}{dt} + \mu A \sigma + B'(y\_u) \cdot \sigma + M\sigma &= h\_1 \text{ in } Q\_{T\_\prime} \\ \sigma \in V\_\prime \\ \sigma(0) = 0 \text{ in } \Omega; \end{cases} \tag{35}$$

*we also have that <sup>σ</sup>* <sup>∈</sup> *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>V</sup>*))<sup>∩</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*;(*H*2(Ω))4) *and <sup>B</sup>* (*yu*)*σL*2(*J*1(Ω)) ) ≤ *cyuσ.*

Let us introduce the definition of locally optimal control.

**Definition 4** (locally optimal control)**.** *A control <sup>u</sup>* <sup>∈</sup> *Uad is said to be locally optimal in <sup>L</sup>*2(Ω)<sup>3</sup> *if there is a constant <sup>β</sup>* > <sup>0</sup> *such that*

$$J(\overline{y}, \overline{u}) \le J(y, u)\_{\prime}$$

*holds for all u* ∈ *Uad with u* − *uL*2(Ω)<sup>4</sup> ≤ *β. Here y and y denote the state of the system associated with u and u, respectively, that is, y* = Υ(*u*) *and y* = Υ(*u*)*.*

The first-order necessary optimization conditions are in many references, but for the related optimization conditions for optimal control problems with elliptic and parabolic partial differential equations (see [6,7,9]), they were the main references that helped us a lot, studying optimal control problems, as well as the study of stratified fluids (see [15–17]).

Now, let us show our main result of this section and show the first-order necessary optimization condition for our control problem given in (6).

**Theorem 4.** *Let <sup>U</sup> be a real Banach space, and let Uad* <sup>⊂</sup> *<sup>L</sup>*2(Ω)<sup>4</sup> *be a non-empty, convex, and closed set in <sup>L</sup>*2(Ω)<sup>4</sup> *and the functional* <sup>Ψ</sup> : *U* −→ <sup>R</sup> *be Gâteaux differentiable on Uad. Let u* ∈ *Uad be a solution of the problem*

$$\min\_{\mu \in \mathcal{U}\_{ad}} \Psi(\mu). \tag{36}$$

*Then the following optimality condition*

$$\Psi'(\overline{u})(u-\overline{u}) \gg 0\tag{37}$$

*holds for all u* ∈ *Uad. If, additionally, u* ∈ *Uad solves the variational inequality above and* Ψ *is convex, then u is the unique solution of* (36)*.*

**Proof.** Let *u* ∈ *Uad* arbitrary. Consider a convex linear combination given by:

$$u(t) = \overline{u} + t(u - \overline{u})\tag{38}$$

for any *t* ∈ [0, 1].

Now since *Uad* is non-empty, convex, and closed in *L*2(Ω)3, then we have that *u*(*t*) = *u* + *t*(*u* − *u*) ∈ *Uad* for all *t* ∈ [0, 1]. Then, from the optimality of *u*, we have that

$$\Psi(\overline{u}) \le \Psi(u(t)) \; \forall t \in [0, 1]. \tag{39}$$

Then, inserting (38) into (39) we obtain:

$$\begin{aligned} \Psi(\overline{u}) &\le \Psi(\overline{u} + t(\mu - \overline{u})) \\ &0 \le \Psi(\overline{u} + t(\mu - \overline{u})) - \Psi(\overline{u}), \end{aligned}$$

We can rewrite the last inequality by:

$$\Psi(\overline{\mathfrak{u}} + t(\mathfrak{u} - \overline{\mathfrak{u}})) - \Psi(\overline{\mathfrak{u}}) \ge 0 \,\,\forall t \in [0, 1].$$

thus, it follows that

$$\lim\_{t \to 0^+} \frac{\Psi(\overline{u} + t(u - \overline{u})) - \Psi(\overline{u})}{t} \ge 0.$$

Now, using the fact that Ψ is Gâteaux differentiable on *Uad*, and taking the limit as *t* −→ 0, we obtain

$$\Psi'(\overline{u})(u-\overline{u}) \ge 0.$$

On the other hand, let *u* ∈ *Uad* be arbitrary and let *u* ∈ *Uad* be a solution of Equation (37). Since Ψ is convex, then we have that:

$$\Psi(\boldsymbol{\mu}) - \Psi(\overline{\boldsymbol{\mu}}) \ge \Psi'(\overline{\boldsymbol{\mu}}) (\boldsymbol{\mu} - \overline{\boldsymbol{\mu}}) \,\,\forall \boldsymbol{\mu} \in \mathcal{U}\_{\text{ad}}.\tag{40}$$

In fact, for all *t* ∈ [0, 1], it follows

$$\Psi(\overline{\mu} + t(\mu - \overline{\mu})) \le (1 - t)\Psi(\overline{\mu}) + t\Psi(\mu),$$

hence,

$$\begin{aligned} \Psi(\boldsymbol{u}) - \Psi(\overline{\boldsymbol{u}}) &= \frac{(1 - t)\Psi(\overline{\boldsymbol{u}}) + t\Psi(\boldsymbol{u}) - \Psi(\overline{\boldsymbol{u}})}{t} \\ &\geq \frac{\Psi(\overline{\boldsymbol{u}} + t(\boldsymbol{u} - \overline{\boldsymbol{u}})) - \Psi(\overline{\boldsymbol{u}})}{t} \\ &\geq \Psi'(\overline{\boldsymbol{u}})(\boldsymbol{u} - \overline{\boldsymbol{u}}) \quad (t \longrightarrow 0^{+}) .\end{aligned}$$

Then, from Equations (37) and (40), we obtain that

$$\begin{aligned} \left| \Psi(u) - \Psi(\overline{u}) \right| &\geq \left| \Psi'(\overline{u}) \right| (u - \overline{u}) \\ &\geq 0 \quad \forall u \in \mathcal{U}\_{ad}. \end{aligned}$$

Therefore, we have that *u* is an optimal solution.

In order to characterize optimal solutions, we introduce the adjoint problem to the equations, which describes the non-linear motion in the dynamics of viscous and incompressible stratified fluids in R3.

**Theorem 5** (Necessary Condition)**.** *Let u be a locally optimal control for* (20) *with associated state y* = Υ(*u*)*. Then there exists a unique solution η* ∈ *V, which is the weak solution of the adjoint equation*

$$\begin{cases} -\overline{\eta}\_t + \mu A \overline{\eta} + B'(\overline{y})^\* \overline{\eta} + M \overline{\eta} = (\hat{y} - y\_d) \text{ with } \mathbf{x} \in \Omega, \ t > 0, \\ \nabla \cdot \overline{\eta} = 0, \text{ with } \mathbf{x} \in \Omega, \ t > 0, \\ \overline{\eta}(\mathbf{x}, t) = 0, \text{ with } \mathbf{x} \in \partial \Omega, \ t > 0, \\ \overline{\eta}(T) = 0 \text{ with } \mathbf{x} \in \Omega. \end{cases} \tag{41}$$

*Moreover, the following inequality*

$$\int\_{0}^{T} \int \left(\overline{\eta} + \lambda \left(\overline{u} - u\_{d}\right)\right) \cdot \left(u - \overline{u}\right) \, \mathrm{d}x \, \mathrm{d}t \geq 0, \text{ for all } u \in \mathsf{U}\_{\mathrm{ad}} \subset L^{2}(\mathbb{Q}\_{T})^{3},\tag{42}$$

*is satisfied.*

**Proof.** First of all, let us work with our reduced energy functional Ψ given in (24), which is given by:

$$J(\mathfrak{u}) := J(\mathfrak{Y}(\mathfrak{u}), \mathfrak{u}) = \mathfrak{Y}(\mathfrak{u}),$$

where we have that

$$\Psi(u) = \frac{1}{2} \int\_0^T \|\mathbf{Y}(u) - y\_d\|\_{L^2(\Omega)}^2 \, \mathrm{d}t + \frac{\lambda}{2} \int\_0^T \|u - u\_d\|\_{L^2(\Omega)}^2 \, \mathrm{d}t.\tag{43}$$

By Banach space optimization principles, we know that the variational inequality

$$\Psi'(\overline{u})(u-\overline{u}) \ge 0,\text{ for all } u \in \mathcal{U}\_{ad}$$

is a necessary condition for local optimality of *u*. It remains to compute Ψ and to derive the adjoint system. Let us write Ψ in the form given by (43).

The first derivative Ψ at *u* is characterized by

$$\begin{split} \Psi'(\overline{u})(u-\overline{u}) &= \int\_{0}^{T} \int\_{\Omega} (y-y\_{d}) \cdot \mathbf{Y}'(u)(u-\overline{u}) \, \mathrm{d}x \, \mathrm{d}t + \int\_{0}^{T} \int\_{\Omega} \lambda (u-u\_{d}) \cdot (u-\overline{u}) \, \mathrm{d}x \, \mathrm{d}t \\ &= \int\_{0}^{T} (y-y\_{d}, \sigma)\_{L^{2}(\Omega)} \, \mathrm{d}x \, \mathrm{d}t + \int\_{0}^{T} (\lambda (u-u\_{d}), (u-\overline{u}))\_{L^{2}(\Omega)} \, \mathrm{d}x \, \mathrm{d}t \text{ for all } u \in \mathcal{U}\_{\mathrm{ad}}, \end{split}$$

where *y* = Υ(*u*) and *σ* = Υ (*u*)(*u* − *u*) denote the weak solution of the equation given by (35). Let *η* be a test function. Now, multiplying by *η* in the weak formulation of (35) and integrating over Ω, we obtain

$$(\sigma, \overline{\eta}) + \mu A(\sigma, \overline{\eta}) + b(\sigma, \overline{\sigma}, \overline{\eta}) + b(\overline{\sigma}, \sigma, \overline{\eta}) + (M\sigma, \overline{\eta}) = (u - \overline{u}, \overline{\eta}),\tag{44}$$

Now, in the same way we can introduce *σ* in the weak formulation of the equation given by (41) and integrate over Ω and we obtain the following

$$(\eta\_{\prime}, \sigma) + \mu A(\eta, \sigma) + b(\eta, \overline{\eta}, \sigma) + b(\overline{\eta}, \eta, \sigma) + (M\eta, \sigma) = (\overline{y} - y\_d, \sigma). \tag{45}$$

From Equations (44) and (45), it follows that

$$\begin{aligned} (\iota - \overline{u}, \overline{\eta})\_{L^2(\Omega)} &= (\overline{y} - y\_{d\prime}\sigma)\_{L^2(\Omega)} = (\overline{y} - y\_{d\prime}Y'(\iota)(\iota - \overline{\iota}))\_{L^2(\Omega)} \\ &= (Y'(\iota)^\*(\overline{y} - y\_d), (\iota - \overline{\iota}))\_{L^2(\Omega)} \\ &= (Y'(\iota)^\*(Y(\iota) - y\_d)\_{\prime}(\iota - \overline{\iota}))\_{L^2(\Omega)}. \end{aligned}$$

Thus, we have that

$$
\overline{\eta} = (\mathbf{Y}'(\overline{\boldsymbol{\pi}})^\*(\overline{\boldsymbol{y}} - \boldsymbol{y}\_d)) = (\mathbf{Y}'(\overline{\boldsymbol{\pi}})^\*(\mathbf{Y}(\overline{\boldsymbol{\pi}}) - \boldsymbol{y}\_d)).
$$

Therefore, it follows that

$$\begin{aligned} &\int\_0^T \int \left(\overline{\eta} + \lambda \left(\overline{u} - u\_d\right)\right) \cdot \left(u - \overline{u}\right) \, \mathrm{d}x \, \mathrm{d}t \\ &= \int\_0^T \left(\overline{\eta} + \lambda \left(\overline{u} - u\_d\right), u - \overline{u}\right) \, \mathrm{d}x \, \mathrm{d}t \geq 0, \quad \text{for all } u \in \mathcal{U}\_{\mathrm{ad}}, \end{aligned}$$

where *η* is the solution for the system (41).

**Author Contributions:** Conceptualization, T.C. and H.C., methodology E.A. and T.C.,writing original draft preparation, H.C., investigation, T.C. and H.C., supervision, E.A. and T.C. The authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** There is no external funding for this research.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

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

#### **References**


## *Article* **Analysis of Solutions, Asymptotic and Exact Profiles to an Eyring–Powell Fluid Modell**

**José Luis Díaz 1,\*, Saeed Ur Rahman 2, Juan Carlos Sánchez Rodríguez 1, María Antonia Simón Rodríguez 1, Guillermo Filippone Capllonch <sup>1</sup> and Antonio Herrero Hernández <sup>1</sup>**


**Abstract:** The aim of this article was to provide analytical and numerical approaches to a onedimensional Eyring–Powell flow. First of all, the regularity, existence, and uniqueness of the solutions were explored making use of a variational weak formulation. Then, the Eyring–Powell equation was transformed into the travelling wave domain, where analytical solutions were obtained supported by the geometric perturbation theory. Such analytical solutions were validated with a numerical exercise. The main finding reported is the existence of a particular travelling wave speed *a* = 1.212 for which the analytical solution is close to the actual numerical solution with an accumulative error of <10−3.

**Keywords:** travelling waves; Eyring–Powell; geometric perturbation; nonlinear reaction–diffusion; unsteady flow

**MSC:** 35Q35; 35B65; 76D05

#### **1. Introduction**

The Eyring–Powell flow is a type of non-Newtonian fluid of paramount relevance in industrial areas, manufacturing, and biological technology. Some trivial examples of non-Newtonian fluids are given by bubbles, boiling, plastic foam processing, columns, toothpaste, mud, honey, and custard. Non-Newtonian fluids are further classified into different classes by virtue of their rheological characteristic conditions. The Eyring–Powell fluid is one such subclass of non-Newtonian fluids with particular features linked with the kinetic theory of liquids. In their seminal paper, Metzner and Otto [1] considered a non-Newtonian fluid focused on the relationship between the speed of flow and shear rate. In 1982, Rajagopal [2] considered the incompressible, unidirectional, and unsteady conditions of a second-grade fluid to obtain solutions for a flow between two rigid plates in which one suddenly starts moving. Later on, with the help of Gupta [3], they established the exact solution for the same kind of fluid between porous plates. These cited seminal works have attracted the attention of the scientific community, leading to further research paths with the same topical background in non-Newtonian fluids. Eldabe et al. [4] obtained results applicable in the field of medicine and the study of blood flow, analysing the effect of coupling forces on an unstable non-Newtonian flow of MHD between two parallel fixed porous plates under a uniform external magnetic field. Another study, carried out by Shao and Lo [5], modelled the hydrodynamics of incompressible particles (SPHs) to simulate Newtonian and non-Newtonian flows with free surfaces. The authors were able to verify the proper functioning of the model in problems such as dam breaks in 2D. Another example of outstanding interest in this regard was the study carried out by Fetecau [6].

**Citation:** Díaz, J.L.; Rahman, S.U.; Sánchez Rodríguez, J.C.; Simón Rodríguez, M.A.; Filippone Capllonch, G.; Herrero Hernández, A. Analysis of Solutions, Asymptotic and Exact Profiles to an Eyring–Powell Fluid Modell. *Mathematics* **2022**, *10*, 660. https:// doi.org/10.3390/math10040660

Academic Editor: Patricia J. Y. Wong

Received: 20 January 2022 Accepted: 17 February 2022 Published: 20 February 2022

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

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Here, solutions were established for unidirectional transient flows of non-Newtonian fluids in pipe-like domains.

Under particular rheological properties describing a non-Newtonian fluid, further applications have been accounted for by the theory of magnetohydrodynamics (MHD). Akbar [7] established the solution for a flow of a two-dimensional fluid under the effect of a magnetic field over stretching surfaces. Hina [8] analysed the heat transfer for the magnetohydrodynamic flow of the Eyring–Powell fluid. Later, Bhatti et al. [9] considered the same MHD fluid over permeable stretching surfaces. In this direction, other relevant studies can be considered (refer to [10–15]).

Further relevant topics in applied sciences involving Eyring–Powell fluids can be mentioned. In [16], the authors analysed the characteristics of the flow of Eyring–Powell nanofluids through a rotating disk subject to various physical phenomena such as a sliding flow and a magnetic field together with homogeneous and heterogeneous reactions. To this end, the proposed equations were solved by a numerical method based on the Runge– Kutta–Fehlberg method of 4th–5th order. Furthermore, in [17], the authors developed a computational technique for a three-dimensional Eyring–Powell fluid with activation energy on a stretched sheet with sliding effects. The resulting nonlinear system of PDEs was transformed into a nonlinear system of ODEs, and a shooting method was explored accordingly. The analysis in [18] discussed the flow and heat transfer of the Eyring–Powell MHD fluid in an infinite circular pipe. The explored solutions of different viscous terms were calculated numerically with the help of an iterative technique.

Note that in all the previously cited references, attention was mainly set on the numerical schemes in search of particular solutions. Analytical conceptions remain within the scope of dimensional analysis.

Further analytical approaches can be found in [19], where a homotopy approach was employed to construct solutions for a boundary layer with natural convection on a permeable vertical plate with thermal radiation. Afterwards, the differential quadrature method (DQM) was used to validate solutions for different parametrical cases involving the local Nusselt number and the local Sherwood number. In [20], the authors used the ADM-Padé approach to study analytical solutions for the deflection and pull-in instability of nanocantilever electromechanical switches, showing the remarkable accuracy compared with the numerical results. The authors claimed the possibility of extending their results to solve a wide range of instability problems. Furthermore, in [21], the authors studied a viscoelastic nanofluid with optimisation techniques subject to the proposal of a certain solution that was progressively optimised. To account for further analytical approaches, in [22], perturbation solutions were obtained for low-Reynolds–Eyring–Powell flow to obtain velocity, temperature, concentration, and stream functions.

After having cited some paramount studies involving analytical conceptions, it shall be noted that in the present study, the intention was to go deeper into the advances of the theory of PDEs to construct profiles of solutions. Unlike the previously cited studies, solutions were explored within the theory of travelling waves. Such a theory was firstly introduced by Kolmogorov, Petrovskii, and Piskunov [23], in combustion theory, and by Fisher [24], to predict the interaction of genes. The main question, introduced by the mentioned authors, was related to the search for an appropriate travelling wave speed for which the analytical travelling wave profile converges to the actual profile (solution of the actual problem, not converted into the travelling domain). Both the travelling profile and the actual one were shown to have the same exponential behaviour. This spirit was kept in our present analysis: indeed, one question to answer is related to the search for an appropriate travelling wave speed for which the analytically obtained solution converges to the actual one (obtained by numerical means) with a certain error tolerance. This was the main target of our analysis, but previously, the regularity, existence, and uniqueness of the solutions were shown. Later, the geometric perturbation theory was employed to support the construction of the analytical profiles of the solutions. These obtained profiles were validated afterwards via a numerical exercise.

#### **2. Mathematical Model**

We consider an incompressible, unsteady, and one-dimensional electrically conducting Eyring–Powell fluid. Under these assumptions, the velocity field is given by **V** = (*u*1(*y*), 0, 0), where *u*1(*y*) refers to the first velocity component. Note that the proposed problem refers to an open geometry not shaped by dedicated containers or stretched by boundary conditions. The continuity and constitutive equations for an Eyring–Powell fluid are generally given by (refer to [25,26] for an additional discussion on the Eyring–Powell governing equations):

$$\operatorname{div} \mathbf{V} = 0,\tag{1}$$

and:

$$
\rho\_f \frac{d\mathbf{V}}{dt} = \operatorname{div} \mathbf{A} + \mathbf{J} \times \mathbf{B} \tag{2}
$$

where *ρ<sup>f</sup>* refers to the density, **J** is the current density, **B** is the magnetic field, which can be split into **B** = **B**<sup>0</sup> + **b** where **B**<sup>0</sup> and **b** are the imposed and induced magnetic fields, respectively, and **A** is given by:

$$\mathbf{A} = -p\mathbf{I} + \tau\_{ij\prime} \tag{3}$$

$$
\delta \upsilon \mathbf{B} = 0, \quad \upsilon \upsilon \mathbf{l} \mathbf{B} = \mu\_1 \mathbf{j}, \quad \upsilon \upsilon \mathbf{l} \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} \tag{4}
$$

$$\mathbf{J} = \sigma(\mathbf{E} + \mathbf{V} \times \mathbf{B}),\tag{5}$$

where *p* is the pressure field, **I** is the identity tensor, *μ*<sup>1</sup> is the magnetic permeability, **E** is the electric field, *σ* is the electric conductivity, and *τij* is the shear stress tensor of an Eyring–Powell fluid [11,13] given by:

$$
\pi\_{ij} = \mu \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{1}{\beta} \sinh^{-1} \left( \frac{1}{d\_1} \frac{\partial u\_i}{\partial \mathbf{x}\_j} \right), \tag{6}
$$

where *μ* is the dynamic viscosity and *β* and *d*<sup>1</sup> are characteristic constants of the Powell-Eyring model. Consider that sinh−<sup>1</sup> <sup>1</sup> *d*1 *∂ui ∂xj* <sup>∼</sup><sup>=</sup> <sup>1</sup> *d*1 *∂ui <sup>∂</sup>xj* <sup>−</sup> <sup>1</sup> 6 <sup>1</sup> *d*1 *∂ui ∂xj* 3 , 1 *d*1 *∂ui ∂xj* <sup>≤</sup> 1. The governing equation, in the absence of an induced magnetic field, can be written as:

$$\frac{\partial u\_1}{\partial t} = -\frac{1}{\rho} \frac{dP}{dx} + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) \frac{\partial^2 u\_1}{\partial y^2} - \frac{1}{2\beta d\_1^3 \rho\_f} \left(\frac{\partial u\_1}{\partial y}\right)^2 \frac{\partial^2 u\_1}{\partial y^2} - \frac{\sigma B\_0^2 u\_1}{\rho\_f}.\tag{7}$$

where *v* = *<sup>μ</sup> <sup>ρ</sup><sup>f</sup>* is the kinematic viscosity. After differentiation in (7) with *x*:

$$-\frac{1}{\rho} \frac{d^2 P}{d\mathfrak{x}^2} = 0, \quad -\frac{1}{\rho} \frac{dP}{d\mathfrak{x}} = A\_1.$$

Using the value of <sup>−</sup><sup>1</sup> *ρ dP dx* in (7), we obtain:

$$\frac{\partial u\_1}{\partial t} = A\_1 + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) \frac{\partial^2 u\_1}{\partial y^2} - \frac{1}{2\beta d\_1^3 \rho\_f} \left(\frac{\partial u\_1}{\partial y}\right)^2 \frac{\partial^2 u\_1}{\partial y^2} - \frac{\sigma B\_0^2 u\_1}{\rho\_f}.\tag{8}$$

with the following initial condition:

$$
u\_1(y,0) = 
u\_0(y) \in L^1\_{loc}(R) \cap L^\infty(R). \tag{9}$$

#### **3. Preliminaries**

The proposed Eyring–Powell model in (8) is expressed making use of a weak formulation to support the analysis of the regularity, existence, and uniqueness of the solutions.

**Definition 1.** *Consider a test function <sup>φ</sup>*<sup>2</sup> <sup>∈</sup> *<sup>C</sup>*∞(*R*) *defined in* (0, *<sup>T</sup>*)*, such that for* <sup>0</sup> <sup>&</sup>lt; *<sup>τ</sup>* <sup>&</sup>lt; *<sup>t</sup>* <sup>&</sup>lt; *T, the following weak formulation of* (8) *holds:*

$$\begin{split} \int\_{R} u\_{1}(t)\phi\_{2}(t)dy &= \int\_{R} u\_{1}(\tau)\phi\_{2}(\tau)dy + \int\_{\tau}^{t} \int\_{R} u\_{1} \frac{\partial \phi\_{2}}{\partial s} dyds \\ &+ A\_{1} \int\_{\tau}^{t} \int\_{R} \phi\_{2} dyds + \left(v + \frac{1}{\beta d\_{1}\rho\_{f}}\right) \int\_{\tau} \int\_{R} u\_{1} \frac{\partial^{2}\phi\_{1}}{\partial y^{2}} dyds \\ &+ \frac{1}{6\beta d\_{1}^{3}\rho\_{f}} \int\_{\tau} \int\_{R} \left(\frac{\partial u\_{1}}{\partial y}\right)^{3} \frac{\partial \phi\_{2}}{\partial y} dyds - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \int\_{\tau} \int\_{R} u\_{1} \phi\_{2} dyds. \end{split}$$

In addition, the following definition holds:

**Definition 2.** *Given a finite spatial location r*0*, admit a ball Br centred in r*<sup>0</sup> *and with radiusr r*0*. In the proximity of the borders <sup>∂</sup>Br and for* <sup>0</sup> <sup>&</sup>lt; *<sup>s</sup>* <sup>&</sup>lt; *<sup>τ</sup>* <sup>&</sup>lt; *<sup>t</sup>* <sup>&</sup>lt; *T, the following equation is defined:*

$$u\_1 \frac{\partial \phi\_2}{\partial s} + A\_1 \phi\_2 + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) u\_1 \frac{\partial^2 \phi\_2}{\partial y^2} + \frac{1}{6 \beta d\_1^3 \rho\_f} \left(\frac{\partial u}{\partial y}\right)^3 \frac{\partial \phi\_2}{\partial y} - \frac{\sigma B\_0^2}{\rho\_f} u\_1 \phi\_2 = 0,\tag{10}$$

*in Br* × (0, *T*), *with the following boundary and initial conditions:*

$$0 < \frac{\partial \phi\_2}{\partial y} = \phi\_2 \ll 1,$$

*and:*

$$
\mu\_1(y,0) = \mu\_0(y) \in L^1\_{loc}(\mathbb{R}) \cap L^\infty(\mathbb{R}).
$$

#### **4. Existence and Uniqueness Analysis**

The following theorem aims to show the existence and bounds of the solutions:

**Theorem 1.** *Given <sup>u</sup>*0(*y*) <sup>∈</sup> *<sup>L</sup>*<sup>1</sup> *loc*(*R*) <sup>∩</sup> *<sup>L</sup>*∞(*R*)*, then the solution is bounded for all* (*y*, *<sup>t</sup>*) <sup>∈</sup> *Br* × [*τ*, *T*) *with r* 1.

**Proof.** Consider a certain value *<sup>η</sup>* <sup>∈</sup> *<sup>R</sup>*<sup>+</sup> such that the following cut-off function is defined (see [27,28]):

$$\begin{aligned} \psi\_{\eta} &\in \mathbb{C}\_{0}^{\infty}(\underline{y}, t), & 0 \le \psi\_{\eta} \le 1, \\ \psi\_{\eta} &= 1 \text{ in } B\_{r - \eta \prime} \quad \psi\_{\eta} = 0 \text{ in } R - B\_{r - \eta \prime} \end{aligned}$$

so that:

$$\left| \frac{\partial \psi\_{\eta}}{\partial \eta} \right| = \frac{B\_{a}}{\eta} \lambda$$

where *Ba* is a suitable constant. Multiplying (10) by *ψη* and integrating in *Br* × [*τ*, *T*), we obtain:

$$\int\_{\tau}^{t} \int\_{B\_{r}} u\_{1} \frac{\partial \phi\_{2}}{\partial s} \psi\_{\eta} \, dyds + A\_{1} \int\_{\tau}^{t} \int\_{B\_{r}} \phi\_{2} \psi\_{\eta} \, dyds + \left(v + \frac{1}{\beta d\_{1} \rho\_{f}}\right) \int\_{\tau}^{t} \int\_{B\_{r}} u\_{1} \frac{\partial^{2} \phi\_{2}}{\partial y^{2}} \psi\_{\eta} \, dyds$$

$$+ \frac{1}{\beta \beta d\_{1}^{3} \rho\_{f}} \int\_{\tau}^{t} \int\_{B\_{r}} \left(\frac{\partial u\_{1}}{\partial y}\right)^{3} \frac{\partial \phi\_{2}}{\partial y} \psi\_{\eta} \, dyds - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \int\_{\tau}^{t} \int\_{B\_{r}} u\_{1} \phi\_{2} \psi\_{\eta} \, dyds = 0. \tag{11}$$

Now, admit an arbitrary *<sup>m</sup>* <sup>&</sup>gt; 1 and some large *<sup>r</sup>*<sup>0</sup> <sup>&</sup>gt; 1 [27,28]:

$$\int\_{\tau} u\_1 ds \le \int\_{\tau} u\_1^m ds \le D\_1(\tau) r^{\frac{2m}{m-1}}.$$

Considering the spatial variable *y* close to *∂Br*, it can be assumed that *y* ∼ *r*. Then, for *m* = 2, it holds that:

$$\int\_{\tau}^{t} \mu\_1 ds \le D\_1(\tau)r^4, \quad \int\_{\tau}^{t} \left(\frac{\partial \mu\_1}{\partial y}\right)^3 ds \le 64 \, D\_1^3(\tau)r^9.$$

The integral for the diffusion term reads:

$$\begin{split} & \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) \int\_{\tau}^{t} \int\_{\mathcal{B}\_r} u\_1 \frac{\partial^2 \phi\_2}{\partial y^2} \, \psi\_\eta dyds \\ & \qquad \le \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) \int\_{\mathcal{B}\_r} D\_1(\tau) r^2 \frac{\partial^2 \phi\_2}{\partial y^2} \, \psi\_\eta dy \\ & \qquad = \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) D\_1(\tau) r^2 \left( \left( \frac{\partial \phi\_2}{\partial y} \psi\_\eta \right)\_{\partial \mathcal{B}\_r} - \int\_{\mathcal{B}\_r} \frac{\partial \phi\_2}{\partial y} \, \frac{\partial \psi\_\eta}{\partial y} dy \right). \end{split}$$

As *<sup>r</sup>* 1 and taking *<sup>φ</sup>*<sup>2</sup> sufficiently small such that *∂φ*<sup>2</sup> *<sup>∂</sup><sup>y</sup> ψη* 1 over *∂Br*, the following holds:

$$\begin{split} & \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) \Big( \int\_{\tau} \mu\_1 \frac{\partial^2 \phi\_2}{\partial y^2} \psi\_{\eta} dy ds \\ & \qquad = - \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) \int\_{B\_r} D\_1(\tau) r^2 \frac{\partial \phi\_2}{\partial y} \frac{\partial \psi\_{\eta}}{\partial y} dy \\ & \qquad \le \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) D\_1(\tau) \int\_{B\_r} r^2 \frac{\partial \phi\_2}{\partial y} \frac{B\_\alpha}{\eta} dy \\ & \qquad = \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) B\_\alpha D\_1(\tau) \int\_{B\_r} r \frac{\partial \phi\_2}{\partial y} dy \, \end{split}$$

and:

$$\frac{1}{6\beta d\_1^3 \rho\_f} \int\_{\tau} \int\_{B\_r} \left(\frac{\partial u\_1}{\partial y}\right)^3 \frac{\partial \phi\_2}{\partial y} \psi\_\eta \,dyds \le \frac{32}{3\beta d\_1^3 \rho\_f} \int\_{B\_r} D\_1^3(\tau) r^9 \frac{\partial \phi\_2}{\partial y} \psi\_\eta \,dyds$$

Now:

$$\begin{split} \frac{1}{6\beta d\_1^3 \rho\_f} \int\_{\mathbf{r}} \int\_{\mathbf{B}\_r} \left(\frac{\partial u\_1}{\partial y}\right)^3 \frac{\partial \rho\_2}{\partial y} \psi\_\eta dyds & \leq \ -\frac{32}{3\beta d\_1^3 \rho\_f} \int\_{\mathbf{B}\_r} D\_1^3(\mathbf{r}) r^9 \phi\_2 \frac{\partial \psi\_\eta}{\partial y} dy \\ & \leq \ \ \frac{32}{3\beta d\_1^3 \rho\_f} \int\_{\mathbf{B}\_r} D\_1^3(\mathbf{r}) r^9 \phi\_2 \frac{\mathbf{B}\_a}{\eta} dy \\ & = \ \ \ \frac{32 D\_1^3(\mathbf{r})}{3\beta d\_1^3 \rho\_f} \int\_{\mathbf{B}\_r} r^8 \phi\_2 dy. \end{split} \tag{12}$$

Using the expressions (12) and (12) in (11), the following holds:

$$
\int\_{\tau}^{t} \int\_{B\_{r}} \eta\_{1} \frac{\partial \phi\_{2}}{\partial s} \psi\_{\eta} dyds + A\_{1} \int\_{\tau}^{t} \int\_{B\_{r}} \phi\_{2} \psi\_{\eta} dyds \leq \left(\upsilon + \frac{1}{\beta d\_{1} \rho\_{f}}\right) B\_{d} D\_{1}(\tau) \int\_{B\_{r}} r \frac{\partial \phi\_{2}}{\partial y} dy
$$

$$
+ \frac{32 B\_{d} D\_{1}^{3}(\tau)}{3 \beta d\_{1}^{3} \rho\_{f}} \int\_{B\_{r}} r^{8} \phi\_{2} dy + \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \int\_{\tau}^{t} \int\_{B\_{r}} u\_{1} \phi\_{2} \psi\_{\eta} dyds. \tag{13}
$$

Next, consider a test function *φ*<sup>2</sup> of the form:

$$
\phi\_2(r,s) = e^{-ks} \left(1 + r^2\right)^{-a}.\tag{14}
$$

We can choose *a* in such a way that (13) is convergent; therefore:

$$\begin{split} & \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) B\_2 B\_1(\tau) \int\_{\mathcal{B}\_r} r \frac{\partial \phi\_2}{\partial y} dy + \\ & \frac{32 B\_d D\_1^3(\tau)}{3 \beta d\_1^3 \rho\_f} \int\_{\mathcal{B}\_r} r^8 \phi\_2 dy + \frac{\sigma B\_0^3 B\_1(\tau)}{\rho\_f} \int\_{\mathcal{B}\_r} r^2 \phi\_2 \psi\_{\mathcal{V}} dy \\ & \le 2a \left( v + \frac{1}{\beta d\_1 \rho\_f} \right) B\_d B\_1(\tau) \int\_{\mathcal{B}\_r} e^{-ks} r^{-2a} dr \\ & \quad + \frac{32 B\_d D\_1^3(\tau)}{3 \beta d\_1^3 \rho\_f} \int\_{\mathcal{B}\_r} r^{8-2a} \phi\_2 dr + \frac{\sigma B\_0^2 B\_1(\tau)}{\rho\_f} \int\_{\mathcal{B}\_r} e^{-ks} r^{2-2a} dr. \end{split} \tag{15}$$

For *<sup>a</sup>* <sup>&</sup>gt; 4 and *<sup>r</sup>* <sup>→</sup> <sup>∞</sup>, the following holds:

$$\int \left(v + \frac{1}{\beta d\_1 \rho\_f} \right) B\_d B\_1(\tau) \int\_{B\_r} r \frac{\partial \phi\_2}{\partial y} dy + \frac{\sigma B\_0^2 B\_1(\tau)}{\rho\_f} \int\_{B\_r} r^2 \phi\_2 \psi\_\eta dy \le 0. \tag{16}$$

Putting (16) into (13):

$$\int\_{\tau}^{t} \int\_{B\_{r}} u\_{1} \frac{\partial \phi\_{2}}{\partial s} \psi\_{\eta} dyds + A\_{1} \int\_{\tau}^{t} \int\_{B\_{r}} \phi\_{2} \psi\_{\eta} dyds \leq 0. \tag{17}$$

As both integrals are finite in *<sup>τ</sup>* < *<sup>s</sup>* < *<sup>t</sup>* < *<sup>T</sup>*, it is possible to conclude the theorem principles related to the bound of the solutions in *R* × (0, *T*).

The next intention is to show the boundness of *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* .

**Theorem 2.** *Given u*1(*y*) *as the solution of* (8)*, then <sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>y</sup> is bounded for*(*y*, *t*) ∈ *R* × (0, *T*)*.*

**Proof.** Multiplying the equation (8) by *u*<sup>1</sup> and using integration by parts:

$$\begin{aligned} \frac{d}{dt} \int\_{\mathcal{R}} |u\_1|^2 dy &= \quad A\_1 \int\_{\mathcal{R}} u\_1 dy - \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) \int\_{\mathcal{R}} \left(\frac{\partial u\_1}{\partial y}\right)^2 dy \\ &+ \frac{1}{6 \beta d\_1^3 \rho\_f} \int\_{\mathcal{R}} \left(\frac{\partial u\_1}{\partial y}\right)^4 dy - \frac{\sigma B\_0^2}{\rho\_f} \int\_{\mathcal{R}} |u\_1|^2 dy \end{aligned}$$

which implies that:

$$\int\_{R} \left(\frac{\partial u\_1}{\partial y}\right)^2 \left(\frac{1}{6\beta d\_1^3 \rho\_f} \left(\frac{\partial u\_1}{\partial y}\right)^2 - \left(v + \frac{1}{\beta d\_1 \rho\_f}\right)\right) dy = \frac{d}{dt} \int\_{R} |u\_1|^2 dy$$

$$-A\_1 \int\_{R} u\_1 dy - \frac{\sigma B\_0^2}{\rho\_f} \int\_{R} |u\_1|^2 dy.$$

After integration on both sides:

$$\int\_{0}^{t} \int\_{R} \left(\frac{\partial u\_{1}}{\partial y}\right)^{2} \left(\frac{1}{6\beta d\_{1}^{3}\rho\_{f}} \left(\frac{\partial u\_{1}}{\partial y}\right)^{2} - \left(v + \frac{1}{\beta d\_{1}\rho\_{f}}\right)\right) dy ds = \int\_{R} |u\_{1}(y,t)|^{2} dy$$

$$-\int\_{R} |u\_{0}(y)|^{2} dy - A\_{1} \int\_{0}^{t} \int\_{R} u\_{1} dyds - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \int\_{0}^{t} \int\_{R} |u\_{1}|^{2} dy ds. \tag{18}$$

From Theorem (1), the right-hand side of (18) is bounded; therefore, we can choose *A*<sup>2</sup> such that:

$$\int\_{0}^{t} \int\_{R} \left(\frac{\partial u\_{1}}{\partial y}\right)^{2} \left(\frac{1}{6\beta d\_{1}^{3}\rho\_{f}} \left(\frac{\partial u\_{1}}{\partial y}\right)^{2} - \left(v + \frac{1}{\beta d\_{1}\rho\_{f}}\right)\right) dy ds \le A\_{2}.\tag{19}$$

which permits concluding that *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* is bounded in *R* × (0, *t*) where we can admit *t* = *T*.

The next intention is to show the uniqueness of the solution.

**Theorem 3.** *Let us admit <sup>u</sup>*<sup>1</sup> <sup>&</sup>gt; <sup>0</sup> *as a minimal solution and <sup>u</sup>*ˆ1 *as a maximal solution for* (8) *in R* × (0, *T*)*, then u*<sup>1</sup> *coincides with the maximal solution u*ˆ1*, i.e., the solution is unique.*

**Proof.** Consider *u*ˆ1 to be the maximal solution of (8) in *R* × (0, *T*) given by:

$$
\mu\_1(y,0) = \mu\_0(y) + \epsilon,\tag{20}
$$

with <sup>&</sup>gt; 0 arbitrarily small. In addition, let us define the minimal solution:

$$
u\_1(y,0) = 
u\_0(y).$$

The maximal and minimal solutions satisfy the following equations:

$$\frac{\partial \hat{u}\_1}{\partial t} = \quad A\_1 + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) \frac{\partial^2 \hat{u}\_1}{\partial y^2} - \frac{1}{2\beta d\_1^3 \rho\_f} \left(\frac{\partial \hat{u}\_1}{\partial y}\right)^2 \frac{\partial^2 \hat{u}\_1}{\partial y^2} - \frac{\sigma B\_0^2 \hat{u}\_1}{\rho\_f},\tag{21}$$

$$\frac{\partial u\_1}{\partial t} = -A\_1 + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) \frac{\partial^2 u\_1}{\partial y^2} - \frac{1}{2\beta d\_1^3 \rho\_f} \left(\frac{\partial u\_1}{\partial y}\right)^2 \frac{\partial^2 u\_1}{\partial y^2} - \frac{\sigma B\_0^2 u\_1}{\rho\_f}.\tag{22}$$

For every test function *<sup>φ</sup>*<sup>2</sup> <sup>∈</sup> *<sup>C</sup>*∞(*R*) and upon subtraction, the following expressions hold:

0 ≤ *R* (*u*ˆ1 − *u*1)*φ*2(*t*)*dy* = *t* 0 *R* (*u*ˆ1 − *u*1) *∂φ*<sup>2</sup> *<sup>∂</sup><sup>s</sup> dyds* + *v* + 1 *βd*1*ρ<sup>f</sup> t* 0 *R* (*u*ˆ1 − *u*1) *∂*2*φ*<sup>2</sup> *<sup>∂</sup>y*<sup>2</sup> *dyds* + 1 6*βd*<sup>3</sup> 1*ρf t* 0 *R ∂u*ˆ1 *∂y* 3 − *∂u*<sup>1</sup> *∂y* 3 *∂*2*φ*<sup>2</sup> *<sup>∂</sup>y*<sup>2</sup> *dyds* <sup>−</sup> *<sup>σ</sup>B*<sup>2</sup> 0 *ρf t* 0 *R* (*u*ˆ1 − *u*1)*φdyds* ≤ *t* 0 *R* (*u*ˆ1 − *u*1) *∂φ*<sup>2</sup> *<sup>∂</sup><sup>s</sup> dyds* <sup>+</sup> *v* + 1 *βd*1*ρ<sup>f</sup> t* 0 *R* (*u*ˆ1 − *u*1) *∂*2*φ*<sup>2</sup> *<sup>∂</sup>y*<sup>2</sup> *dyds* + 1 6*βd*<sup>3</sup> 1*ρf t* 0 *R ∂u*ˆ1 *<sup>∂</sup><sup>y</sup>* <sup>−</sup> *<sup>∂</sup>u*<sup>1</sup> *∂y ∂u*ˆ1 *∂y* 2 + *∂u*ˆ1 *∂y ∂u*<sup>1</sup> *∂y* + *∂u*<sup>1</sup> *∂y* 2 *∂φ*<sup>2</sup> *<sup>∂</sup><sup>y</sup> dyds* <sup>−</sup> *<sup>σ</sup>B*<sup>2</sup> 0 *ρf t* 0 *R* (*u*ˆ1 − *u*1)*φdyds* (23)

Based on Theorem 2's results, we can choose *<sup>A</sup>*<sup>3</sup> such that *<sup>A</sup>*<sup>3</sup> <sup>=</sup> sup{ *<sup>∂</sup>u*ˆ1 *<sup>∂</sup><sup>y</sup>* , *<sup>∂</sup>u*<sup>1</sup> *<sup>∂</sup><sup>y</sup>* }, so that the following holds:

$$\int\_{R} (\hat{\mu}\_{1} - u\_{1}) \phi\_{2}(t) dy \le \int\_{0}^{t} \int\_{R} (\hat{\mu}\_{1} - u\_{1}) \frac{\partial \phi\_{2}}{\partial s} dy ds + \left(v + \frac{1}{\beta d\_{1} \rho\_{f}}\right) \int\_{0}^{t} \int\_{R} (\hat{\mu}\_{1} - u\_{1}) \frac{\partial^{2} \phi\_{2}}{\partial y^{2}} dy ds$$

$$+ \frac{A\_{3}}{6 \beta d\_{1}^{3} \rho\_{f}} \int\_{0}^{t} \int\_{R} \left(\frac{\partial \mathbf{\hat{a}}\_{1}}{\partial y} - \frac{\partial \mathbf{u}\_{1}}{\partial y}\right) \frac{\partial \phi\_{2}}{\partial y} dy ds - \frac{\sigma B\_{0}^{3}}{\rho\_{f}} \int\_{0}^{t} \int\_{R} (\mathbf{\hat{a}}\_{1} - u\_{1}) \phi dyds$$

$$= \int\_{0}^{t} \int\_{R} (\hat{\mu}\_{1}^{\*} - u\_{1}) \frac{\partial \phi\_{2}}{\partial s} dy ds + \left(v + \frac{1}{\beta d\_{1} \rho\_{f}} - \frac{A\_{3}}{6 \beta d\_{1}^{3} \rho\_{f}}\right) \int\_{0}^{t} \int\_{R} (\hat{\mu}\_{1}^{\*} - u\_{1}) \frac{\partial^{2} \phi\_{2}}{\partial y^{2}} dy ds\tag{24}$$

$$- \frac{\sigma B\_{0}^{3}}{\rho\_{f}} \int\_{0}^{t} \int\_{R} (\mathbf{\hat{a}}\_{1} - u\_{1}) \phi dyds.\tag{25}$$

Now, consider the test function given by:

$$\phi\_2(|y|,s) = e^{A\_4(T-s)} \left(1 + |y|^2\right)^{-b} \,\_{\prime} \tag{26}$$

where *A*<sup>4</sup> and *b* are constants. Making the differentiation of *φ*<sup>2</sup> with regards to *s* and *y*, the following holds:

$$\frac{\partial \phi\_2}{\partial s} = -A\_4 \phi\_2(|y|, s)\_{\prime} \quad \frac{\partial^2 \phi\_2}{\partial y^2} \le A\_5(b) \phi\_2(|y|, s)\_{\prime}$$

then:

=

−*A*<sup>4</sup> +

*v* +

$$(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1})\frac{\partial \boldsymbol{\phi}\_{2}}{\partial \mathbf{s}} + \left(\boldsymbol{\upsilon} + \frac{1}{\beta \boldsymbol{d}\_{1} \rho\_{f}} - \frac{A\_{3}}{6\beta \boldsymbol{d}\_{1}^{3} \rho\_{f}}\right)(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1})\frac{\partial^{2} \boldsymbol{\phi}\_{2}}{\partial \boldsymbol{y}^{2}} - \frac{\sigma B\_{0}^{2}}{\rho\_{f}}(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1})\boldsymbol{\phi}\_{2}$$

$$\leq \quad -A\_{4}\boldsymbol{\phi}\_{2}(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1}) + \left(\boldsymbol{\upsilon} + \frac{1}{\beta \boldsymbol{d}\_{1} \rho\_{f}} - \frac{A\_{3}}{6\beta \boldsymbol{d}\_{1}^{3} \rho\_{f}}\right)A\_{5}(\boldsymbol{b})\boldsymbol{\phi}\_{2}(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1}) - \frac{\sigma B\_{0}^{2}}{\rho\_{f}}(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1})\boldsymbol{\phi}\_{2}$$

$$= \quad \left(-A\_{4} + \left(\boldsymbol{\upsilon} + \frac{1}{6\beta \boldsymbol{d}\_{1} \rho\_{f}} - \frac{A\_{3}}{6\beta \boldsymbol{d}\_{1}^{3} \rho\_{f}}\right)A\_{5}(\boldsymbol{b}) - \frac{\sigma B\_{0}^{2}}{\rho\_{f}}\right)(\boldsymbol{\vartheta}\_{1} - \boldsymbol{\mu}\_{1})\boldsymbol{\phi}\_{2}. \tag{27}$$

0 *ρf*

(*u*ˆ1 − *u*1)*φ*2. (27)

*βd*1*ρ<sup>f</sup>*

6*βd*<sup>3</sup> 1*ρf* Using (27) in (24), we obtain:

$$\begin{split} \int\_{\mathbb{R}} (\mathfrak{a}\_{1} - u\_{1}) \mathfrak{\phi}\_{2}(t) dy &\leq \left( -A\_{4} + \left( v + \frac{1}{\beta d\_{1} \rho\_{f}} - \frac{A\_{3}}{6 \beta d\_{1}^{3} \rho\_{f}} \right) A\_{5}(b) - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \right) \\ &\times \int\_{0}^{t} \int\_{\mathbb{R}} (\mathfrak{A}\_{1} - u\_{1}) \mathfrak{\phi}\_{2} dy ds \leq \left| -A\_{4} + \left( v + \frac{1}{\beta d\_{1} \rho\_{f}} - \frac{A\_{3}}{6 \beta d\_{1}^{3} \rho\_{f}} \right) A\_{5}(b) - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \right| \\ &\times \int\_{0}^{t} \int\_{\mathbb{R}} (\mathfrak{A}\_{1} - u\_{1}) \mathfrak{\phi}\_{2} dy ds. \end{split} \tag{28}$$

Making the differentiation with regard to *t*:

$$\begin{split} \frac{d}{dt} \int\_{R} (\mathfrak{d}\_{1} - u\_{1}) \mathfrak{d}\_{2}(t) dy &\leq \left| -A\_{4} + \left( v + \frac{1}{\beta d\_{1} \rho\_{f}} - \frac{A\_{3}}{6 \beta d\_{1}^{3} \rho\_{f}} \right) A\_{5}(b) - \frac{\sigma B\_{0}^{2}}{\rho\_{f}} \right| \\ &\times \int\_{R} (\mathfrak{d}\_{1} - u\_{1}) \mathfrak{d}\_{2}(t) dy. \end{split} \tag{29}$$

Now, let us define:

$$h(t) = \int\_{\mathbb{R}} (\mathfrak{A}\_1 - \mathfrak{u}\_1) \mathfrak{d}\_2(t) d\mathfrak{y}.\tag{30}$$

Putting (30) into (29), the following holds:

$$\frac{d\mathbf{h}}{dt} \le \left| -A\_4 + \left( v + \frac{1}{\beta d\_1 \rho\_f} - \frac{A\_3}{6\beta d\_1^3 \rho\_f} \right) A\_5(b) - \frac{\sigma B\_0^2}{\rho\_f} \right| h(t), \tag{31}$$

with:

$$h(0) = \epsilon \to 0.$$

After solving (31) by standard means, we obtain *h*(*t*) = 0, i.e., *u*ˆ1 = *u*1, which shows the uniqueness of the solutions, as was intended to be proven.

#### **5. Travelling Waves' Existence and Regularity**

The travelling wave profiles are described as *u*1(*y*, *t*) = *k*(*ζ*), where *ζ* = *y* − *at* ∈ *R*, *a* refers to the travelling wave speed and *<sup>k</sup>* : *<sup>R</sup>* <sup>→</sup> (0, <sup>∞</sup>) belongs to *<sup>L</sup>*∞(*R*).

The equation (8) is transformed into the travelling wave domain as follows:

$$-ak'(\zeta) = A\_1 + \left(v + \frac{1}{\beta d\_1 \rho\_f}\right) k''(\zeta) - \frac{1}{2\beta d\_1^3 \rho\_f} \left(k'(\zeta)\right)^2 k''(\zeta) - \frac{\sigma B\_0^2}{\rho\_f} k(\zeta). \tag{32}$$

with *k* (*ζ*) <sup>&</sup>lt; 0 in the hypothesis of a purely decreasing travelling wave (this assumption is further discussed later). Now, let us consider the following new variables:

$$X = k(\mathbb{Z}), \quad \mathcal{Y} = k'(\mathbb{Z}), \tag{33}$$

such that the following system holds:

$$\begin{array}{rcl} X' &=& \Upsilon, \\ Y' &=& \frac{2\beta d\_1^3 \rho\_f}{2v\beta d\_1^3 \rho\_f + 2d\_1^2 - Y^2} \left( -aY - A\_1 + \frac{\sigma B\_0^2}{\rho\_f} X \right). \end{array} \tag{34}$$

To analyse the suggested system in the proximity of the critical point, admit *X* = 0 and *Y* = 0, yielding:

$$\mathbf{X} = \frac{A\_1 \rho\_f}{\sigma B\_0^2}.$$

Therefore, *<sup>A</sup>*1*ρ<sup>f</sup> σB*<sup>2</sup> 0 , 0 represents the system critical point.

Our intention in the coming sections was to make use of the geometric perturbation theory to characterise the existing critical point and to explore solution orbits close to such a critical point.

#### *5.1. Geometric Perturbation Theory*

In this section, we use the singular geometric perturbation theory to show the asymptotic behaviour of an appropriately defined manifold close to the critical point. Afterwards, the obtained results are used to derive a dedicated travelling wave profile.

For this purpose, admit the following manifold as:

$$N\_0 = \left\{ X, Y \;/\ X' = Y; \; Y' = \frac{2\beta d\_1^3 \rho\_f}{2v\beta d\_1^3 \rho\_f + 2d\_1^2 - Y^2} \left( -aY - A\_1 + \frac{\sigma B\_0^2}{\rho\_f} X \right) \right\},\tag{35}$$

with critical point *<sup>A</sup>*1*ρ<sup>f</sup> σB*<sup>2</sup> 0 , 0 . The perturbed manifold *N* close to *N*<sup>0</sup> in the critical point *<sup>A</sup>*1*ρ<sup>f</sup> σB*<sup>2</sup> 0 , 0 is defined as:

$$N\_{\varepsilon} = \left\{ X, Y \;/\ X' = \varepsilon Y; \; Y' = F\varepsilon \left( X - \frac{A\_1 \rho\_f}{\sigma B\_0^2} \right) \right\},\tag{36}$$

where denotes a perturbation parameter close to equilibrium (*X*1, 0) and *F* is a suitable constant, which is found after root factorisation. Firstly, admit *<sup>X</sup>*<sup>3</sup> <sup>=</sup> *<sup>X</sup>* <sup>−</sup> *<sup>A</sup>*1*ρ<sup>f</sup> σB*<sup>2</sup> 0 . Our intention was to apply the Fenichel invariant manifold theorem [29] as formulated in [30]. For this purpose, we have to show that *N*<sup>0</sup> is a normally hyperbolic manifold, i.e., the eigenvalues of *N*<sup>0</sup> in the linearised frame close to the critical point, and transversal to the tangent space, have non-zero real part. This is shown based on the following equivalent flow associated with *<sup>N</sup>*<sup>0</sup> : *X*

$$
\begin{pmatrix} X'\_3 \\ Y' \end{pmatrix} = \begin{pmatrix} 0 & \epsilon \\ F\epsilon & 0 \end{pmatrix} \begin{pmatrix} X\_3 \\ Y \end{pmatrix}.
$$

The associated eigenvalues are both real ± <sup>√</sup>*F* , which shows that *N*<sup>0</sup> is a hyperbolic manifold. Now, we want to show that the manifold *N* is locally invariant under the flow (34), so that the manifold *N*<sup>0</sup> can be shown as an asymptotic approach to *N* and vice versa. On this basis, we consider the functions:

$$\begin{array}{rcl} \psi\_1 &=& \epsilon Y\_{\prime} \\ \psi\_2 &=& FeX\_{3\prime} \end{array}$$

which are *C<sup>i</sup>* (*<sup>R</sup>* <sup>×</sup> [0, *<sup>δ</sup>*)), *<sup>i</sup>* > 0, in the proximity of the critical point *<sup>A</sup>*1*ρ<sup>f</sup> σB*<sup>2</sup> 0 , 0 . In this case, *δ* is determined based on the following flows that are considered to be measurable a.e. in *R*:

$$\left\| \left| \psi\_1^{M\_0} - \psi\_1^{M\_\epsilon} \right| \right\| \le F\epsilon \left\| X\_3 \right\| \le \delta \epsilon.$$

Since the solutions are bounded, we conclude that *δ* = *FX*3 is finite; therefore, the distance between the manifolds holds the normal hyperbolic condition for *δ* ∈ (0, ∞) and sufficiently small close to the critical point *<sup>A</sup>*1*<sup>ρ</sup> <sup>f</sup> σB*<sup>2</sup> 0 , 0 .

#### *5.2. Travelling Waves' Profiles*

Based on the normal hyperbolic condition shown for the manifold *N*<sup>0</sup> under the flow (34), asymptotic TW profiles can be obtained. For this purpose, let us consider firstly (34) such that the following family of trajectories in the phase plane (*X*,*Y*) holds:

$$\frac{dY}{dX} = \frac{2\beta d\_1^3 \rho\_f}{\left(2v\beta d\_1^3 \rho\_f + 2d\_1^2 - Y^2\right)Y} \left(-aY - A\_1 + \frac{\sigma B\_0^2}{\rho\_f}X\right) = H(X,Y). \tag{37}$$

As *H*(*X*,*Y*) is continuous and is changing the sign character if we take *X* sufficiently large and sufficiently small, it is possible to conclude the existence of a critical trajectory of the form:

$$-aX' - A\_1 + \frac{\sigma B\_0^2}{\rho\_f}X = 0\_\prime$$

which implies that:

$$X' = \frac{\sigma B\_0^2}{a\rho\_f} \left( X - \frac{A\_1 \rho\_f}{\sigma B\_0^2} \right). \tag{38}$$

.

Solving (38), we obtain:

$$X = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{\frac{\sigma B\_0^2}{4\rho\_f} \xi\_f}$$

After using the value of *X*, we obtain:

$$k(\zeta) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{\frac{\sigma B\_0^2}{\alpha\_f}\zeta}$$

which implies that:

$$u\_1(y,t) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{\frac{\sigma B\_0^2}{4\rho\_f}(y-at)}.$$

This last expression shows the existence of an exponential profile along the travelling wave frame. This is not a trivial result for the nonlinear reaction under the Eyring–Powell fluid.

Note that the solution holds by the symmetry (*ζ* → −*ζ*) of travelling wave profiles. It suffices to admit *ζ* = *y* + *at*, so that:

$$k(\zeta) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{-\frac{\sigma B\_0^2}{4\rho\_f}\zeta}, \quad u\_1(y, t) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{-\frac{\sigma B\_0^2}{4\rho\_f}(y + at)}.\tag{39}$$

Now, it is the aim to show that the defined supporting manifold *N* preserves the exponential behaviour close to the critical points. For this purpose, the expression (36) is re-written as:

$$\frac{dY}{dX} = \frac{F}{Y} \left( X - \frac{A\_1 \rho\_f}{\sigma B\_0^2} \right). \tag{40}$$

After solving (40):

$$Y = F\left(X - \frac{A\_1 \rho\_f}{\sigma B\_0^2}\right). \tag{41}$$

From the expression (36), the equation (41) becomes:

$$X' = F\epsilon \left( X - \frac{A\_1 \rho\_f}{\sigma B\_0^2} \right). \tag{42}$$

After solving (42), we have:

$$X = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + \epsilon^{Fc\mathbb{Z}\_p}.\tag{43}$$

From (33), the expression (43) becomes:

$$k(\zeta) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{Fc\tilde{\zeta}}, \quad u\_1(y, t) = \frac{A\_1 \rho\_f}{\sigma B\_0^2} + e^{Fc(y - at)}.$$

This last expression permits showing the conservation of the exponential profile close to the critical points defined by the asymptotic manifolds *N*.

#### **6. Numerical Validation Assessments**

The aim in this section is to develop a numerical simulation to determine an appropriate travelling wave velocity (*a*) for which the approximated analytical solution (39) and the exact one, obtained numerically, in (34) behave similarly. This exercise can be seen as a validation process of the obtained analytical paths presented in the previous sections. This validation was explored for certain combinations of the fluid properties. Note that other combinations do not have an impact on the analytical ending in the exponential kind of solutions.

The numerical exploration was performed as per the following principles:


The results are compiled in Figures 1–3. The existence of a critical travelling wave speed *a*∗ = 1.212 for which the analytical solution in (39) is close to the numerically exact one of (34) with an accumulative error of <10−<sup>3</sup> was concluded. This numerical exploration permits accounting for the validation of the analytical exponential profile obtained.

**Figure 1.** *a* = 0.1 (**left**), *a* = 1 (**right**). The blue line is the exact numerical profile of the set of Equations (34). The red line is the analytical solution obtained in (39) up to *ζ* = 5 (beyond such values, it is required to change the scale). Solutions on the left are provided for *a* = 1 and solutions on the right for *a* = 1.5. For increasing values of the travelling speed, the solutions behave similarly in their exponential tail.

**Figure 2.** *a* = 1.212 (**left**), *a* = 1.5 (**right**). The blue line is the exact numerical profile of the set of Equations (34). The red line is the analytical solution obtained in (39). The approximated solution and the exact profile closely match an accumulative error (as the integration of the difference of both solutions) of < <sup>10</sup>−<sup>3</sup> for *<sup>a</sup>* <sup>=</sup> 1.212. Solutions on the right are given for *<sup>a</sup>* <sup>=</sup> 1.5. The approximated solution is above the numerical one.

**Figure 3.** *a* = 2 (**left**), *a* = 3 (**right**). The blue line is the exact numerical profile of the set of Equations (34). The red line is the analytical solution obtained in (39). Solutions on the left are provided for *a* = 2 and solutions on the right for *a* = 3. Note that for increasing values of the travelling wave speed, both profiles diverge.

#### **7. Conclusions**

The presented analysis in this article permitted accounting for the regularity, existence, and uniqueness of solutions to an Eyring–Powell fluid flow. Solutions were explored in the travelling wave domain, and asymptotic approaches were provided making use of the singular geometric perturbation theory. Afterwards, the obtained analytical solution was validated for a certain combination of fluid constants and making use of a numerical exercise. The existence of a travelling wave speed of *a* = 1.212 for which the analytical solution is close to the actual numerical solution with an accumulative error of <10−<sup>3</sup> was concluded. The existence of an exponential travelling wave tail together with a certain minimizing error critical speed constituted the main novelty reported by the present study.

**Author Contributions:** Conceptualisation, J.L.D.; methodology, J.L.D. and S.U.R.; validation, J.L.D.; formal analysis, J.L.D., S.U.R., J.C.S.R., M.A.S.R., G.F.C. and A.H.H.; investigation, J.L.D., S.U.R., J.C.S.R., M.A.S.R., G.F.C. and A.H.H.; resources, J.L.D.; data curation, J.L.D.; writing—original draft preparation, J.L.D., S.U.R., J.C.S.R., M.A.S.R., G.F.C. and A.H.H.; writing—review and editing, J.L.D., S.U.R., J.C.S.R., M.A.S.R., G.F.C. and A.H.H.; supervision, J.L.D.; project administration, J.L.D.; funding acquisition, J.L.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the University Francisco de Victoria School of Engineering.

**Data Availability Statement:** This research has no associated data.

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

#### **References**

