Mathematical analysis and numerical resolution of a heat transfer problem arising in water recirculation

TL;DR

This paper presents a mathematical model combining heat transfer and fluid dynamics for water recirculation, proves solution existence, and offers a numerical algorithm validated with realistic examples.

Contribution

It introduces a novel coupled heat transfer and Navier-Stokes model with turbulence, providing theoretical analysis and a computational algorithm.

Findings

Proved existence of solutions for the coupled model

Developed a numerical algorithm for simulation

Validated with realistic water recirculation examples

Abstract

This work is devoted to the analysis and resolution of a well-posed mathematical model for several processes involved in the artificial circulation of water in a large waterbody. This novel formulation couples the convective heat transfer equation with the modified Navier-Stokes system following a Smagorinsky turbulence model, completed with a suitable set of mixed, nonhomogeneous boundary conditions of diffusive, convective and radiative type. We prove several theoretical results related to existence of solution, and propose a full algorithm for its computation, illustrated with some realistic numerical examples.

Figures

Tables

Table 1. Table 1: Physical parameters for the numerical example.

Parameters	Values	Units
$\nu$	${\mathrm{1.3\; 10}}^{-3}$	$m^2{s}^{-1}$
${\nu }_{tur}$	${\mathrm{5.0\; 10}}^{-2}$	$m^2$
$K$	${\mathrm{1.4\; 10}}^{-5}$	$m^2{s}^{-1}$
$h^N$, $h^S$	${\mathrm{3.0\; 10}}^2$	$W{m}^{-2}{K}^{-1}$
$\rho$	${\mathrm{9.9\; 10}}^2$	$g{m}^{-3}$
$c_p$	$4.2$	$Ws{g}^{-1}{K}^{-1}$
${\theta }^0$	$283.0$	K
${\theta }_S$	$286.0$	K
${\theta }_N$	$283.0$	K
${\alpha }_0$	${\mathrm{8.7\; 10}}^{-7}$	$K^{-1}$

Equations

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\theta}{\partial t}+\mathbf{v}\cdot\nabla\theta-\nabla\cdot(K\nabla\theta)=0\quad\mbox{in}\;\Omega\times]0,T[,\\ \displaystyle\theta=\phi_{\theta}\quad\mbox{on}\;\Gamma_{T}\times]0,T[,\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=0\quad\mbox{on}\;\Gamma_{C}\times]0,T[,\vspace{0.1cm}\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=b_{1}^{N}(\theta_{N}-\theta)\quad\mbox{on}\;\Gamma_{N}\times]0,T[,\vspace{0.1cm}\\ \displaystyle K\frac{\partial\theta}{\partial\mathbf{n}}=b_{1}^{S}(\theta_{S}-\theta)+b_{2}^{S}(T_{r}^{4}-|\theta|^{3}\theta)\quad\mbox{on}\;\Gamma_{S}\times]0,T[,\\ \theta(0)=\theta^{0}\quad\mbox{in}\;\Omega,\end{array}\right.$$

$${\varphi }_{\theta }(\mathbf{x},t)=\sum _{k=1}^{N_{CT}}{\phi }^k(\mathbf{x}){\int }_{-T}^T{\rho }_{ϵ}(t-ϵ-s){\gamma }_{\theta }^k(s)ds$$

$$\gamma_{\theta}^{k}(s)=\left\{\begin{array}[]{lcl}\displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}\theta^{0}\,d\gamma&\mbox{ if }&s\leq 0,\\ \displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}\theta(s)\,d\gamma&\mbox{ if }&s>0,\displaystyle\end{array}\right.$$

$$\rho_{\epsilon}(t)=\left\{\begin{array}[]{lcl}\displaystyle\frac{c}{\epsilon}\,\exp\left(\frac{t^{2}}{t^{2}-\epsilon^{2}}\right)&\mbox{ if }&|t|<\epsilon,\\ \displaystyle 0&\mbox{ if }&|t|\geq\epsilon,\end{array}\right.$$

$${\int }_{\mathbb{R}}{\rho }_1(t)dt=1.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\mathbf{v}}{\partial t}+\nabla\mathbf{v}\,\mathbf{v}-\nabla\cdot\Xi(\mathbf{v})+\nabla p=\alpha_{0}(\theta-\theta^{0})\mathbf{a}_{g}\quad\mbox{in}\;\Omega\times]0,T[,\\ \displaystyle\nabla\cdot\mathbf{v}=0\quad\mbox{in}\;\Omega\times]0,T[,\\ \displaystyle\mathbf{v}=\boldsymbol{\phi}_{\mathbf{g}}\quad\mbox{on}\;\partial\Omega\times]0,T[,\\ \displaystyle\mathbf{v}(0)=\mathbf{v}^{0}\quad\mbox{in}\;\Omega,\end{array}\right.$$

$${\mathbf{\varphi }}_{\mathbf{g}}(t,\mathbf{x})=\sum _{k=1}^{N_{CT}}{g}^k(t)\left[\frac{{\phi }^k(\mathbf{x})}{\mu (T^k)}-\frac{{\tilde{\phi }}^k(\mathbf{x})}{\mu (C^k)}\right]\mathbf{n}.$$

$$\Xi (\mathbf{v})={\frac{\partial D(e)}{\partial e}|}_{e=e(\mathbf{v})},\text{mbox}withe(\mathbf{v})=\frac{1}{2}\left(\nabla \mathbf{v}+\nabla {\mathbf{v}}^t\right),$$

$$D(e)=\nu \left[e:e\right]+\frac{2}{3}{\nu }_{tur}{\left[e:e\right]}^{3/2},$$

$$\begin{array}[]{rcl}\displaystyle\Xi(\mathbf{v})&=&\displaystyle\left.\frac{\partial D(e)}{\partial e}\right|_{e=\epsilon(\mathbf{v})}=\displaystyle 2\nu\,\epsilon(\mathbf{v})+2\nu_{tur}\left[\epsilon(\mathbf{v}):\epsilon(\mathbf{v})\right]^{1/2}\epsilon(\mathbf{v})\\ &=&\displaystyle\left(2\nu+2\nu_{tur}\left[\epsilon(\mathbf{v}):\epsilon(\mathbf{v})\right]^{1/2}\right)\epsilon(\mathbf{v})=\displaystyle\beta(\epsilon(\mathbf{v}))\,\epsilon(\mathbf{v}),\end{array}$$

$$\begin{array}[]{rcl}X_{1}&=&\{\theta\in H^{1}(\Omega):\;\theta_{|_{\Gamma_{S}}}\in L^{5}(\Gamma_{S})\},\\ \widetilde{X}_{1}&=&\{\theta\in X_{1}:\;\theta_{|_{\Gamma_{T}}}=0\},\end{array}$$

$$\Vert \theta {\Vert }_{X_1}=\Vert \theta {\Vert }_{H^1(\Omega )}+\Vert \theta {\Vert }_{L^5({\Gamma }_S)}.$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{X}_{2}&=&\displaystyle\left\{\mathbf{v}\in[W^{1,3}(\Omega)]^{3}:\;\nabla\cdot\mathbf{v}=0,\;\mathbf{v}_{|_{\Gamma_{S}\cup\Gamma_{N}}}=\boldsymbol{0}\right\},\\ \displaystyle\widetilde{\mathbf{X}}_{2}&=&\displaystyle\left\{\mathbf{v}\in[W^{1,3}(\Omega)]^{3}:\;\nabla\cdot\mathbf{v}=0,\;\mathbf{v}_{|_{\partial\Omega}}=\boldsymbol{0}\right\}.\end{array}$$

$$W^{1,p,q}(0,T;V_1,V_2)=\left\{u\in L^p(0,T;V_1):\frac{du}{dt}\in L^q(0,T;V_2)\right\},$$

$$\begin{array}[]{rcl}\displaystyle W_{1}&=&\displaystyle\{\theta\in W^{1,2,5/4}(0,T;X_{1},X_{1}^{\prime}):\\ &&\displaystyle\theta_{|_{\Gamma_{S}}}\in L^{5}(0,T;L^{5}(\Gamma_{S}))\}\cap L^{\infty}(0,T;L^{2}(\Omega)),\\ \widetilde{W}_{1}&=&\displaystyle\{\theta\in W^{1,2,5/4}(0,T;\widetilde{X}_{1},\widetilde{X}_{1}^{\prime}):\\ &&\displaystyle\theta_{|_{\Gamma_{S}}}\in L^{5}(0,T;L^{5}(\Gamma_{S}))\}\cap L^{\infty}(0,T;L^{2}(\Omega)),\end{array}$$

$$\begin{array}[]{rcl}\displaystyle\mathbf{W}_{2}&=&\displaystyle W^{1,\infty,2}(0,T;\mathbf{X}_{2},[L^{2}(\Omega)]^{3})\cap\mathcal{C}([0,T];\mathbf{X}_{2}),\\ \displaystyle\widetilde{\mathbf{W}}_{2}&=&\displaystyle W^{1,\infty,2}(0,T;\widetilde{\mathbf{X}}_{2},[L^{2}(\Omega)]^{3})\cap\mathcal{C}([0,T];\widetilde{\mathbf{X}}_{2}).\end{array}$$

$$\begin{array}[]{rcl}R_{\mathbf{h}}:[L^{2}(0,T)]^{N_{CT}}&\rightarrow&W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))\\ \mathbf{h}&\rightarrow&R_{\mathbf{h}}(\mathbf{h})=\zeta_{\mathbf{h}},\end{array}$$

$${\zeta }_{\mathbf{h}}(\mathbf{x},t)=\sum _{k=1}^{N_{CT}}{\beta }_0({\phi }^k(\mathbf{x})){\int }_{-T}^T{\rho }_{ϵ}(t-ϵ-s){\gamma }_{\mathbf{h}}^k(s)ds,$$

$$\gamma^{k}_{\mathbf{h}}(s)=\left\{\begin{array}[]{lcl}\displaystyle\frac{1}{\mu(C^{k})}\int_{C^{k}}\theta^{0}\,d\gamma&\mbox{ if }&s\leq 0,\\ \displaystyle{h^{k}}(s)&\mbox{ if }&s>0,\end{array}\right.$$

$$\Vert {\zeta }_{\mathbf{h}}{\Vert }_{W^{1,2,2}(0,T;H^2(\Omega ),H^2(\Omega ))}\le C_1({\theta }^0)(1+\Vert \mathbf{h}{\Vert }_{[L^2(0,T){]}^{N_{CT}}}).$$

$$b_{\mathbf{h}}^k(t)={\int }_{-T}^T{\rho }_{ϵ}(t-ϵ-s){\gamma }_{\mathbf{h}}^k(s)ds,$$

$$\Vert {\beta }_0({\phi }^k){\Vert }_{H^2(\Omega )}\le C\Vert {\phi }^k{\Vert }_{H^{3/2}(\partial \Omega )},k=1,\dots ,N_{CT}.$$

$$\Vert b_{\mathbf{h}}^k{\Vert }_{H^1(0,T)}\le C\Vert {\gamma }_{\mathbf{h}}^k{\Vert }_{L^2(-T,T)}.$$

$$\begin{array}[]{rcl}\displaystyle\|\gamma_{\mathbf{h}}^{k}\|_{L^{2}(-T,T)}^{2}&=&\displaystyle\int_{-T}^{T}\left(\gamma_{\mathbf{h}}^{k}(s)\right)^{2}\,ds\\ &=&\displaystyle\frac{1}{\mu(C^{k})^{2}}\int_{-T}^{0}\left(\int_{C^{k}}\theta^{0}\,d\gamma\right)^{2}+\int_{0}^{T}\left(h^{k}(s)\right)^{2}\,ds\\ &\leq&\displaystyle\frac{T}{\mu(C^{k})}\|\theta^{0}\|_{L^{2}(\partial\Omega)}^{2}+\|h^{k}\|^{2}_{L^{2}(0,T)}.\end{array}$$

$$\Vert {\zeta }_{\mathbf{h}}{\Vert }_{W^{1,2,2}(0,T;H^2(\Omega ),H^2(\Omega ))}$$

$$W^{1,2,2}(0,T;H^2(\Omega ),H^2(\Omega ))\subset \subset W_1.$$

$$h^k(s)=\frac{1}{\mu (C^k)}{\int }_{C^k}\theta (s)d\gamma ,k=1,2,\dots ,N_{CT}.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\xi}{\partial t}+\mathbf{v}\cdot\nabla\xi-\nabla\cdot(K\nabla\xi)\\ \displaystyle\quad=-\frac{\partial\zeta_{\mathbf{h}}}{\partial t}-\mathbf{v}\cdot\nabla\zeta_{\mathbf{h}}+\nabla\cdot(K\nabla\zeta_{\mathbf{h}})\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle\xi=0\quad\mbox{on}\;T^{k}\times(0,T),\quad\mbox{for}\;k=1,\ldots,N_{CT},\\ \displaystyle K\frac{\partial\xi}{\partial\mathbf{n}}=-K\frac{\partial\zeta_{\mathbf{h}}}{\partial\mathbf{n}}\quad\mbox{on}\;C^{k}\times(0,T),\quad\mbox{for}\;k=1,\ldots,N_{CT},\\ \displaystyle K\frac{\partial\xi}{\partial\mathbf{n}}=b_{1}^{N}\big{(}\theta_{N}-\zeta_{\mathbf{h}}-\frac{K}{b_{1}^{N}}\frac{\partial\zeta_{\mathbf{h}}}{\partial\mathbf{n}}-\xi\big{)}\quad\mbox{on}\;\Gamma_{N}\times(0,T),\\ \displaystyle K\frac{\partial\xi}{\partial\mathbf{n}}=b_{1}^{S}\big{(}\theta_{S}-\zeta_{\mathbf{h}}-\frac{K}{b_{1}^{S}}\frac{\partial\zeta_{\mathbf{h}}}{\partial\mathbf{n}}-\xi\big{)}\\ \displaystyle\quad+b_{2}^{S}\big{(}T_{r}^{4}-|\xi+\zeta_{\mathbf{h}}|^{3}(\xi+\zeta_{\mathbf{h}})\big{)}\quad\mbox{on}\;\Gamma_{S}\times(0,T),\\ \displaystyle\xi(0)=\theta^{0}-\zeta_{\mathbf{h}}(0)\quad\mbox{in}\;\Omega.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\mathbf{z}}{\partial t}+\nabla(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})\mathbf{z}+\nabla\mathbf{z}\boldsymbol{\zeta}_{\mathbf{g}}\\ \displaystyle\quad-\text{div}\left(2\nu\epsilon(\mathbf{z})+2\nu_{tur}\int_{\Omega}\left[\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}):\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})\right]^{1/2}\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})\right)\\ \displaystyle\quad+\nabla p=\alpha_{0}(\theta-\theta^{0})\,\mathbf{a}_{g}-\frac{\partial\boldsymbol{\zeta}_{\mathbf{g}}}{\partial t}-\nabla\boldsymbol{\zeta}_{\mathbf{g}}\boldsymbol{\zeta}_{\mathbf{g}}\\ \quad\displaystyle+2\nu\nabla\cdot\epsilon(\boldsymbol{\zeta}_{\mathbf{g}})\quad\mbox{in}\;\Omega\times(0,T),\\ \displaystyle\mathbf{z}=\mathbf{0}\quad\mbox{on}\;\partial\Omega\times(0,T),\\ \displaystyle\mathbf{z}(0)=\mathbf{v}^{0}\quad\mbox{in}\;\Omega.\end{array}\right.$$

$$h^k(s)=\frac{1}{\mu (C^k)}{\int }_{C^k}\theta (s)d\gamma ,k=1,2,\dots ,N_{CT}.$$

Full text

,

,

1 Introduction

Artificial circulation in large waterbodies is a management technique aimed to disrupt stratification of temperature and, consequently, to minimize the development of stagnant zones that may be subject to water quality problems (for instance, low levels of dissolved oxygen or high concentrations of phytoplankton). For its operation, a set of flow pumps take water from the upper layers by means of collectors and inject it into the bottom layers, setting up a recirculation pattern that prevents stratification by means of a forced mixing of water. One of the main problems of the temperature stratification is related to algal blooms produced in the upper layers due to high temperature and solar radiation. However, if we circulate water from the bottom layers (where the temperature is lower) to the upper layers, we can mitigate this negative effect. Further details and remarks on several issues related to the optimal design and control of water artificial circulation techniques have been analyzed by the authors in their recent work [19].

Convective heat transfer has been the subject of an intensive mathematical research in last five decades (ranging, for instance, from the pioneering works on the Boussinesq system of Joseph [15] in the 1960s to the present). Among the recent contributions we must mention, for instance, some papers devoted to study related problems in the steady case [4, 16], the analysis a time-dependant case, but not including convective phenomena, [20], and some numerical approaches [3, 18]. Nevertheless, after an exhaustive search we have not been able to find in the mathematical literature the analysis of the particular problem arising in the setting of our water recirculation model: a coupled problem linking a heat equation with mixed nonlinear boundary conditions to a modified Navier-Stokes equation following the Smagorinsky model of turbulence. Thus, the present work deals with the mathematical analysis and the numerical resolution of this heat transfer problem with specific boundary conditions related to water artificial circulation in a body of water (for instance, a lake or a reservoir). The main difficulties in the study of this problem lie in the nonlinear boundary condition related to the solar irradiation on the surface, the relations between the water temperature in the collectors and the injectors, and the coupling between water temperature and water velocity due to convective effects. We use the Smagorinsky model of turbulence instead of other approaches, like the celebrated $k-\epsilon$ system, due to the fact that the modified Navier-Stokes equations following the Smagorinsky model of turbulence present very interesting properties from a mathematical viewpoint, in particular the uniqueness of solution and its additional regularity.

The organization of this paper is as follows: First we introduce a well-posed formulation of the physical problem and present a rigorous definition of a solution for the problem. In the central part of the paper we prove the existence of this solution, and in the final part we propose a numerical algorithm for its resolution, showing several computational tests for a realistic example. At the end of the paper we include an appendix with several results for a general heat equation with an advective term and mixed boundary conditions of diffusive, convective and radiative type. These results, as far as we know, are new since we use techniques that will allow us to treat the low regularity of the time derivative of the solution. This lack of regularity represents, together with the nonlinearity of the boundary conditions, the main difficulty.

2 Mathematical formulation of the problem

In this section we present in detail the three-dimensional mathematical model under study. So, we consider a convex domain $\Omega\subset\mathbb{R}^{3}$ (representing the waterbody) whose boundary surface $\partial\Omega$ can be split into four smooth enough, disjoint sections: $\Gamma_{S}$, $\Gamma_{C}$, $\Gamma_{T}$ and $\Gamma_{N}$, in such a way that $\partial\Omega=\Gamma_{S}\cup\Gamma_{C}\cup\Gamma_{T}\cup\Gamma_{N}$. Subset $\Gamma_{S}$ represents the part of the boundary in contact with air, $\Gamma_{C}$ is the part of the boundary where the collectors are located, $\Gamma_{T}$ is the part of the boundary where the injectors are located, and $\Gamma_{N}$ stands for the rest of the boundary. We will suppose that each collector is linked to an injector by means of a pumped pipeline, and we also assume that there exist $N_{CT}$ collector/injector pairs $\{(C^{k},T^{k})\}_{k=1}^{N_{CT}}$. Therefore, $\Gamma_{C}=\cup_{k=1}^{N_{CT}}C^{k}$, and $\Gamma_{T}=\cup_{k=1}^{N_{CT}}T^{k}$. In Fig. 1 we can see a schematic geometrical configuration of a rectangular domain $\Omega$ for a particular case of $N_{CT}=4$ collector/injector pairs.

As above commented, we suppose the boundary $\partial\Omega$ regular enough to assure the existence of elements $\varphi^{k},\,\widetilde{\varphi}^{k}\in H^{3/2}(\partial\Omega)$, $k=1,\ldots,N_{CT}$, satisfying the following assumptions (corresponding to suitable regularizations of the indicator functions of $T^{k}$ and $C^{k}$, respectively):

• •

$\varphi^{k}(\mathbf{x}),\,\widetilde{\varphi}^{k}(\mathbf{x})\geq 0$, a.e. $\mathbf{x}\in\partial\Omega$,

• •

$\varphi^{k}(\mathbf{x})=0$, a.e. $\mathbf{x}\in\partial\Omega\setminus T^{k}$, and $\displaystyle\int_{T^{k}}\varphi^{k}(\mathbf{x})\,d\gamma=\mu(T^{k})$,

• •

$\widetilde{\varphi}^{k}(\mathbf{x})=0$, a.e. $\mathbf{x}\in\partial\Omega\setminus C^{k}$, and $\displaystyle\int_{C^{k}}\widetilde{\varphi}^{k}(\mathbf{x})\,d\gamma=\mu(C^{k})$,

where $\mu(S)$ represents the area measure of any set $S\subset\partial\Omega$.

We denote by $\theta(\mathbf{x},t)$ (measured in K) the solution of the following convection-diffusion partial differential equation with nonhomogeneous, nonlinear, mixed boundary conditions:

[TABLE]

where Dirichlet boundary condition $\phi_{\theta}$ is given by expression:

[TABLE]

with, for each $k=1,\ldots,N_{CT}$,

[TABLE]

representing the mean temperature of water in the collector $C_{k}$, and with the weight function $\rho_{\epsilon}$ defined by:

[TABLE]

for $c\in\mathbb{R}$ the positive constant satisfying the unitary condition:

[TABLE]

In other words, we are assuming that the mean temperature of water at each injector $T_{k}$ is a weighted average in time of the mean temperatures of water at its corresponding collector $C_{k}$. In order to obtain the mean temperature at each injector, we convolute the mean temperature at the collector with a smooth function with support in $(t-2\epsilon,t)$. In this way, we have that the temperature in the injector only depends on the mean temperature in the collector in the time interval $(t-2\epsilon,t)$. Parameter $0<\epsilon<T$ represents, in a certain sense, the technical characteristics of the pipeline that define the stay time of water in the pipe. We also suppose that there is not heat transfer thought the walls of the pipelines (that is, they are isolated).

Moreover,

• •

$T>0\ ({\rm s})$ is the length of the time interval.

• •

$\mathbf{n}$ is the unit outward normal vector to the boundary $\partial\Omega$.

• •

$K>0\ ({\rm m}^{2}\,{\rm s}^{-1})$ is the thermal diffusivity of the fluid: $K=\frac{\alpha}{\rho\,c_{p}}$, where $\alpha\ ({\rm W}\,{\rm m}^{-1}\,{\rm K}^{-1})$ is the thermal conductivity, $\rho\ ({\rm g}\,{\rm m}^{-3})$ is the density, and $c_{p}\ ({\rm W}\,{\rm s}\,{\rm g}^{-1}\,{\rm K}^{-1})$ is the specific heat capacity of water.

• •

$b_{1}^{K}\geq 0\ ({\rm m}\,{\rm s}^{-1})$, for $K=N,S$, are the coefficients related to convective heat transfer through the boundaries $\Gamma_{N}$ and $\Gamma_{S}$, obtained from the relation $\rho\,c_{p}\,b_{1}^{K}=h^{K}$, where $h^{K}\geq 0\ ({\rm W}\,{\rm m}^{-2}\,{\rm K}^{-1})$ are the convective heat transfer coefficients on each surface. These coefficients are relevant in the convective heat transfer flux through the frontiers $\Gamma_{S}$ and $\Gamma_{N}$, $b_{1}^{K}(\theta_{K}-\theta)$, $K=N,S$.

• •

$b_{2}^{S}>0\ ({\rm m}\,{\rm s}\,{\rm K}^{-3})$ is the coefficient related to radiative heat transfer through the boundary $\Gamma_{S}$, given by $b_{2}^{S}=\frac{\sigma_{B}\,\varepsilon}{\rho\,c_{p}}$, where $\sigma_{B}\ ({\rm W}\,{\rm m}^{-2}\,{\rm K}^{-4})$ is the Stefan-Boltzmann constant and $\varepsilon$ is the emissivity. This coefficient is fundamental in the radiative flux through the frontier $\Gamma_{S}$ (see, for instance, the classical reference [5] for a complete description of this type of boundary conditions).

• •

$\theta^{0}\geq 0\ ({\rm K})$ is the initial temperature.

• •

$\theta_{S},\,\theta_{N}\geq 0\ ({\rm K})$ are the temperatures related to convection heat transfer on the surfaces $\Gamma_{S}$ and $\Gamma_{N}$.

• •

$T_{r}\geq 0\ ({\rm K})$ is the radiation temperature on the surface $\Gamma_{S}$, derived from expression $\sigma_{B}\,\varepsilon\,T_{r}^{4}=(1-a)R_{sw,net}+R_{lw,down}$, where $a$ is the albedo, $R_{sw,net}\ ({\rm W}\,{\rm m}^{-2})$ denotes the net incident shortwave radiation on the surface $\Gamma_{R}$, and $R_{lw,down}\ ({\rm W}\,{\rm m}^{-2})$ denotes the downwelling longwave radiation.

Finally, $\mathbf{v}(\mathbf{x},t)\ ({\rm m}\,{\rm s}^{-1})$ is the water velocity, solution of a modified Navier-Stokes equations following a Smagorinsky model of turbulence:

[TABLE]

where $\mathbf{a}_{g}=g\,\mathbf{e}_{3}$ $({\rm m}\,{\rm s}^{-1})$ is the gravity acceleration, $\alpha_{0}=-\frac{1}{\rho}\frac{\partial\rho}{\partial\theta}$ $({\rm K}^{-1})$ is the thermic expansion coefficient. We must remark here that we are assuming the thermodynamic process to be close to an initial equilibrium state that we denote with the zero subscript, so $\rho=\rho_{0}-\alpha_{0}\rho_{0}(\theta-\theta_{0})+k_{0}\rho_{0}(p-p_{0})+o(\theta-\theta_{0})+o(p-p_{0})$, where $k_{0}=\frac{1}{\rho}\frac{\partial\rho}{\partial p}$ is the isothermal compressibility coefficient. The details of this approach (known as the Boussinesq model for natural convection) can be consulted, for instance, in Section 10.7 of [6]. Finally, $\mathbf{v}^{0}$ is the initial velocity, and boundary field $\boldsymbol{\phi}_{\mathbf{g}}$ is the element given by:

[TABLE]

with, for each $k=1,\ldots,N_{CT}$, $g^{k}(t)\in H^{1}(0,T)$, representing the volumetric flow rate by pump $k$ at each time $t$ ($g^{k}(t)>0$, $\forall t\in]0,T[$, and $g^{k}(0)=0$). The turbulence term $\Xi(\mathbf{v})$ is given by:

[TABLE]

where $D$ is a potential function (for instance, in the particular case of the classical Navier-Stokes equations, $D(e)=\nu\left[e:e\right]$, with $\nu$ $({\rm m}^{2}\,{\rm s}^{-1})$ the kinematic viscosity of the water, and, consequently, $\Xi(\mathbf{v})=2\nu\,e(\mathbf{v})$). However, in our case, the Smagorinsky model, the potential function is defined as [17]:

[TABLE]

where $\nu_{tur}$ $({\rm m}^{2})$ is the turbulent viscosity. So, for the Smagorinsky case,

[TABLE]

with $\beta(\epsilon(\mathbf{v}))=2\nu+2\nu_{tur}\left[\epsilon(\mathbf{v}):\epsilon(\mathbf{v})\right]^{1/2}$.

System (6) has been recently studied by the authors in [9], where it is demonstrated the existence and the uniqueness of solution for a recirculation model based in the modified Navier-Stokes equations. In the present work we will use some of the results shown in [9] in order to prove the existence of solution for the coupled problem (1) and (6). The main difficulties of this work relies in the coupling of a heat equation with nonlinear boundary value terms and the modified Navier-Stokes system. The nonlinear terms do not allow us to work with regular solutions for the heat equation, which forces us to use more sophisticated techniques in order to demonstrate existence and uniqueness of solution.

3 The concept of solution

We start this section defining the functional spaces used in the definition of solution for the system (1) and (6). So, for the water temperature we consider:

[TABLE]

and we define the following norm associated to above space $X_{1}$:

[TABLE]

We have that $X_{1}$ is a reflexive separable Banach space (cf. Lemma 3.1 of [8]) and $\widetilde{X}_{1}\subset L^{2}(\Omega)\subset\widetilde{X}_{1}^{\prime}$ is an evolution triple. For the water velocity we consider:

[TABLE]

In order to define an appropriate space for the solution of problems (1) and (6), we consider, for a Banach space $V_{1}$ and a locally convex space $V_{2}$ such that $V_{1}\subset V_{2}$, the following Sobolev-Bochner space (cf. Chapter 7 of [21]), for $1\leq p,q\leq\infty$:

[TABLE]

where $\frac{du}{dt}$ denotes the derivative of $u$ in the sense of distributions. It is well known that, if both $V_{1}$ and $V_{2}$ are Banach spaces, then $W^{1,p,q}(0,T;V_{1},V_{2})$ is also a Banach space endowed with the norm $\|u\|_{W^{1,p,q}(0,T;V_{1},V_{2})}=\|u\|_{L^{p}(0,T;V_{1})}+\left\|\frac{du}{dt}\right\|_{L^{q}(0,T;V_{2})}$.

Then, we define the following spaces that will be used in the mathematical analysis of system (1):

[TABLE]

and, for the system (6), we define

[TABLE]

Hypothesis 1

We will assume the following hypotheses for coefficients and data of the problem:

1. (a)

$\theta^{0}\in X_{2}$** 2. (b)

$\theta_{S}\in L^{2}(0,T;L^{2}(\Gamma_{S}))$** 3. (c)

$\theta_{N}\in L^{2}(0,T;L^{2}(\Gamma_{N}))$** 4. (d)

$T_{r}\in L^{5}(0,T;L^{5}(\Gamma_{S}))$** 5. (e)

$\mathbf{v}^{0}\in\left[H_{\sigma}^{2}(\Omega)\right]^{3}=\{\mathbf{v}\in[H^{2}(\Omega)]^{3}:\;\nabla\cdot\mathbf{v}=0,\ \mathbf{v}_{|_{\partial\Omega}}=\boldsymbol{0}\}\subset\widetilde{\mathbf{X}}_{2}$** 6. (f)

$g^{k}\in H^{1}(0,T)$* with $g^{k}(0)=0$, $\forall k=1,\ldots,N_{CT}$*

Remark 2

In order to define in a rigorous way the concept of solution, we will need to extend Dirichlet conditions of $\theta$ and $\mathbf{v}$ to the whole domain $\Omega$.

So, for water velocity $\mathbf{v}$, thanks to Lemma 2 of [9], for each $\mathbf{g}\in[H^{1}(0,T)]^{N_{CT}}$, there exists an element $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;\left[H_{\sigma}^{2}(\Omega)\right]^{3},\left[H_{\sigma}^{2}(\Omega)\right]^{3})$ such that ${\boldsymbol{\zeta}_{\mathbf{g}}}_{|_{\partial\Omega}}=\boldsymbol{\phi}_{\mathbf{g}}$, with $\boldsymbol{\phi}_{\mathbf{g}}$ defined by (7). Besides, by Lemma 3 of [9], $W^{1,2,2}(0,T;\left[H_{\sigma}^{2}(\Omega)\right]^{3},$ $\left[H_{\sigma}^{2}(\Omega)\right]^{3})\subset W^{1,\infty,2}(0,T;[H^{2}(\Omega)]^{3},[H^{2}(\Omega)]^{3})\cap\mathcal{C}([0,T];[H^{2}(\Omega)]^{3})$ and, then, we can use this element to reformulate the original problem for $\mathbf{v}$ as an homogeneous Dirichlet boundary condition one.

For water temperature $\theta$ we can proceed in an analogous way and prove that there exists an extension that allows us to reformulate the problem for $\theta$ as one with homogeneous boundary conditions.

Lemma 3

We have that the following operator is compact

[TABLE]

where:

[TABLE]

with $\gamma^{k}_{\mathbf{h}}(s)\in L^{2}(-T,T)$, for $k=1,2,\ldots,N_{CT}$, defined by:

[TABLE]

and $\beta_{0}:u\in H^{3/2}(\partial\Omega)\rightarrow\beta_{0}(u)\in H^{2}(\Omega)$ the right inverse of the classical trace operator $\gamma_{0}$ (that is, $(\gamma_{0}\circ\beta_{0})(u)=u$.)

We also have that there exists a constant $C_{1}$, that depends continuously on the space-time configuration of our computational domain and $\theta^{0}$, such that:

[TABLE]

{@proof}

[Proof.] Let $\{\mathbf{h}_{n}\}_{n\in\mathbb{N}}$ be a bounded sequence in $[L^{2}(0,T)]^{N_{CT}}$. Then, taking subsequences if necessary, we have that $\mathbf{h}_{n}\rightharpoonup\mathbf{h}$ weakly in $[L^{2}(0,T)]^{N_{CT}}$. We also have that $\gamma_{\mathbf{h}_{n}}\rightharpoonup\gamma_{\mathbf{h}}$ weakly in $[L^{2}(-T,T)]^{N_{CT}}$ and, if we denote by

[TABLE]

we obtain that $b_{\mathbf{h}_{n}}^{k}(t)\rightarrow b_{\mathbf{h}}^{k}(t)$ pointwise a.e. $t\in[0,T]$, $\forall k=1,\ldots,N_{CT}$. Thus, the sequence $\{\mathbf{b}_{\mathbf{h}_{n}}\}_{n\in\mathbb{N}}$ is bounded by a function in $[L^{2}(0,T)]^{N_{CT}}$, so we have the strong convergence in $[L^{2}(0,T)]^{N_{CT}}$. We can repeat the same argument with the time derivative of $\{\mathbf{b}_{\mathbf{h}_{n}}\}_{n\in\mathbb{N}}$, obtaining that $\mathbf{b}_{\mathbf{h}_{n}}\rightarrow\mathbf{b}_{\mathbf{h}}$ strongly in $H^{1}(0,T)$, thus $\zeta_{\mathbf{h}_{n}}\rightarrow\zeta_{\mathbf{h}}$ strongly in $W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$. Finally, by the properties of the operator $\beta_{0}$ we have that

[TABLE]

So, thanks to the regularity of function $\rho_{\epsilon}$, it is clear that $b_{\mathbf{h}}^{k}\in H^{1}(0,T)$, and

[TABLE]

In the other hand,

[TABLE]

Thereby, we have the following inequality:

[TABLE]

where $C(\theta^{0})$ is a positive constant than depends continuously on the spatial-time configuration of our computational domain and on the initial temperature. Moreover, it is worthwhile remarking here that we can make this constant as small as we want by considering data appropriately. $\blacksquare$

Finally, the following technical lemma will be necessary in order to guaranty that the sum of an element of $W_{1}$ plus an element of $W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ makes sense.

Lemma 4

We have that the following inclusion is compact:

[TABLE]

Now, we define the concept of solution for coupled system (1) and (6) in terms of homogeneous Dirichlet systems, so we establish the following notations:

• •

$\xi=\theta-\zeta_{\mathbf{h}}\in\widetilde{W}_{1}$, with $\displaystyle\zeta_{\mathbf{h}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ the extension obtained from Lemma 3, where:

[TABLE]

• •

$\mathbf{z}=\mathbf{v}-\boldsymbol{\zeta}_{\mathbf{g}}\in\widetilde{\mathbf{W}}_{2}$, with $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$ the extension of the trace given in Lemma 2 of [9].

Thus, using above notations, we can reformulate the state system (1) and (6) in the following way:

[TABLE]

[TABLE]

It is worthwhile remarking here that above system shows homogeneous Dirichlet boundary conditions and, consequently, we will be able to define the concept of solution of the original state systems (1) and (6) in terms of the modified state systems (27) and (28). It should be also noted that, in the case of equation (1), the coupling terms in the Dirichlet boundary conditions are now transferred to the partial differential equation in system (27).

Definition 5

A pair $(\theta,\mathbf{v})\in W_{1}\times\mathbf{W}_{2}$ is said to be a solution of problem (1) and (6) if there exist elements $(\xi,\mathbf{z})\in\widetilde{W}_{1}\times\widetilde{\mathbf{W}}_{2}$ such that:

$\mathbf{v}=\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}$, with $\displaystyle\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H^{2}(\Omega)]^{3},[H^{2}(\Omega)]^{3})$ the reconstruction of the trace given in Lemma 2 of **[9]**, and $\theta=\zeta_{\mathbf{h}}+\xi$, with $\displaystyle\zeta_{\mathbf{h}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ the extension obtained in Lemma 3, where:

[TABLE] 2. 2.

$\mathbf{z}(0)=\mathbf{v}^{0}$, and $\xi(0)=\theta^{0}-\zeta_{\mathbf{h}}(0)$, a.e. $\mathbf{x}\in\Omega$. 3. 3.

$(\xi,\mathbf{z})$* verifies the following variational formulation:*

[TABLE]

[TABLE]

where

[TABLE]

4 Existence of solution

We will prove now that, under certain hypotheses over coefficients and data, there exists a unique solution for the system (1) and (6) in the sense of Definition 5. The procedure used here for demonstrating the existence of solution is based in the Schauder fixed point Theorem (cf. section 9.5 of [7]) and is similar to one employed by the authors, for instance, in [10]. The main difficulties in the present case lie in the coupling of the Dirichlet conditions for the water temperature, and in the nonlinear radiation terms. To overcome these difficulties we will need to define a compact extension for the nonhomogeneous Dirichlet conditions and to prove novel results for the heat equation with radiation boundary conditions (see Appendix A).

So, we consider the following operator:

[TABLE]

where:

• •

$\mathbf{v}\in\mathbf{W}_{2}$ is such that $\mathbf{z}=\mathbf{v}-\boldsymbol{\zeta}_{\mathbf{g}}\in\widetilde{\mathbf{W}}_{2}$ is the solution of the following problem:

[TABLE]

with $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$ the extension obtained from the Lemma 2 of [9].

• •

$\theta\in W_{1}$ is such that $\xi=\theta-\zeta_{\mathbf{h}^{*}}\in\widetilde{W}_{1}$ is the solution of:

[TABLE]

with $\zeta_{\mathbf{h}^{*}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ defined as in Lemma 3.

• •

$\mathbf{h}\in[L^{2}(0,T)]^{N_{CT}}$ is such that:

[TABLE]

The following technical results are necessary to prove that the operator $M$ defined in (33) is well defined. The first one corresponds to the existence of solution for problem (34), and the second one is related to the existence of solution for problem (35).

Theorem 6

Within the framework stablished in Hypothesis 1, given elements $\boldsymbol{\zeta}_{\mathbf{g}}\in W^{1,2,2}(0,T;[H_{\sigma}^{2}(\Omega)]^{3},[H_{\sigma}^{2}(\Omega)]^{3})$ and $\theta^{*}\in L^{2}(0,T;L^{2}(\Omega))$, there exists an element $\mathbf{v}\in\mathbf{W}_{2}$ such that $\mathbf{z}=\mathbf{v}-\boldsymbol{\zeta}_{\mathbf{g}}\in\widetilde{\mathbf{W}}_{2}$ is the unique solution of problem (34) in the following sense:

[TABLE]

with $\mathbf{z}(0)=\mathbf{v}^{0}$, a.e. $\mathbf{x}\in\Omega$, where:

[TABLE]

Besides, we have the following estimates:

[TABLE]

[TABLE]

where $C_{2}$ and $C_{3}$ are positive constants that depend continuously with respect to the space-time configuration of our computational domain, $\mathbf{v}^{0}$, $\theta^{0}$ and $\mathbf{g}$.

{@proof}

[Proof.] It is a direct consequence of Theorem 8 in [9]. In this theorem the authors prove the wellposedness (existence, uniqueness and regularity of solution) for a modified Navier-Stokes system with non-homogeneous Dirichlet boundary conditions, by building a continuous extension of the Dirichlet condition and a Galerking approximation of the corresponding problem with homogeneous boundary conditions. $\blacksquare$

Theorem 7

Within the framework stablished in Hypothesis 1, given elements $\zeta_{\mathbf{h}^{*}}\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$ and $\mathbf{v}\in L^{10/3}(0,T;[L^{3}_{\sigma}(\Omega)]^{3})$, there exits an element $\theta\in W_{1}$ such that $\xi=\theta-\zeta_{\mathbf{h}^{*}}\in\widetilde{W}_{1}$ is the unique solution of problem (35) in the following sense:

[TABLE]

with $\xi(0)=\theta^{0}-\zeta_{\mathbf{h}^{*}}(0)$, a.e. $\mathbf{x}\in\Omega$.

Besides, we have the following estimates:

[TABLE]

[TABLE]

where $C_{4}$ and $C_{5}$ are positive constants that depend continuously with respect to the space-time configuration of our computational domain, $\theta^{0}$, $\theta_{N}$, $\theta_{S}$ and $T_{r}$.

{@proof}

[Proof.] It is a straightforward consequence of Theorem 14 that we will prove in Appendix A. In our case, we have to choose there $\Gamma_{R}=\Gamma_{S}$, $\Gamma_{L}=\Gamma_{T}$ and $\Gamma_{A}=\Gamma_{C}\cup\Gamma_{N}$. $\blacksquare$

Lemma 8

The operator $M$ defined in (33) is well defined and compact.

{@proof}

[Proof.] Thanks to above Theorems 6 and 7, it is straightforward that the operator $M$ is well defined. Let us check now its compactness.

So, given a bounded sequence $\{(\theta_{n}^{*},\mathbf{h}_{n}^{*})\}_{n\in\mathbb{N}}\subset L^{2}(0,T;L^{2}(\Omega))\times[L^{2}(0,T)]^{N_{CT}}$, we have, using estimates (39) and (40), that the corresponding sequence $\{\mathbf{z}_{n}\}_{n\in\mathbb{N}}\subset\widetilde{\mathbf{W}}_{2}$ of solutions for the problem (37) is bounded in $\widetilde{\mathbf{W}}_{2}$. Then, we have, taking subsequences if necessary, that:

• •

$\theta_{n}^{*}\rightharpoonup\theta$ in $L^{2}(0,T;L^{2}(\Omega))$,

• •

$\mathbf{h}_{n}^{*}\rightharpoonup\mathbf{h}$ in $[L^{2}(0,T)]^{N_{CT}}$,

• •

$\mathbf{z}_{n}\rightarrow\mathbf{z}$ strongly in $L^{p}(0,T;[L^{q}(\Omega)]^{3})$, for all $1<p<\infty$, and $2\leq q<\infty$,

• •

$\displaystyle\mathbf{z}_{n}\rightharpoonup\mathbf{z}$ weakly in $L^{3}(0,T;\widetilde{\mathbf{X}})$,

• •

$\displaystyle\frac{d\mathbf{z}_{n}}{dt}\rightharpoonup\frac{d\mathbf{z}}{dt}$ weakly in $L^{2}(0,T;[L^{2}(\Omega)]^{3})$,

• •

$\displaystyle\nabla\mathbf{z}_{n}\rightharpoonup^{*}\nabla\mathbf{z}$ weakly-$*$ in $L^{\infty}(0,T;[L^{3}(\Omega)]^{3})$,

• •

$\displaystyle\beta(\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}_{n}))\,\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}_{n})\rightharpoonup\beta(\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z}))\,\epsilon(\boldsymbol{\zeta}_{\mathbf{g}}+\mathbf{z})$ weakly in $L^{3/2}(0,T;\widetilde{\mathbf{X}}^{\prime})$,

where $\mathbf{z}\in\widetilde{\mathbf{W}}_{2}$ is the solution of (37) associated to $\theta$. The last convergence is a consequence of the monotony of operator (see [9] for more details)

[TABLE]

where, for any $\boldsymbol{\xi}\in\widetilde{\mathbf{X}}_{2}$,

[TABLE]

Now, using the results proved in Lemma 3, we have that the corresponding sequence $\{\zeta_{\mathbf{h}_{n}^{*}}\}$ converges to $\boldsymbol{\zeta}_{\mathbf{h}}$ strongly in $W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$, and, consequently,

• •

$H_{\mathbf{h}_{n}^{*}}\rightarrow H_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\Omega))$,

• •

$g^{C}_{\mathbf{h}_{n}^{*}}\rightarrow g^{C}_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\Gamma_{C}))$,

• •

$g^{N}_{\mathbf{h}_{n}^{*}}\rightarrow g^{N}_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\Gamma_{N}))$,

• •

$g^{S}_{\mathbf{h}_{n}^{*}}\rightarrow g^{S}_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\Gamma_{S}))$.

Finally, thanks to estimates (42) and (43), the corresponding sequence $\{\xi_{n}\}_{n\in\mathbb{N}}$ $\subset\widetilde{W}_{1}$ is bounded. Thus, taking subsequences if necessary, we have thanks to a straightforward adaptation of Lemma 19 in Appendix A, that

• •

$\displaystyle\xi_{n}\rightharpoonup\xi$ in $L^{2}(0,T;\widetilde{X})$,

• •

$\displaystyle\xi_{n}\rightharpoonup^{*}\xi$ in $L^{\infty}(0,T;L^{2}(\Omega))$,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{10/3-\epsilon}(0,T;L^{10/3-\epsilon}(\Omega))$,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{2}(0,T;L^{2}(\Gamma_{C}))$,

• •

$\displaystyle\xi_{n}\rightarrow\xi$ in $L^{4}(0,T;L^{4}(\Gamma_{R}))$.

Using the same techniques that we present in Appendix A for the demonstration of Theorem 14, we can pass to the limit in the variational formulation of $\xi_{n}$ and prove that $\xi\in\widetilde{W}_{2}$ is the solution of (7) associated to $\theta$ and $\mathbf{h}$. Thus, we have that

• •

$\displaystyle\theta_{n}=\xi_{n}+\zeta_{\mathbf{h}_{n}^{*}}\rightarrow\theta=\xi+\zeta_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\Omega))$,

• •

$\displaystyle\theta_{n}=\xi_{n}+\zeta_{\mathbf{h}_{n}^{*}}\rightarrow\theta=\xi+\zeta_{\mathbf{h}}$ in $L^{2}(0,T;L^{2}(\partial\Omega))$,

and, consequently, $M(\theta_{n}^{*},\mathbf{h}_{n}^{*})\rightarrow M(\theta,\mathbf{h})$ in $L^{2}(0,T;L^{2}(\Omega))\times[L^{2}(0,T)]^{N_{CT}}$, which concludes the proof. $\blacksquare$

Theorem 9

Given positive constants $\widehat{C}^{1}$ and $\widehat{C}^{2}$, there exist coefficients and data small enough such that the operator $M$ defined in (33) has a fixed point in the space $\{(\theta,\mathbf{h})\in L^{2}(0,T;L^{2}(\Omega))\times[L^{2}(0,T)]^{N_{CT}}:\;\|\theta\|_{L^{2}(0,T;L^{2}(\Omega))}\leq\widehat{C}_{1},\;\|\mathbf{h}\|_{[L^{2}(0,T)]^{N_{CT}}}\leq\widehat{C}_{2}\}$. Moreover, the corresponding $(\theta,\mathbf{v})\in W_{1}\times\mathbf{W}_{2}$ is a solution for the system (1) and (6) in the sense of Definition 5.

{@proof}

[Proof.] The existence is a direct consequence of the Schauder fixed point Theorem. Given an element $(\theta^{*},\mathbf{h}^{*})\in L^{2}(0,T;L^{2}(\Omega))\times[L^{2}(0,T)]^{N_{CT}}$, we have, thanks to (20), (39), (40), (42) and (43), the following estimates for $(\theta,\mathbf{h})=M(\theta^{*},\mathbf{h}^{*})$:

[TABLE]

where $C_{6}$, $C_{7}$, $C_{8}$ are positive constants that depend continuously on the coefficients and data. If we take the first inequality to the second and third ones, we obtain that:

[TABLE]

So, if we suppose that $\|\theta^{*}\|_{L^{2}(0,T;L^{2}(\Omega))}\leq\widehat{C}_{1}$ and $\|\mathbf{h}^{*}\|_{[L^{2}(0,T)]^{N_{CT}}}^{2}\leq\widehat{C}_{2}$, we have that

[TABLE]

Thus, we are led to solve the following inequality:

[TABLE]

However, it is obvious that, given $\widehat{C}_{1}$ and $\widehat{C}_{2}$, we can consider small enough data $\mathbf{v}^{0}$, $\mathbf{g}$, $\theta^{0}$, $\theta_{N}$, $\theta_{S}$ and $T_{r}$, such that

[TABLE]

Then, choosing suitable coefficients and data that verify (48)-(50), we have that $M$ maps elements of the set $\{(\theta^{*},\mathbf{h}^{*})\in L^{2}(0,T;L^{2}(\Omega))\times[L^{2}(0,T)]^{N_{CT}}:\;\|\theta^{*}\|_{L^{2}(0,T;L^{2}(\Omega))}\leq\widehat{C}_{1},\;\|\mathbf{h}^{*}\|_{[L^{2}(0,T)]^{N_{CT}}}\leq\widehat{C}_{2}\}$ into itself. Thus, thanks to Schauder fixed point Theorem, there exists a fixed point $(\theta,\mathbf{h})$ of operator $M$, such that the corresponding $(\theta,\mathbf{v})$ is a solution of the coupled system (1) and (6) in the sense of Definition 5. $\blacksquare$

5 Numerical resolution

Once proved in above section that the coupled system (1) and (6) admits a solution, we will introduce here a full numerical algorithm in order to compute it, and show several computational test for a realistic example.

We must recall here that our main aim is related to understanding which is the best strategy for reducing the water temperature in the upper layers. In order to achieve this objective, and for the sake of completeness, we will consider an algorithm able to deal with more general states than those we have presented in previous mathematical analysis of the problem. In particular, in the numerical resolution proposed here we will also take into account the possibility that $g^{k}(t)$ takes non-positive values. To be exact, if $g^{k}(t)>0$ we will say that the pump $k$ is turbinating (water enters by the collector $C^{k}$ and is turbinated by the corresponding pipeline to the injector $T^{k}$), and if $g^{k}(t)<0$ we will say that the pump $k$ is pumping (water enters by injector $T^{k}$ and is pumped to the collector $C^{k}$). As it is evident, the situation $g^{k}(t)=0$ corresponds to the case in which the pump $k$ is off. In addition, we will also suppose that the parameter $\epsilon$ used in the definition (2) tends to cero, that is, the mean temperature in the injectors is equal to the mean temperature in the collectors. It is essential emphasizing here that it is also possible to perform a similar mathematical analysis for the general case (with the obvious embarrassing notations), under the only assumption of the existence of a partition of the time interval verifying that the groups do not change their state within any element of the partition (as can be seen in following section devoted to the numerical examples).

5.1 Space-time discretization

For the discretization of the problem, let us consider a regular partition $0=t_{0}<t_{1}<\ldots<t_{N}=T$ of the time interval $[0,T]$ such that $t_{n+1}-t_{n}=\Delta t=\frac{1}{\alpha}$, $\forall n=0,\ldots,N-1$, and a family of meshes $\tau_{h}$ for the domain $\Omega$ with characteristic size $h$. Associated to this family of meshes, we also consider three compatible finite element spaces $Z_{h}$, $\mathbf{W}_{h}$ and $M_{h}$ corresponding, respectively, to the water temperature, velocity and the pressure of water. From the computational viewpoint, for the generation of the mesh associated to the domain and for the numerical resolution of the system, we propose the use of FreeFem++ [14]. Finally, we have employed an Uzawa algorithm [12] for computing the solution of the Stokes problems that appears after the discretization, and a fixed point algorithm for solving the nonlinearities.

So, we consider the following space-time discretization for system (1) and (6):

Dirichlet condition for the water velocity: We consider the following approximation for a function $\mathbf{g}=(g^{1},\ldots,g^{N_{CT}})\in[H^{1}(0,T)]^{N_{CT}}$,

[TABLE]

where, for all $n=0,\ldots,N$, $e^{n}\in\mathcal{C}([0,T])$ is such that $e_{|_{[t^{n-1},t_{n}]}}\in\mathcal{P}_{1}([t^{n-1},t^{n}])$, $n=1,\ldots,N$ and $e^{n}(t^{k})=\delta_{k,n}$, $k,n=0,\ldots,N$. It is well known that the linear closure of the functions of the basis $\{e^{n}\}_{n=0}^{N}$ is a vector subspace of $H^{1}(0,T)$, so if we suppose that $g^{k,0}=0$, $k=1,\ldots,N_{CT}$, we can consider the following coordinate vector in the basis of the corresponding subspace of $[H^{1}(0,T)]^{N_{CT}}$:

[TABLE]

with $g^{k,n}\in[-M,M]$, $k=1,\ldots,N_{CT}$, $n=1,\ldots,N$, with $M>0$ a technical bound related to mechanical characteristics of pumps. 2. 2.

Coupling of temperature in collectors and injectors: We denote by $\theta^{n}\in Z_{h}$ the water temperature at time step $n=0,\ldots,N$. Then, we can consider the following approximation in the case of $g^{k,n}>0$, $k=1,\ldots,N_{CT}$, $n=1,\ldots,N$, for (3) by functions $\gamma_{\theta}^{k}$, $k=1,\ldots,N_{CT}$, :

[TABLE]

Moreover, if we assume the value $\epsilon=\frac{\Delta t}{2}$ in the definition (4) of function $\rho_{\epsilon}$ we have that the support of $\rho_{\Delta t/2}(t^{n}-\frac{\Delta t}{2}-s)$ is contained in $(t^{n}-\Delta t,t^{n})=(t^{n-1},t^{n})$, for all $n=1,\ldots,N$, and then:

[TABLE]

Finally, we approximate each element $\varphi^{k}$ by the indicator function of the injector $T^{k}$, $k=1,\ldots,N_{CT}$, and each element $\tilde{\varphi}^{k}$ by the indicator function of the collector $C^{k}$, $k=1,\ldots,N_{CT}$. Thus, the temperature in each injector at time step $t_{n}$ is the mean temperature in the corresponding collector at time step $t_{n-1}$. The previous approximation is still valid in the general case with obvious modifications. For instance, if we suppose that $g^{k,n}<0$, $k=1,\ldots,N_{CT}$, $n=1,\ldots,N$, then, we consider the Dirichlet condition in the collectors:

[TABLE] 3. 3.

Water temperature: Given $\theta^{0}\in Z_{h}$, $\theta^{1}\in Z_{h}$ is the solution of:

[TABLE]

where the discrete characteristic $X^{0}(\mathbf{x})=\mathbf{x}-\Delta t\,\mathbf{v}^{0}(\mathbf{x})$.

Then, for each $n=1,\ldots,N$, $\theta^{n+1}\in Z_{h}$, with

[TABLE]

for all $k=1,\ldots,N_{CT}$, is the solution of:

[TABLE]

where the discrete characteristic $X^{n}(\mathbf{x})=\mathbf{x}-\Delta t\,\mathbf{v}^{n}(\mathbf{x})$, for $n=1,\dots,N$, and where the functional space $Z_{h}^{n}$ is given by:

[TABLE]

with $\textrm{sign}(y)$ denoting the sign function:

[TABLE]

For instance, in the case of $g^{k,n}>0$ for all $k=1,\ldots,N_{CT}$, $n=1,\ldots,N$, then $\left(\cup_{k=1}^{N_{CT}}\frac{1}{2}(1-\textrm{sign}(g^{k,n}))C^{k}\right)|\textrm{sign}(g^{k,n})|$ $\cup\left(\cup_{k=1}^{N_{CT}}\frac{1}{2}(1+\textrm{sign}(g^{k,n}))T^{k}\right)$ $=\cup_{k=1}^{N_{CT}}T^{k}$, so

[TABLE]

That is, we are considering the Dirichlet condition on the injectors. In the oposite case, $g^{k,n}<0$ for all $k=1,\ldots,N_{CT}$, $n=1,\ldots,N$, we have that

[TABLE]

and we are considering the Dirichlet condition on the collectors. In the general case, we can have alternating Dirichlet conditions on the collectors and injectors, so the definition for the space $Z_{h}^{n}$, $n=1,\ldots,N$, given in (58) covers all the possibilities. 4. 4.

Water velocity and pressure: Given $\mathbf{v}^{0}\in\mathbf{V}_{h}$, for each $n=0,1,\ldots,N-1$, the pair velocity/pressure $(\mathbf{v}^{n+1},p^{n+1})\in\mathbf{W}_{h}\times M_{h}$, with:

[TABLE]

is the solution of:

[TABLE]

where the functional space $\mathbf{V}_{h}=\{\mathbf{z}\in\mathbf{W}_{h}:\;\mathbf{z}_{|{\cup_{k=1}^{N_{CT}}(T^{k}\cup C^{k})}}=\mathbf{0}\}$.

Remark 10

It is worthwhile noting here that in above scheme we have to compute one additional time step for the water temperature. This shift is motivated by the dependence diagram shown in Fig. 2. We observe that, due to the time discretization proposed here, the Dirichlet boundary condition for the hydrodynamic model begins to have influence from the second step of time for the water temperature. In the first time step $n=1$, water circulation does not affect the water temperature, so for this first time step we consider $\Gamma_{N}=\partial\Omega\setminus\Gamma_{S}$. However, for $n\geq 2$, we impose at each of the $N_{CT}$ collector/injector pairs one of the boundary conditions (54), (55) or (56), depending on the sign of the Dirichlet condition $g^{k,n-1}$ in previous time step.

5.2 Numerical results

This final subsection is devoted to present some numerical results that we have obtained using realistic coefficients and data. Nevertheless, for the sake of clarity and comprehensibility, we will show results for a simplified 2D domain. For this purpose, we have considered a rectangular domain $\Omega=[0,16]\times[0,19]$ (measured in meters), corresponding a reservoir, in which we have distributed $N_{CT}=4$ collector/injector pairs with a symmetrical configuration similar to that shown in Fig. 1. For the time discretization we have chosen a time step of $\Delta t=1800$ seconds with $N=96$ time steps (which represents a time period of 2 days), and for the space discretization we have used a regular mesh formed by triangles of characteristic size $h=0.5$ meters (corresponding to $1599$ vertices). Finally, the finite element spaces employed for space discretizations have been the Taylor-Hood element $\mathcal{P}_{2}/\mathcal{P}_{1}$ for the hydrodynamic model, and the Lagrange $\mathcal{P}_{2}$ element for the water temperature. In Fig. 3 we can observe the evolution of the radiation temperature:

[TABLE]

along the whole period of 2 days ($17.28\ 10^{4}$ seconds), considering $a=0.1$, $R_{sw,net}=1000\,{\rm W}\,{\rm m}^{-2}$, $R_{lw,down}=350\,{\rm W}\,{\rm m}^{-2}$ (typical values in mediterranean countries during the summer) and multiplied by a sinusoidal function in order to simulate the effects of day and night. The parameters used for the numerical resolution of the coupled system can be seen in Table 1.

In order to analyze the influence of water artificial circulation in the thermal behavior of top $1.5$ meters from water upper layer we have solved the problem in five different scenarios:

NNNN: In this configuration we take $g^{k,n}=0$, for all $k=1,\ldots,N_{CT}$ and $n=1,\ldots,N$ (reference configuration with all the groups off). 2. 2.

TTTT: In this configuration we take $g^{k,n}=2.0\ 10^{-3}$, for all $k=1,\ldots,N_{CT}$ and $n=1,\ldots,N$ (all the groups are turbinating). 3. 3.

PPPP: In this configuration we take $g^{k,n}=-2.0\ 10^{-3}$, for all $k=1,\ldots,N_{CT}$ and $n=1,\ldots,N$ (all the groups are pumping). 4. 4.

TPTP: In this configuration we take $g^{1,n}=g^{3,n}=2.0\ 10^{-3}$ and $g^{2,n}=g^{4,n}=-2.0\ 10^{-3}$, for all $k=1,\ldots,N_{CT}$ and $n=1,\ldots,N$ (groups $1$ and $3$ are turbinating, and groups $2$ and $4$ are pumping). 5. 5.

PTPT: In this configuration we take $g^{1,n}=g^{3,n}=-2.0\ 10^{-3}$ and $g^{2,n}=g^{4,n}=2.0\ 10^{-3}$, for all $k=1,\ldots,N_{CT}$ and the groups $2$ and $4$ are turbinating).

In Fig. 4 we present the evolution of the mean temperature in the top $1.5$ meters upper layer along the whole time interval corresponding to two days. We can clearly distinguish here that the best configurations correspond, in a very evident manner, to the second scenario (TTTT) and to the fourth one (TPTP). Moreover, we can notice how third and fifth scenarios (PPPP and PTPT, respectively) do not improve in a significant way the reference configuration (NNNN).

Finally, we show in Fig. 5 water temperatures and velocities at last time step for NNNN and TTTT configurations and, in Fig. 6, the behavior of water at same time step for configurations TPTP and PTPT. (In all of the cases, velocities have been multiplied by an amplifying factor to make their graphic representations more perceptible). As we can easily notice, the best strategies correspond to evacuating the excess of temperature in the upper layers to the bottom layers instead of refrigerating the upper layers with cold water from the bottom ones.

Although we only present here one realistic example of application of our approach to understand the behaviour of water velocity and temperature in the upper section of the domain, we have developed many other numerical experiences for different choices of parameters and data (that will not be presented here for the sake of conciseness). However, from these computational tests we can derive two important consequences: For the modified Navier-Stokes equations, the second member -corresponding to the thermic term- shows less influence in the numerical resuls than the Smagorinsky turbulence term. For the convective heat equation, the radiation term in the nonlinear boundary condition affects in a significative way the final results. Finally, for a better resolution of the numerical examples, it would be possible to use, instead of an uniform mesh like the one employed in previous example, a finer mesh in the neighbourhoods of collectors and injectors. Nevertheless, this approach would mean a significant increase in the computational time, already quite high in the current case (especially in the part referring to the resolution of the hydrodynamic problem).

Acknowledgments

The authors thank the funding from project MTM2015-65570-P of MINECO/ FEDER (Spain).

Appendix A A radiation heat transfer problem with nonhomogeneous mixed boundary conditions

In this appendix we mathematically analyze a heat equation with an advective term and mixed boundary conditions of diffusive, convective and radiative type. The essential difficulty for demonstrating the existence and uniqueness of solution lies, in one hand, in the presence of nonlinear boundary conditions of radiative type and, in the other hand, the lack of regularity of the time derivative of the solutions. This lack of regularity does not allow us to take the solution itself as a test function in the variational formulation, and we are led to use more refined techniques (similar to those employed, for instance, in [1]).

So, we suppose that we have a convex domain $\Omega\subset\mathbb{R}^{3}$, whose boundary can be split into three disjoint, smooth enough parts: $\Gamma_{L}$, $\Gamma_{A}$ and $\Gamma_{R}$, with $\partial\Omega=\Gamma_{L}\cup\Gamma_{A}\cup\Gamma_{R}$. We denote by $\theta$ the solution of the following initial-boundary value problem:

[TABLE]

where $K>0\ ({\rm m}^{2}\,{\rm s}^{-1})$ is the thermal diffusivity, $b_{1}^{H}\geq 0\ ({\rm m}\,{\rm s}^{-1})$, for $H=A,R$, are the coefficients related to convective heat transfer through the boundaries $\Gamma_{A}$ and $\Gamma_{R}$, and $b_{2}^{R}>0\ ({\rm m}\,{\rm s}\,{\rm K}^{-3})$ is the coefficient related to radiative heat transfer through the boundary $\Gamma_{R}$, $\theta^{0}\geq 0\ ({\rm K})$ is the initial temperature, $\theta_{L}\geq 0\ ({\rm K})$ is Dirichlet temperature on $\Gamma_{L}$, $\varphi\geq 0\ ({\rm K})$ is the radiation temperature on $\Gamma_{R}$, and $\theta_{A},\,\theta_{R}\geq 0\ ({\rm K})$ are the temperatures related to convection heat transfer in surfaces $\Gamma_{A}$ and $\Gamma_{R}$.

Remark 11

In this work we will suppose that $\Gamma_{L}$, $\Gamma_{A}$ and $\Gamma_{R}$ are nonempty, but all the results can be easily extended to the case where $\Gamma_{A}$ and/or $\Gamma_{R}$ are empty sets. The only drawback is when $\Gamma_{L}=\emptyset$ because in this case we cannot use Poincare type inequalities, and we should apply another type of techniques for obtaining energy estimates in the Galerkin approximation. A related problem with $\Gamma_{L}=\emptyset$ was studied, for instance, in [20].

We consider the following spaces

[TABLE]

and

[TABLE]

Hypothesis 12

We will assume the following hypotheses for the coefficients and data:

$\theta^{0}\in L^{2}(\Omega)$. 2. 2.

$f\in L^{2}(0,T;L^{2}(\Omega))$. 3. 3.

$\mathbf{v}\in L^{10/3}(0,T;[L^{3}_{\sigma}(\Omega)]^{3})$. 4. 4.

$\theta_{A}\in L^{2}(0,T;L^{2}(\Gamma_{A}))$. 5. 5.

$\theta_{R}\in L^{2}(0,T;L^{2}(\Gamma_{R}))$. 6. 6.

$\varphi\in L^{5}(0,T;L^{5}(\Gamma_{R}))$. 7. 7.

$\theta_{L}(\mathbf{x},t)=g(t)\theta_{D}(\mathbf{x})$, a.e. $(\mathbf{x},t)\in\Gamma_{L}\times]0,T[$, with $g\in H^{1}(0,T)$ and $\theta_{D}={\widehat{\theta_{D}}}_{|_{\Gamma_{L}}}$, where $\widehat{\theta_{D}}\in H^{3/2}(\partial\Omega)$.

Definition 13

Within the framework established in Hypothesis 12, we say that an element $\theta\in W$ is a solution of the system (64) if there exists $\xi\in\widetilde{W}$ such that:

• •

$\theta=\zeta_{D}+\xi$, with $\zeta_{D}=g\,\beta_{0}(\widehat{\theta_{D}})\in W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))$, where $\beta_{0}$ is the right inverse of trace operator $\gamma_{0}:H^{2}(\Omega)\rightarrow H^{3/2}(\partial\Omega)$.

• •

$\xi(0)=\theta^{0}-\zeta_{D}(0)$, a.e. $\mathbf{x}\in\Omega$.

• •

$\xi\in\widetilde{W}$* is the solution of the following variational formulation:*

[TABLE]

where some of previous integrals must be understood as duality pairs, and

[TABLE]

We have the following result that we will prove in the following subsections:

Theorem 14

Within the framework established in Hypothesis 12, there exists a unique solution $\xi\in\widetilde{W}$ of equation (67) in the sense of Definition 13. Moreover, there exists a constant $C>0$, such that this solution satisfies the following inequalities:

[TABLE]

[TABLE]

In order to better understand the proof of the previous result, we will divide it in five parts: in the first part we will prove some technical results related to the space where we look for the solution. In the second part we will obtain the Galerkin approximation of problem (64). In the third part we will analyze the differential equation obtained from the Galerkin discretization. In the fourth part we will derive the convergence of the Galerkin approximation in a suitable space. Finally, in fifth part we will prove Theorem 14.

A.1 Part 1: Some technical results

We have the following lemma that we will use in following subsections:

Lemma 15

Let $Z\subset\tilde{X}$ a Banach space such that $L^{5}(\Gamma_{R})\subset Z^{\prime}$. The following inclusion is compact:

[TABLE]

{@proof}

[Proof.] The proof is a direct consequence of Aubin and Lions Lemma (see, for instance, Lemma 7.7 of [21]), the compactness of $\widetilde{X}$ in $L^{2}(\Gamma_{R})$ and the fact that $L^{4}(0,T;L^{4}(\Gamma_{R}))$ is is an interpolant between $L^{2}(0,T;L^{2}(\Gamma_{R}))$ and $L^{5}(0,T;L^{5}(\Gamma_{R}))$. Indeed, given a bounded sequence $\{\theta_{n}\}_{n\in\mathbb{N}}$ in $\{\theta\in W^{1,2,5/4}(0,T;\widetilde{X},\widetilde{X}^{\prime}):\theta_{|_{\Gamma_{R}}}\in L^{5}(0,T;L^{5}(\Gamma_{R}))\}$ we have that there exists a subsequence, that we will still denote in the same way, such that ${\theta_{n}}_{|_{\Gamma_{R}}}\rightarrow{\theta}_{|_{\Gamma_{R}}}$ in $L^{2}(0,T;L^{2}(\Gamma_{R}))$ and ${\theta_{n}}_{|_{\Gamma_{R}}}\rightharpoonup{\theta}_{|_{\Gamma_{R}}}$ in $L^{5}(0,T;L^{2}(\Gamma_{R}))$. We realize that ${\theta_{n}}_{|_{\Gamma_{R}}}\rightarrow\theta_{|_{\Gamma_{R}}}$ in $L^{4}(0,T;L^{4}(\Gamma_{R}))$ because:

[TABLE]

when $n\to\infty$. $\blacksquare$

Remark 16

We have that $W^{1,2,2}(0,T;H^{2}(\Omega),H^{2}(\Omega))\subset\mathcal{C}([0,T];H^{2}(\Omega))$ (in fact, it is well known that if $p,q\geq 1$ and $V_{1}\subset V_{2}$ continuously, then $W^{1,p,q}(0,T;$ $V_{1},V_{2})\subset\mathcal{C}([0,T];V_{2})$ continuously), and that $\mathcal{C}([0,T];H^{2}(\Omega))\subset W$. Then, the sum $\theta_{D}+\xi$ makes sense in the space $W$. Moreover, $\mathbf{v}\cdot\nabla\zeta_{D}\in L^{10/3}(0,T;L^{2}(\Omega))$.

We can also consider as a Dirichlet condition the restriction to $\Gamma_{L}$ of one element of the space $W^{1,2,2}(0,T;H^{s-1/2}(\partial\Omega),H^{s-5/2}(\partial\Omega))$, with $s\geq 2$. In this case, we can obtain an extension in the space $W^{1,2,2}(0,T;H^{s}(\Omega),H^{s-2}(\Omega))$ (cf. Theorem 3.2 of [11]) and, if we want to ensure that $W^{1,2,2}(0,T;H^{s}(\Omega),H^{s-2}(\Omega))$ $\subset W$, we can take, for instance, $s\geq 3$.

A.2 Part 2: Galerkin approximation

In this part we will construct a sequence of approximations that will converge to a solution of problem (64). So, let $\{\omega_{n}\}_{n\in\mathbb{N}}\subset\widetilde{X}$ be a dense subset of independent vectors of $\widetilde{X}$, such that ${\omega_{n}}_{|_{\Gamma_{R}}}\in L^{\infty}(\Gamma_{R})$, $\forall n\in\mathbb{N}$, which we can assume orthonormal in $L^{2}(\Omega)$. We also assume that the projection

[TABLE]

is selfadjoint and $\|{P_{N}}_{|_{Z}}\|_{\mathcal{L}(Z,Z)}\leq 1$, $\forall N\in\mathbb{N}$, where $Z$ is a Banach space, as given in Lemma 15. Then, for $N\in\mathbb{N}$, we denote by:

[TABLE]

where the coefficients $\xi_{n}^{N}(t)$, $n=1,\ldots,N$, are such that $\xi_{N}$ is the solution of the following differential equation:

[TABLE]

which can be rewritten in the following standard formulation:

[TABLE]

where:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $k=1,\ldots,N$, and $\xi_{0}^{N}=P_{N}(\theta^{0}-\zeta_{D}(0))$ is the projection of $\theta^{0}-\zeta_{D}(0)$ onto $\widetilde{X}_{N}=\langle\{\omega_{1},\ldots,\omega_{N}\}\rangle$.

Thus, we say than an element $\xi_{N}\in W^{1,5/4}(0,T;\widetilde{X}_{N})$ is a solution of system (75) if it satisfies the ordinary differential equation problem (76). We must recall here that $\|\xi_{0}^{N}\|_{L^{2}(\Omega)}\leq\|\theta^{0}-\zeta_{D}(0)\|_{L^{2}(\Omega)}$, $\forall N\in\mathbb{N}$, and that, if $\xi_{N}$ is a solution of problem (75), then $\mathbf{y}\in\mathcal{C}([0,T];\mathbb{R}^{N})$.

Lemma 17

Within the framework established in Hypothesis 12, there exits a constant $C>0$ independent of $N$ such that:

[TABLE]

[TABLE]

{@proof}

[Proof.] Multiplying (75) by $\xi_{k}^{N}(t)$, summing in $k$ and adding to both sides the term

[TABLE]

we have:

[TABLE]

As a consequence of the continuity of trace operator and of the inequalities of Young and Poincare, we obtain:

[TABLE]

where $\epsilon_{k}$, $k=1,\ldots,5$, are arbitrary strictly positive numbers, and $C_{k}$, $k=1,\ldots,5$, are constants that may depend on the trace operator and Young and Poincare’s inequalities. If we take $\epsilon_{k}$, $k=1,2,3$, such that $K-\sum_{k=1}^{3}\epsilon_{k}=K/2$ and $\epsilon_{5}=b_{2}^{R}/2$, then (renaming the constants if necessary):

[TABLE]

Adding to both sides $\frac{b^{S}_{2}}{2}\|\zeta_{D}(s)\|_{L^{5}(\Gamma_{R})}^{5}$, using the inequality $(a+b)^{5}\leq 16a^{5}+16b^{5}$ (for $a,b\geq 0$), and taking $\epsilon_{5}=b^{R}_{2}/64$:

[TABLE]

Finally, integrating over the time interval $[0,t]$ (renaming again the constants):

[TABLE]

where we have used that $\|\xi_{0}^{N}\|_{L^{2}(\Omega)}\leq\|\theta^{0}-\zeta_{D}(0)\|_{L^{2}(\Omega)}$, $\forall N\in\mathbb{N}$. Finally, applying Gronwall’s Lemma, we obtain that there exists a positive constant $C$ independent of $N$ such that (83) is satisfied.

For obtaining (84) it is sufficient to apply Holder inequality and bear in mind the fact that the projection operator $P_{N}$ onto $\widetilde{X}_{N}$ is bounded independently of $N$:

[TABLE]

$\blacksquare$

A.3 Part 3: Existence of solution for the Galerkin approximation

Now we will demonstrate that there exists, for each $N\in\mathbb{N}$, a unique absolutely continuous solution $\xi_{N}$ of equation (75).

Lemma 18

Within the framework established in Hypothesis 12, there exists a unique absolutely continuous solution defined on the whole time interval $[0,T]$ of Cauchy problem (76).

{@proof}

[Proof.] To state this lemma we can apply the Caratheodory theorem for ordinary differential equations (cf., for example, Theorem 5.2 of [13]). Indeed, $\mathbf{F}(\cdot,t)$ is continuous for any $t\in[0,T]$, and $\mathbf{F}(\mathbf{y},\cdot)\in L^{5/4}(0,T)$ for any $\mathbf{y}\in\mathbb{R}^{N}$. Then, given an open ball $B$ in $\mathbb{R}^{N}$, if we prove that there exist two functions $m_{B},\,l_{B}\in L^{1}(0,T)$ such that:

[TABLE]

we can conclude that problem (76) has a unique absolutely continuous solution, which can be extended to the boundary of $]0,T[\times B$. It is worthwhile mentioning here that, if

[TABLE]

then the solution $\mathbf{y}$ cannot reach the boundary of $B$ because of the a priori estimate (83) and the fact that $\|\xi_{N}(t)\|^{2}_{L^{2}(\Omega)}=\|\mathbf{y}(t)\|^{2}$. The first of the Caratheodory conditions (90) can be obtained by applying Holder inequality (in fact, we obtain that $m_{B}\in L^{5/4}(0,T)$). In the other hand, for the second Caratheodory condition, we can use the following inequality (straightforward consequence of the mean value Theorem):

[TABLE]

for $c=\lambda b+(1-\lambda)a$ and $\lambda\in(0,1)$. So, we obtain the following inequality:

[TABLE]

Therefore, we can conclude the existence of a function $l_{B}\in L^{1}(0,T)$ (in fact, $l_{B}\in L^{5/3}(0,T))$ such that the second of the Caratheodory conditions (90) is achieved. $\blacksquare$

A.4 Part 4: Convergence of the Galerkin approximation

In previous subsections we have seen that there exists a bounded sequence $\{\xi_{N}\}_{N\in\mathbb{N}}\subset\widetilde{W}$ of solutions of problem (76). In this subsection we will pass to the limit and obtain a solution of equation (67).

Lemma 19

There exists a subsequence of $\{\xi_{N}\}_{N\in\mathbb{N}}$, still denoted in the same way, such that:

$\displaystyle\xi_{N}\rightharpoonup\xi$* in $L^{2}(0,T;\widetilde{X})$,* 2. 2.

$\displaystyle\xi_{N}\rightharpoonup^{*}\xi$* in $L^{\infty}(0,T;L^{2}(\Omega))$,* 3. 3.

$\displaystyle\xi_{N}\rightarrow\xi$* in $L^{10/3-\epsilon}(0,T;L^{10/3-\epsilon}(\Omega))$,* 4. 4.

$\displaystyle\xi_{N}\rightarrow\xi$* in $L^{2}(0,T;L^{2}(\Gamma_{A}))$,* 5. 5.

$\displaystyle\xi_{N}\rightarrow\xi$* in $L^{4}(0,T;L^{4}(\Gamma_{R}))$**.*

{@proof}

[Proof.] Thanks to the boundedness in $\widetilde{W}$ of the sequence $\{\xi_{N}\}_{N\in\mathbb{N}}$, we obtain the first two convergences. The third and fourth limits are a direct consequence of Aubin, and Lions Lemma (cf. Lemma 7.7 of [21]) and the compactness of $H^{1}(\Omega)$, respectively, in $L^{6-\epsilon}(\Omega)$ and $L^{4-\epsilon}(\partial\Omega)$, $\forall\epsilon>0$. Finally, the fifth convergence is a consequence of Lemma 15. $\blacksquare$

Lemma 20

If $\{\xi_{N}\}_{N\in\mathbb{N}}$ is a bounded sequence in $\widetilde{W}$, then there exists a subsequence of $\{\xi_{N}\}_{N\in\mathbb{N}}$, still denoted in the same way, such that, for all $\phi\in L^{\infty}(0,T;L^{\infty}(\Gamma_{R}))$:

[TABLE]

{@proof}

[Proof.] Using the same technique that we have employed in the proof of Lemma 18, from the strong convergence of $\{\xi_{N}\}_{N\in\mathbb{N}}$ to $\xi$ in $L^{4}(0,T;L^{4}(\Gamma_{R}))$ we have:

[TABLE]

$\blacksquare$

A.5 Part 5: Proof of the main result of the Appendix

Now, we can demonstrate the Theorem 14:

{@proof}

[Proof.] We will divide the proof into three parts, in the first part we will pass to the limit in the Galerkin approximation in order to obtain a solution for the system (67), in the second part, we will derive the estimates (69) and (70) and, finally, in the third part we will prove the uniqueness of solution.

First, for a fixed index $k\in\mathbb{N}$, if we multiply (75) by a scalar function $\psi$ continuously differentiable on $[0,T]$, such that $\psi(T)=0$, integrate with respect to $t$, and integrate by parts, we have, $\forall N\geq k$:

[TABLE]

The passage to the limit for $N\to\infty$ in the integrals of the left-hand side is due to the Lemmas 19 and 20. We observe also that $\xi_{0}^{N}\to\theta^{0}-\zeta_{D}(0)$ in $L^{2}(\Omega)$. Hence, we find in the limit:

[TABLE]

for each $\eta\in\widetilde{X}$ which is a finite lineal combination of elements $\omega_{k}$. Since each term of above expression depends linearly and continuously on $\eta$, for the norm of $\widetilde{X}$, previous equality remains still valid, by continuity, for each $\eta\in\widetilde{X}$. Now, writing in particular (95) for $\psi=\phi\in\mathcal{D}(0,T)$, we obtain the variational formulation (67). Finally, we can prove that $\xi(0)=\theta^{0}-\zeta_{D}(0)$ multiplying (67) by the same $\psi$ as before, integrating by parts with respect to $t$, and comparing with (95).

Then, multiplying inequality (89) by $\phi$, with $\phi\in\mathcal{D}(0,T)$, $\phi\geq 0$, and integrating in $[0,T]$ we have:

[TABLE]

from which, taking into account that the norm of a reflexive Banach space is weakly lower semicontinuous, we can pass to the inferior limit thanks to convergences of Lemma 19:

[TABLE]

Thus, we obtain the following energy inequality for a.e. $t\in(0,T)$:

[TABLE]

Finally, (69) can be derived from above expression thanks to the Gronwall’s Lemma, and estimate (70) is a direct consequence of Holder inequality.

Now, we will prove the uniqueness of solution. Let us assume the existence of two solutions $\xi_{1}$ and $\xi_{2}$ for problem (67), and define $\xi_{12}=\xi_{1}-\xi_{2}$. We have that $\xi_{12}\in\widetilde{W}$, $\xi_{12}(0)=0$, a.e. $\mathbf{x}\in\Omega$, and that satisfies the following variational formulation:

[TABLE]

From a direct computation, we observe that $\frac{d\xi_{12}}{dt}=\alpha+\beta+\gamma$, where $\alpha\in L^{2}(0,T;\widetilde{X}^{\prime})$, $\beta\in L^{5/4}(0,T;L^{5/4}(\Gamma))$ and $\gamma\in L^{4/3}(0,T;L^{5/6}(\Omega))$, defined by:

[TABLE]

$\forall z\in\widetilde{X}$, \ \ \beta=-b_{2}^{R}\big{[}|{\xi_{1}}_{|_{\Gamma_{R}}}+{\zeta_{D}}_{|_{\Gamma_{R}}}|^{3}({\xi_{1}}_{|_{\Gamma_{R}}}+{\zeta_{D}}_{|_{\Gamma_{R}}})-|{\xi_{2}}_{|_{\Gamma_{R}}}+{\zeta_{D}}_{|_{\Gamma_{R}}}|^{3}({\xi_{2}}_{|_{\Gamma_{R}}}+{\zeta_{D}}_{|_{\Gamma_{R}}})\big{]}, and $\gamma=-\mathbf{v}\cdot\nabla\xi_{12}$.

Now, for each $\delta>0$, we define the following function:

[TABLE]

and its primitive,

[TABLE]

We can extend the results proved in Corollary 9 of [2] to our case (we just need to add the term in $L^{5/4}(0,T;L^{5/4}(\Gamma_{R}))$ in the proof of Lemma 6 of [2] and then we can use this result to prove Corollary 9) and then, we have:

[TABLE]

By a direct evaluation of previous expression, we have:

[TABLE]

a.e. $t\in]0,T[$. Moreover, it is obvious that:

[TABLE]

for $c(s)=\lambda(s)\xi_{1}(s)+(1-\lambda(s))\xi_{2}(s)$ and $\lambda(s)\in(0,1)$, a.e. $s\in]0,T[$. Thus, we can deduce that:

[TABLE]

Finally, from the definition (101) of $K_{\delta}$ it is straightforward that $0\leq|r|-K_{\delta}(r)\leq\frac{\delta}{2}$, $\forall r\in\mathbb{R}$, and then:

[TABLE]

with $|\Omega|$ denoting the volume of $\Omega$, which implies that $\xi_{12}(\mathbf{x},t)=0$, a.e. $(\mathbf{x},t)\in\Omega\times]0,T[$, and, consequently, $\xi_{1}=\xi_{2}$. $\blacksquare$