\documentclass[reqno]{amsart} \usepackage{graphicx} \usepackage{hyperref} \AtBeginDocument{{\noindent\small {\em Electronic Journal of Differential Equations}, Vol. 2006(2006), No. 30, pp. 1--17.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu (login: ftp)} \thanks{\copyright 2006 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2006/30\hfil Existence of weak solutions] {Existence of weak solutions for a non-homogeneous solidification \\ mathematical model} \author[M. M. Dur\'an, E. Ortega-Torres\hfil EJDE-2006/30\hfilneg] {Mario M. Dur\'an, Elva Ortega-Torres} % in alphabetical order \address{Mario Manuel Dur\'an\hfill\break Facultad de Ingenier\'{\i}a, Pontificia Universidad Cat\'olica de Chile, Casilla 306, Santiago 22, Chile} \email{mduran@ing.puc.cl} \address{Elva Ortega-Torres \hfill\break Departamento de Matem\'aticas, Universidad de Antofagasta, Casilla 170, Antofagasta, Chile} \email{eortega@uantof.cl} \date{} \thanks{Submitted April 21, 2005. Published March 16, 2006.} \thanks{This work supported by research grants Fondecyt\# 1040205 and ECOS/Conicyt\# C03-E08} \subjclass[2000]{35D05, 35Q35} \keywords{Weak solution; non-homogeneous solidification} \begin{abstract} This article studies the existence of weak solutions for a stationary non-homogeneous system of equations describing the solidification process of a binary alloy confined to a regular bounded domain in $\mathbb{R}^3$ and having mixed boundary conditions. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{definition}[theorem]{Definition} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} In a binary alloy an element, called {\it solute}, is dissolved in another substance, called {\it solvent}. The procedure for moulding the alloy consists of pouring the compound in its liquid state at a hight temperature in a mould. Once incandescent has been poured, it will loose heat and begin to solidify. At the beginning of the process a solidification front appears, which consists of solid material dendrites. Around these, the material in a liquid state adheres until reaching a completely solid phase. This process of solidification is fundamentally controlled by the thermodynamical variables of solute concentration and temperature in the mold. Thus the mathematical model takes into account the diffusion process, the heat transfer in a multi-component mixture and the transport of solute in the binary mixture. The reader is referred to \cite{hi-lo-ro} for the analysis of solidification models from the thermodynamic point of view. In \cite{ra-vo} solidification models are studied using the properties of thermodynamic averages. For a review of bifurcation phenomena in solidification models the reader is referred to \cite{ba-du}, whereas for numerical analysis of solidification models \cite{ah,pr-vo,co-le,ba-du-or,du-or-st} may be consulted. Additional analysis may be found in \cite{be-sa,don,luc-vi}. Preliminary studies regarding the existence and uniqueness of solutions for {\it stationary models} may be found \cite{bl-gas,bl-gas-ra,gai,gas,du-or-ra}. The aim of this article is to generalize the results established in the above-mentioned works to a less restricted model, resorting to Galerkin's technique and considering a non-homogeneous fluid. The outline of the paper is as follows. In Section 2 we state the mathematical model of the solidification process for a binary alloy. In Section 3 we state some notations, define the notion of weak solution, and prove a maximum principle. In Section 4 we study a regularized problem, define the notion of a regularized weak solution, and prove some a priori estimates. The proof of existence of a regularized weak solution is done in Section 5. Section 6 is devoted to proving the main existence theorem. \section{Model for binary alloy solidification} In this section we derive the mathematical model for studing the existence of weak solutions. Denoting by $\Omega$ the mould into which the melted matter was poured, we call $\Omega_l, \Omega_m$ and $\Omega_s$ the regions occupied by matter in liquid, mixture (coexistence of solid and liquid) and solid state, respectively. Furthermore, we identify the interfaces between these three sub-domains, calling $\Gamma_{ml}$ the mixture-liquid interface, and $\Gamma_{sm}$ the one for solid-mixture. With the purpose of considering boundary conditions, we set three regions of the exterior frontier $\Gamma= \partial \Omega$ which are denoted by $\Gamma_t, \Gamma_b$ and $\Gamma_v$ (see Figure 1). \begin{figure}[ht] \begin{center} \includegraphics[width=0.5\textwidth]{fig1} % fig11.eps \end{center} \caption{Solidification mould} \end{figure} During a binary alloy solidification process, the variables for solute concentration and temperature in the mould, determine the state of the matter through the phase diagram (see Figure 2) where $\theta_F$ and $\theta_E$ are the fusion and eutectic temperature respectively. As a fundamental reference we can point out, among others, the text-books by Milne-Thomson \cite{mi-th}, Shapiro \cite{sh}, the works by Blanc \& Gasser \cite{bl-gas}, Blanc \& al. \cite{bl-gas-ra} and the doctoral theses by Ahmad \cite{ah}, Gaillard \cite{gai} and Gasser \cite{gas}. To determine the state of the alloy, in the presence of thermodynamic equilibrium, we resort to the {\it phase diagram} (see Figure 2) which is formed by two regular curves: the {\it liquidus} $\gamma_\ell(\theta)$ and the {\it solidus} $\gamma_s(\theta)$. Then, we identify three zones: The liquid $\Theta_l$, the mixture $\Theta_m$ and the solid $\Theta_s$ which are defined by \begin{gather*} \Theta_l=\{(c, \theta) \in \mathbb{R}^2: 0 < c < c_{\max},\; \theta > \gamma_\ell^{-1}(c) \}\,, \\ \Theta_m=\{(c, \theta) \in \mathbb{R}^2 : \theta_E < \theta < \theta_F ,\; \gamma_s(\theta) < c < \gamma_\ell(\theta) \},\\ \Theta_s=\{(c, \theta) \in \mathbb{R}^2 : 0 < c < c_{\max},\; \theta < \mbox{Max}\{\gamma_s^{-1}(c), \theta_E\}\}\,, \end{gather*} where $c_{\max}=\gamma_\ell(\theta_E)$. Consequently the sub-domains $\Omega_l, \Omega_m$ and $\Omega_s$ can be described - using the phase diagram- by the thermodynamic variables for concentration and temperature, $(c,\theta)$. In fact, considering $\Omega$ an open domain subset of $\mathbb{R}^3$ with regular boundary $\Gamma$, we have \begin{gather*} \Omega_a =\{x \in \Omega : (c(x), \theta(x)) \in \Theta_a\}, \quad a = l, m, s, \\ \Omega_{ml} = \Omega_m \cup \Gamma_{ml} \cup \Omega_l. \end{gather*} \begin{figure}[ht] \begin{center} \includegraphics[width=0.6\textwidth]{fig2} % fig12.eps \end{center} \caption{Phase Diagram} \end{figure} In the solid-liquid mixture region, we define a function that will indicate the fraction of matter in solid state $f_s(c,\theta)$ and in liquid state $f_l (c, \theta) = 1-f_s(c,\theta)$, per unit of volume. Thus, a precise definition of $f_s(c, \theta)$ is given by (see Figure 2) \[ f_s(c,\theta)=\begin{cases} 0 & \mbox{if } (c,\theta) \in \Theta_l,\\ \frac{a} {a+b} & \mbox{if } (c,\theta) \in \Theta_m,\\ 1 & \mbox{if } (c,\theta) \in \Theta_s, \end{cases} \] where $a= \gamma_{\ell}(\theta)- c$ and $ b= c -\gamma_s(\theta)$. To model the transport phenomenon in the mixture region, we denote by $c_l$ and $c_s$ the concentration of solute in liquid and solid phases, which can be precisely determined through the {\it phase diagram} (see e.g. \cite{ah,ba-du-or,du-or-st,gai}). Thus, we deduce \begin{equation} \label{e2.1} c_l(c,\theta)=\begin{cases} c & \mbox{if } (c,\theta) \in \Theta_l,\\ \gamma_l(\theta) & \mbox{if } (c,\theta) \in \Theta_m,\\ 0 &\mbox{in }\{(c,\theta)\in\mathbb{R}^2 : 0 < c < \gamma_s(\theta_E),\; \theta < \gamma_s^{-1}(c)\}\,, \\ \gamma_{\ell}(\theta_E)&\mbox{in }\{(c,\theta)\in\mathbb{R}^2 : \gamma_s(\theta_E)\leq c\leq\gamma_{\ell}(\theta_E),\; \theta\leq\theta_E\}\,, \\ \gamma_{\ell}(\theta_E) & \mbox{in } \{(c,\theta)\in\mathbb{R}^2: c\ge\gamma_{\ell}(\theta_E)\}\,, \\ 0 & \mbox{in }\{(c,\theta)\in\mathbb{R}^2: c\leq 0\}\,. \end{cases} \end{equation} For a technical reason, to obtain the Maximum Principle, in \eqref{e2.1} we have modified the natural definition of the function $c_l(\cdot, \cdot)$. This fact does not change either the physical meaning or the mathematical sense of the model, since $f_s=1$ and $\mathbf{v}={\bf 0}$ on $\Theta_s$, and because $c_l(\cdot, \cdot)$ was extended outside $\overline{\Theta_l \cup \Theta_m}$ by continuity. In the sequel, the unknown function of concentration $c(\cdot)$ that intervenes in the equations is a weighted average between $c_l(\cdot, \cdot)$ and $c_s(\cdot, \cdot)$, that is, $ c= c_l (1-f_s)+c_s f_s$. It is clear that, as time passes, the zones $\Omega_a, a\in\{l,m,s\}$, will appear, evolve, and disappear within $\Omega$, which added to the force of gravity acting over the alloy, can develop a movement of velocity $\mathbf{v}$ and pressure field $p$ in the interior of the mixture $\Omega_{ml}$ (see \cite{ba-du-or} or \cite{gai-ra}). Then, using the incompressible Navier-Stokes equations coupled with terms of force (external and internal), and considering the diffusion processes of the concentration and the heat transfer in a multi-component mixture, the stationary nonhomogeneous solidification model for a binary alloy is given by (see \cite{ah,gai}) \begin{gather} - \mathop{\rm div} ( \eta (c, \theta) \nabla c )+\mathbf{v} \cdot \nabla c_l (c, \theta) = 0 \quad \text{in } \Omega, \label{e2.2}\\ \frac{\partial c}{\partial n} = 0 \quad \text{on } \Gamma, \label{e2.3}\\ \int_\Omega c (x) dx = c_g, \label{e2.4}\\ - \mathop{\rm div} ( \kappa (c, \theta) \nabla \theta )+\rho C_p \mathbf{v} \cdot \nabla \theta = 0 \quad \text{in } \Omega, \label{e2.5}\\ \theta = \theta_\delta \quad \text{on }\Gamma_b \cup \Gamma_t, \label{e2.6}\\ \frac{\partial \theta}{\partial n} = 0 \quad \text{on } \Gamma_v, \label{e2.7}\\ - 2 \mathop{\rm div} ( \nu (c, \theta) e (\mathbf{v}))+ \rho (\mathbf{v} \cdot \nabla )\mathbf{v} + F_i (c, \theta) \mathbf{v} +\nabla p = \mathbf{F}_e (c, \theta) \quad \text{in }\Omega_{ml}, \label{e2.8}\\ \mathop{\rm div} \mathbf{v} = 0 \quad \text{in }\Omega_{ml}, \label{e2.9}\\ \mathbf{v} = {\bf 0} \quad \text{in }\Omega_s, \label{e2.10}\\ \mathbf{v} = {\bf 0} \quad \text{on }\Gamma_b \cup \Gamma_t, \label{e2.11}\\ \mathbf{v} \cdot {\bf n} = 0 \quad \text{on }\Gamma_ v, \label{e2.12}\\ \sigma (\mathbf{v}, p) {\bf n} \wedge {\bf n} = {\bf 0} \quad \text{on }\Gamma_v, \label{e2.13} \end{gather} where $\eta(\cdot, \cdot), \kappa (\cdot, \cdot), \nu (\cdot, \cdot) $ are strictly positive real functions that represent {\it concentration diffusion}, {\it heat diffusion} and {\it dynamic viscosity}, respectively. $\rho$ is the {\it mean density} in the mixture considered constant, $C_p$ is the {\it caloric capacity}, and $ e (\mathbf{v})$ is the {\it linearized stress tensor} of $\mathbf{v}$. In equation \eqref{e2.4}, $ c_g$ denotes the total quantity of solute in the mixture and is a strictly positive real constant that satisfies the compatibility relation \begin{equation} \label{e2.14} 0 \leq c_g \leq |\Omega| \gamma_{\ell} (\theta_E). \end{equation} Furthermore, $ \theta_\delta$ is a known distribution of temperature on $\Gamma_b \cup \Gamma_ t $, ${\bf n}$ denotes the normal unit vector outward to the surface $\Gamma $ and $\sigma(\mathbf{v}, p)$ represents the stress tensor given by $$ \sigma(\mathbf{v}, p)=-p I + 2 \nu e (\mathbf{v})\quad \mbox{(Stokes' law).} $$ The external force $\mathbf{F}_e(\cdot,\cdot)$ which acts over the alloy is given by the Boussinesq approximation \begin{equation} \label{e2.15} \mathbf{F}_e(c, \theta)= \rho \mathbf{g} (\alpha (\theta-\theta_r)+ \beta (c_l(c, \theta)-c_r)), \end{equation} where $\alpha, \beta$ are known real constants, $\theta_r$ and $c_r$ are respectively a temperature and a concentration of reference, and $\mathbf{g}= (0, 0, - g)^t$ is the gravitational force. The term $F_i(\cdot,\cdot)$ represents the {\it internal force} within a porous material which opposes the movement of the fluid on $\Omega_{ml}$. It is known as the {\it Carman-Kozeny law}, and it is given by \begin{equation} F_i(c, \theta)= C_0 \frac{f_s^2(c, \theta)}{(1-f_s(c, \theta))^3} \quad \mbox{with } C_0>0. \end{equation} It is evident that in the liquid region $F_i(\cdot, \cdot) = 0$ and that $F_i(\cdot, \cdot)$ is not a $L^2(\Omega_{ml})$ function (this fact deals with a non classical problem, because we can not integrate the equation directly (2.8)). Note that the model \eqref{e2.2}-\eqref{e2.13} is more general than the one studied in \cite{bl-gas-ra} and \cite{du-or-ra}, where $\eta, \kappa, \nu$ constants were considered. Furthermore, in \cite{bl-gas-ra} the conditions $\theta=\theta_\delta$ and $\mathbf{v}={\bf 0}$ were considered on whole $\Gamma\ (\Gamma_v=\emptyset)$ which leads to a $H^2$-regular solution and a simpler model than the treated here. \section{Preliminaries, Variational Formulation and Maximum Principle} The vector-valued functions in $\mathbb{R}^3$ are denoted by $\{\mathbf{u}, \mathbf{v}, \mathbf{w}\}$, while the scalar functions are written $\{p, \varphi, \psi, \phi, c, \theta \}$. However, for the sake of simplicity we do not make any difference in denoting spaces of vector-valued or scalar-valued functions. We denote by $C$ a generic real constant. For all integer $m \geq 0$, and all real $1\leq p \leq \infty$, with $\|\cdot\|_{m,p,\Omega}$ we denote the standard norms of the scalar and vector-valued Sobolev spaces $W^{m,p}(\Omega)$ (see \cite{ad}). When there is no ambiguity about the domain $\Omega$, we simply write $ \|\cdot \|_{m,p}$. As usual, we will denote $ W^{m,2}(\Omega) \equiv H^m(\Omega)$ and $ W^{0,2}(\Omega) \equiv L^2(\Omega)$, which are Hilbert spaces with the usual inner product, and the standard norms are denoted by $\|\cdot\|_m$ and $ \| \cdot \|$, respectively. To study the existence of at least one weak solution of problem \eqref{e2.2}-\eqref{e2.13}, we defined the following function spaces: \begin{gather} H_\sigma(\Omega)=\{\mathbf{w} \in L^2(\Omega) : \mathop{\rm div} \mathbf{w}=0 \textrm{ in } \Omega,\; \mathbf{w} \cdot \mathbf{n}=0 \textrm{ on } \Gamma\}\,,\\ H_\theta(\Omega)=\{\psi \in H^1(\Omega) : \psi=0 \textrm{ on } \Gamma_b\cup\Gamma_t\}\,, \\ H_\mathbf{v}(\Omega)=\{\mathbf{w} \in H^1(\Omega) : \mathop{\rm div} \mathbf{w}=0 \textrm{ in } \Omega, \mathbf{w} \cdot \mathbf{n}=0 \textrm{ on } \Gamma_v \mbox{ and } \mathbf{w} =\mathbf{0} \textrm{ on }\Gamma_b\cup\Gamma_t\}\,. \end{gather} It is straightforward to prove that they are Hilbert spaces with the corresponding $L^2(\Omega)$ or $H^1(\Omega)$ norms (see \cite{Lions}). Moreover, it follows from Poincare's inequality that the gradient seminorm is, in fact, a norm equivalent to $H^1(\Omega)$-norm in the spaces $H_\theta(\Omega)$ and $H_\mathbf{v}(\Omega)$. Also, over the space $H_\mathbf{v}(\Omega)$, we consider the norm defined by \begin{equation} \label{e3.4} | e (\mathbf{v})|_{0,\Omega} = ( 2 \int_{\Omega} |e (\mathbf{v})|^2 dx ) ^{1/2} \end{equation} which by Korn's inequality, is a well known equivalent to the norm $\| \nabla \mathbf{v} \|$. Let $1 < p,q < \infty$ be such that $p^{-1}+q^{-1}=1$ (or $p=1$ and $q=\infty$), then for $\mathbf{u} \in L^p(\Omega), \mathbf{v} \in L^q(\Omega)$, we write as usual \[ (\mathbf{u},\mathbf{v})=\int_\Omega \mathbf{u}\cdot\mathbf{v} dx= \int_\Omega \sum_{i=1}^3 u_i v_i dx. \] We denote the classical bilinear and trilinear forms by \begin{equation} \label{e3.5} (\nabla \mathbf{u}, \nabla \mathbf{v}) = \sum_{i,j=1}^{3} \int_\Omega \frac{\partial u_j}{\partial x_i} \frac{\partial v_j}{\partial x_i} dx, \quad (\mathbf{u} \cdot \nabla \mathbf{v}, \mathbf{w}) = \sum_{i,j=1}^{3} \int_\Omega u_j \frac{\partial v_i}{\partial x_j} w_i dx, \end{equation} for all vector-valued functions $ \mathbf{u}, \mathbf{v}, \mathbf{w}$ for which the integrals make sense, and it is known that the trilinear form has the properties \begin{equation} \label{e3.6} (\mathbf{u} \cdot \nabla \mathbf{v}, \mathbf{v}) = 0, \quad (\mathbf{u} \cdot \nabla \mathbf{v}, \mathbf{w}) = - (\mathbf{u} \cdot \nabla \mathbf{w}, \mathbf{v}) \quad \forall \mathbf{u} \in H_\mathbf{v}, \ \forall \mathbf{v}, \mathbf{w} \in H^1(\Omega). \end{equation} Now, we state our problem rigorously establishing regularity assumptions on boundary $\Gamma $ and some technical assumptions on external data. \begin{itemize} \item[(S1)] $\Omega$ is a bounded, connected, and with non-empty interior set of $\mathbb{R}^3$; \item[(S2)] $\Gamma \in C^2$ such that $\mathop{\rm meas }(\Gamma_b)$ and $\mathop{\rm meas}(\Gamma_t)$ are strictly positive; \item[(S3)] The given external field of temperature $\theta_\delta$ belongs to $L^\infty(\Gamma_b\cup\Gamma_t)$ and has an extension in $H^1(\Omega)$; \item[S4)] The function of concentration in liquid phase $c_l(\cdot,\cdot)$ belongs to $W^{1,\infty}(\mathbb{R}^2)$ \item[(S5)] The functions $\eta(\cdot,\cdot), \kappa(\cdot, \cdot), \nu(\cdot, \cdot)$ belong to $C^0(\mathbb{R}^2)$ and satisfy \begin{gather*} \inf_{(c, \theta) \in \mathbb{R}^2} \{ \eta (c, \theta), \kappa (c, \theta), \nu (c, \theta) \} = \chi_0 > 0,\\ \sup_{(c, \theta) \in \mathbb{R}^2} \{ \eta (c, \theta), \kappa (c, \theta), \nu (c, \theta) \} = \chi_1 < \infty, \end{gather*} where $ \chi_0$ and $\chi_1$ are positive real constants. \end{itemize} In what follows, we define the meaning of a weak solution for \eqref{e2.2}-\eqref{e2.13} and state the main result. \begin{definition} \label{def3.1} \rm We say that the triplet $(c (x), \theta (x), \mathbf{v} (x)) \in H^1(\Omega)\times H^1(\Omega) \times H_\mathbf{v}(\Omega)$ is a weak solution of the stationary solidification problem \eqref{e2.2}-\eqref{e2.13} if it satisfies the variational formulation \begin{gather} (\eta (c,\theta)\nabla c, \nabla \varphi)+(\mathbf{v} \cdot \nabla c_l (c, \theta), \varphi)=0\,, \label{e3.7}\\ (\kappa (c, \theta)\nabla \theta, \nabla \psi)+\rho C_p (\mathbf{v}\cdot \nabla \theta, \psi )=0\,, \label{e3.8}\\ 2 (\nu (c,\theta) e (\mathbf{v}), e (\mathbf{w})) +\rho ((\mathbf{v} \cdot\nabla ) \mathbf{v}, \mathbf{w}) + (F_i (c, \theta) \mathbf{v}, \mathbf{w}) = (\mathbf{F}_e (c, \theta), \mathbf{w}), \label{e3.9} \end{gather} for all $\varphi\in H^1(\Omega), \psi\in H_\theta(\Omega), \mathbf{w}\in H_\mathbf{v}(\Omega)$ with $\mathop{\rm supp}\mathbf{w} \subset {\Omega_{ml}}$, and such that \begin{gather} \int_\Omega c (x) dx=c_g\,,\label{e3.10}\\ \theta-\theta_\delta \in H_\theta\,, \label{e3.11}\\ \mathbf{v} = \mathbf{0} \quad \ in \ \Omega_s\,. \label{e3.12} \end{gather} \end{definition} \begin{remark} \label{rm3.2} \rm It is not difficult to prove that any solution of the variational formulation \eqref{e3.7}-\eqref{e3.12} is also a solution of the problem \eqref{e2.2}-\eqref{e2.13} in the distribution sense, and so, both problems are equivalent. \end{remark} \begin{theorem} \label{thm3.3} Under the hypotheses {\rm (S1)--(S5)}, the problem \eqref{e2.2}-\eqref{e2.13} admits at least one weak solution. \end{theorem} We finish this section by proving that any weak solution of the system \eqref{e2.2}-\eqref{e2.13} satisfies some a priori estimates. \begin{proposition}[Maximum Principle] \label{prop3.4} Let $(c, \theta, \mathbf{v})$ be a weak solution of \eqref{e2.2}-\eqref{e2.13}. Then \begin{gather} 0\leq c (x) \leq \gamma_\ell(\theta_E) \quad \mbox{ a.e. in } \Omega, \label{e3.13}\\ \mathop{\rm ess\,inf}_{x \in \Gamma_b \cup \Gamma_t} \theta_\delta(x) \leq \theta (x) \leq \mathop{\rm ess\,sup} {x \in \Gamma_b \cup \Gamma_t} \theta_\delta(x) \quad \mbox{a.e. in } \Omega. \label{e3.14} \end{gather} \end{proposition} \begin{proof} Given a function $f$ defined in $\Omega$, we denote by $\Omega_+$ the set where $f$ is strictly positive and by $\Omega_-$ where it is not positive, then $\Omega=\Omega_+\cup\Omega_-$ and $\Omega_+\cap\Omega_-=\emptyset$. To show that the right side of inequality \eqref{e3.13} holds, we choose $\varphi=[c-\gamma_\ell(\theta_E)]_{+} \in H^1(\Omega)$ as a test function in \eqref{e3.7}, then by \eqref{e2.1} we have \[ (\eta (c, \theta) \nabla c, \nabla [c-\gamma_\ell(\theta_E)]_{+}) = -\int_{\Omega_{+}} \mathbf{v}(x) \cdot \nabla c_l (c,\theta)(c(x) -\gamma_\ell(\theta_E)) dx=0\,. \] A simple calculation gives \begin{align*} (\eta (c,\theta) \nabla c, \nabla [c-\gamma_\ell(\theta_E)]_{+}) &= (\eta (c,\theta)\nabla [c-\gamma_\ell(\theta_E)]_{+}, \nabla [c-\gamma_\ell(\theta_E)]_{+})\\ &= \|\eta^{1/2}(c, \theta) \nabla [c-\gamma_\ell(\theta_E)]_{+}\|^2\,, \end{align*} then, from (S5) and the above equalities, we obtain $\nabla [c-\gamma_\ell(\theta_E)]_{+} = 0$. Thus, $[c-\gamma_\ell(\theta_E)]_{+} = C $ is a non negative constant function in $\Omega$. By the definition of $[c-\gamma_\ell(\theta_E)]_{+}$, we have that $C =0$ or $C=c-\gamma_\ell(\theta_E)$. If $C =c-\gamma_\ell(\theta_E)$, then $c= C+ \gamma_\ell(\theta_E)\quad \textrm{a.e. in } \Omega$ and by using \eqref{e3.10}, we deduce \[ c_g=\int_\Omega c(x) dx = C |\Omega|+ \gamma_\ell(\theta_E) |\Omega| > \gamma_\ell(\theta_E) |\Omega|, \] which is in contradiction with the compatibility condition \eqref{e2.14}. Therefore, $C =0$ and the right side of \eqref{e3.13} has been proved. To verify the inequality of the left side of \eqref{e3.13}, we take $\varphi=[c]_-$ as a test function in \eqref{e3.7}. Then \[ ( \eta ( c, \theta) \nabla c, \nabla [c]_-) = -(\mathbf{v} \cdot \nabla c_l (c, \theta), [c]_-) = \int_{\Omega_-} \mathbf{v}(x) \cdot \nabla c_l(c(x), \theta(x))\ c(x) dx=0. \] As above, a simple calculation gives \[ ( \eta ( c, \theta) \nabla c, \nabla [c]_{-})= - ( \eta (c, \theta)\nabla [c]_{-},\nabla [c]_{-})= - \|\eta^{1/2}( c, \theta) \nabla [c]_{-}\|^2, \] thus, from (S5) and above inequalities, we deduce that $[c]_{-}= C $ is a non negative constant function in $\Omega$. By definition of $[c]_{-}$, we have that $C =0$ or $C= -c$. If $C= - c$, by integrating and using \eqref{e3.10}, we obtain \[ c_g=\int_\Omega c (x) dx = - C |\Omega| < 0\,, \] which implies $C =0$ and the left side of \eqref{e3.13} is obtained. To prove inequality \eqref{e3.14}, let us define $w=\mathop{\rm ess\,sup}\{\theta_\delta (x) : x \in \Gamma_b\cup\Gamma_t\}$ and consider $\psi=[\theta-w]_+ \in H_\theta(\Omega)$ as a test function in \eqref{e3.8}. Then \[ (\kappa (c,\theta )\nabla \theta, \nabla [\theta-w]_+) = - \rho C_p {\int_{\Omega_{+}}}\mathbf{v}(x) \cdot \nabla (\theta(x)-w)\ (\theta(x)-w) dx=0, \] and since $(\kappa (c,\theta)\nabla \theta, \nabla [\theta-w]_+)= \| \kappa^{1/2} (c, \theta) \nabla [\theta-w]_+\|^2$, we conclude that $[\theta -w]_+$ is a constant function in $\Omega$. From \eqref{e3.11} we get \[ \theta \leq \mathop{\rm ess\,sup} _{ x \in \Gamma_b \cup \Gamma_t} \theta_\delta(x) \text{ a.e. in } \Omega. \] To prove the statement of left-hand side of \eqref{e3.14}, we proceed in a similar way, defining $w =\mathop{\rm ess\,inf}\{\theta_\delta(x): x \in \Gamma_b\cup\Gamma_t\}$ and taking $[\theta-w]_- \in H_\theta(\Omega)$ as a test function in \eqref{e3.8}. In this way, we obtain \[ (\kappa (c, \theta) \nabla \theta,\nabla [\theta-w]_-) = - \rho C_p (\mathbf{v} \cdot \nabla \theta, [\theta-w]_-)=0, \] and since $(\kappa ( c, \theta)\nabla \theta,\nabla [\theta-w]_-)= -\| \kappa^{1/2}(c, \theta) \nabla [\theta-w]_-\|^2$, then we have $\nabla [\theta-w]_-= 0$ and consequently $[\theta-w]_-$ is a non negative constant function. Finally, we deduce that \[ \mathop{\rm ess\,inf}_{x \in \Gamma_b \cup \Gamma_t} \theta_\delta(x)\leq\theta \quad\textrm{ a.e. in } \Omega. \] \end{proof} \section{The Regularized Mathematical Problem} It is important to observe that due to its dependence on concentration c and temperature $\theta$, the set $\Omega_{ml}$ and the interfaces $\Gamma_{sm},\Gamma_{ml}$ are unknown a priori, and this is an additional difficulty of the problem. Moreover, as we said before the Carman-Kozeny term given in (2.16), becomes discontinuous when $x\in \Omega_s$ since $f_s(c(x), \theta(x))= 1$. Then, to avoid this singularity and to consider the system defined as a whole by the domain $\Omega$, we define a family of regularized problems, which are obtained by modifying the Carman-Kozeny term, as follows: \begin{equation} \label{e4.1} F_i^\epsilon (c_\epsilon, \theta_\epsilon)=C_0 \frac{f_s^2(c_\epsilon, \theta_\epsilon)}{(1-f_s(c_\epsilon, \theta_\epsilon)+\epsilon)^3}\quad\forall \epsilon \in (0,1]\,. \end{equation} More precisely, the family of regularized associated problems is defined as: Find the functions $ ( c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon , p_\epsilon ) : \Omega \to \mathbb{R}^6$, such that \begin{gather} - \mathop{\rm div} ( \eta (c_\epsilon, \theta_\epsilon) \nabla c_\epsilon )+\mathbf{v}_\epsilon \cdot \nabla c_l (c_\epsilon, \theta_\epsilon) = 0 \quad \text{in }\Omega, \label{e4.2}\\ \frac{\partial c_\epsilon}{\partial n} = 0 \quad \text{on }\Gamma, \label{e4.3}\\ \int_\Omega c_\epsilon (x) dx = c_g, \label{e4.4}\\ - \mathop{\rm div} ( \kappa (c_\epsilon, \theta_\epsilon) \nabla \theta_\epsilon ) +\rho C_p \mathbf{v}_\epsilon \cdot \nabla \theta_\epsilon = 0 \quad \text{in }\Omega, \label{e4.5}\\ \theta_\epsilon = \theta_\delta \quad \text{on }\Gamma_b \cup \Gamma_t, \label{e4.6}\\ \frac{\partial \theta_\epsilon}{\partial n} = 0 \quad \text{on }\Gamma_v, \label{e4.7}\\ - 2 \mathop{\rm div} ( \nu (c_\epsilon, \theta_\epsilon) e (\mathbf{v}_\epsilon))+ \rho (\mathbf{v}_\epsilon \cdot \nabla ) \mathbf{v}_\epsilon +F_i^\epsilon (c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon +\nabla p_\epsilon = \mathbf{F}_e (c_\epsilon, \theta_\epsilon) \quad \text{in }\Omega, \label{e4.8}\\ \mathop{\rm div} \mathbf{v}_\epsilon = 0 \quad \text{in }\Omega, \label{e4.9}\\ \mathbf{v}_\epsilon = {\bf 0} \quad \text{on }\Gamma_b \cup \Gamma_t, \label{e4.10}\\ \mathbf{v}_\epsilon \cdot {\bf n} = 0 \quad \text{on }\Gamma_ v, \label{e4.11}\\ \sigma (\mathbf{v}_\epsilon, p_\epsilon) {\bf n} \wedge {\bf n} = {\bf 0} \quad \text{on } \Gamma_v. \label{e4.12} \end{gather} In the same way as in Definition \ref{def3.1}, we give the following definition of a regularized weak solution. \begin{definition}[Regularized Weak Solution] \label{def4.1} We say that the triplet of functions $(c_\epsilon(x), \theta_\epsilon(x), \mathbf{v}_\epsilon(x))\in H^1(\Omega)\times H^1(\Omega) \times H_\mathbf{v} (\Omega)$ is a regularized weak solution of the stationary problem \eqref{e2.2}-\eqref{e2.13} if it satisfies the variational equations \begin{gather} (\eta (c_\epsilon, \theta_\epsilon) \nabla c_\epsilon, \nabla \varphi)+(\mathbf{v}_\epsilon \cdot \nabla c_l (c_\epsilon, \theta_\epsilon ), \varphi)= 0, \label{e4.13}\\ (\kappa (c_\epsilon, \theta_\epsilon) \nabla \theta_\epsilon, \nabla \psi)+ \rho C_p ( \mathbf{v}_\epsilon \cdot \nabla \theta_\epsilon, \psi ) = 0, \label{e4.14}\\ 2 ( \nu (c_\epsilon, \theta_\epsilon) e (\mathbf{v}_\epsilon), e (\mathbf{w})) + \rho ((\mathbf{v}_\epsilon \cdot \nabla ) \mathbf{v}_\epsilon, \mathbf{w}) +(F_i^\epsilon (c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon, \mathbf{w})= (\mathbf{F}_e (c_\epsilon, \theta_\epsilon), \mathbf{w}), \label{e4.15} \end{gather} for all $\varphi \in H^1(\Omega)$, $\psi \in H_\theta(\Omega)$, $\mathbf{w} \in H_\mathbf{v}(\Omega)$, and such that \begin{gather} \int_\Omega c_\epsilon (x)dx = c_g, \label{e4.16}\\ \theta_\epsilon - \theta_\delta \in H_\theta(\Omega). \label{e4.17} \end{gather} \end{definition} \begin{remark} \label{rm4.2} \rm It is not difficult to prove that any solution of the variational formulation \eqref{e4.13}-\eqref{e4.17} is a solution of the problem \eqref{e4.2}-\eqref{e4.12} in the distribution sense, and thus both problems are equivalent. \end{remark} In the sequel, we refer to the regularized weak solution of the stationary problem \eqref{e2.2}-\eqref{e2.13} solely as the weak solution of \eqref{e4.2}-\eqref{e4.12}. Also, in this case we have the Maximum Principle property for the regularized weak solution. \begin{proposition}[$\epsilon$ Maximum Principle] \label{prop4.3} The weak solution $( c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon )$ of problem \eqref{e4.2}-\eqref{e4.12}, satisfies the following inequalities \begin{gather} 0 \leq c_\epsilon (x) \leq \gamma_{\ell}(\theta_E) \quad \mbox{ a.e. in } \Omega,\label{e4.18}\\ \mathop{\mathop{\rm ess\,inf}}_{ x \in \Gamma_b \cup \Gamma_t} \theta_\delta (x) \leq \theta_\epsilon(x) \leq \mathop{\rm ess\, sup}_{ x \in \Gamma_b \cup \Gamma_t} \theta_\delta (x) \quad \mbox{ a.e. in } \Omega. \label{e4.19} \end{gather} \end{proposition} The proof of the above proposition is similar to the proof of Proposition \ref{prop3.4}, so we omit it. Using Proposition \ref{prop4.3} we can easily show the following a priori estimate. \begin{proposition} \label{prop4.4} Let $( c_\epsilon,\theta_\epsilon, \mathbf{v}_\epsilon )$ be a weak solution of the problem \eqref{e4.2}-\eqref{e4.12}. Then \begin{equation} \label{e4.20} \|\nabla c_\epsilon\|^2+\|\nabla \theta_\epsilon\|^2 +| e (\mathbf{v}_\epsilon)|_{0,\Omega}^2 \leq C, \end{equation} where $C$ does not depend on $\epsilon$. \end{proposition} \begin{proof} We deduce the following inequalities \begin{gather} \|\nabla c_\epsilon\| \leq C | e (\mathbf{v}_\epsilon)|_{0,\Omega}, \label{e4.21} \\ \|\nabla \theta_\epsilon\| \leq C (1+| e (\mathbf{v}_\epsilon)|_{0,\Omega} ), \label{e4.22}\\ | e (\mathbf{v}_\epsilon)|_{0,\Omega} \leq C\,, \label{e4.23} \end{gather} from which \eqref{e4.20} is an immediate consequence. To prove inequality \eqref{e4.21}, we take $\varphi=c_\epsilon$ as a test function in \eqref{e4.13}. Then, by using the H\"{o}lder inequality, \eqref{e3.6} and (S5), we have \[ \chi_0 \|\nabla c_\epsilon\|^2 \leq \|\eta^{1/2}(c_\epsilon, \theta_\epsilon) \nabla c_\epsilon\|^2 = ( \mathbf{v}_\epsilon \cdot \nabla c_\epsilon, c_l (c_\epsilon, \theta_\epsilon)) \leq C \|\mathbf{v}_\epsilon\| \|\nabla c_\epsilon\| \|c_l (c_\epsilon, \theta_\epsilon)\|_{0,\infty}. \] Thus, from (S4) and \eqref{e3.4}, the inequality \eqref{e4.21} follows. Since $\theta_\delta \in H^1 (\Omega)$, we put $\psi=\theta_\epsilon-\theta_\delta \in H_\theta(\Omega)$ as a test function in \eqref{e4.14}. Considering (S5), we have \[ \chi_0 \|\nabla \theta_\epsilon\|^2 \leq \| \kappa^{1/2}(c_\epsilon, \theta_\epsilon) \nabla \theta_\epsilon\|^2 = ( \kappa (c_\epsilon, \theta_\epsilon)\nabla \theta_\epsilon, \nabla \theta_\delta)+\rho C_p (\mathbf{v}_\epsilon \cdot \nabla \theta_\epsilon, \theta_\delta). \] Using the H\"older inequality, (S5) and the continuous Sobolev embedding $H^1(\Omega) \hookrightarrow L^4(\Omega)$, we obtain \begin{align*} \chi_0 \|\nabla \theta_\epsilon\|^2 &\leq \| \kappa (c_\epsilon,\theta_\epsilon)\|_{0,\infty} \|\nabla \theta_\epsilon\| \|\nabla \theta_\delta\| +C \|\mathbf{v}_\epsilon \|_{0,4} \|\nabla \theta_\epsilon\| \| \theta_\delta\|_{0,4}\\ &\leq \chi_1 \|\nabla \theta_\epsilon\| \|\nabla \theta_\delta\| + C \|\nabla \mathbf{v}_\epsilon\| \|\nabla \theta_\epsilon\| \|\theta_\delta\|_1\\ &\leq C \|\nabla \theta_\epsilon\| \|\theta_\delta\|_1 ( 1 + \|\nabla \mathbf{v}_\epsilon\| ) \end{align*} which by \eqref{e3.4}, is none other than inequality \eqref{e4.22}. Now, considering $\mathbf{w}=\mathbf{v}_\epsilon$ as a test function in \eqref{e4.15}, we have \begin{equation} \label{e4.24} \| \nu^{1/2} (c_\epsilon, \theta_\epsilon) e(\mathbf{v}_\epsilon)\|^2 + (F_i^\epsilon (c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon, \mathbf{v}_\epsilon) = (\mathbf{F}_e(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon). \end{equation} From the definition of $\mathbf{F}_e(\cdot, \cdot)$ given in \eqref{e2.15}, we obtain \begin{equation} \label{e4.25} \begin{aligned} |(\mathbf{F}_e(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon)| &\leq \rho \alpha |(\mathbf{g} \theta_\epsilon, \mathbf{v}_\epsilon)| + \rho | \alpha \theta_r + \beta c_r | |(\mathbf{g}, \mathbf{v}_\epsilon)| + \rho \beta |(\mathbf{g} c_l(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon)| \\ &\leq C \|\mathbf{g}\| \|\theta_\epsilon\|_{0,3} \| \mathbf{v}_\epsilon\|_{0,6} + C \|\mathbf{g}\| \| \mathbf{v}_\epsilon\|+ C \|\mathbf{g}\| \|c_l(c_\epsilon, \theta_\epsilon) \|_{0,\infty}\| \mathbf{v}_\epsilon\|\,. \end{aligned} \end{equation} Since $ c_l$ satisfies (S4), $\mathbf{g} \in L^\infty(\Omega)$ and $H^1(\Omega) \hookrightarrow L^6(\Omega)$ continuously, from \eqref{e4.25} we get \begin{equation} \label{e4.26} |(\mathbf{F}_e(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon)| \leq C \|\theta_\epsilon\|_{0,3} \| \nabla \mathbf{v}_\epsilon\|+ C \|\nabla \mathbf{v}_\epsilon\| \leq C ( 1+\|\theta_\epsilon\|_{0,3} )\|\nabla \mathbf{v}_\epsilon\|. \end{equation} substituting this inequality in \eqref{e4.24}, and taking into account \eqref{e3.4}, \eqref{e4.19}, (S5) and the fact that $(F_i^\epsilon (c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon, \mathbf{v}_\epsilon) \geq 0$, we get \[ \chi_0 | e (\mathbf{v}_\epsilon)|^2_{0,\Omega} \leq \| \nu^{1/2} (c_\epsilon, \theta_\epsilon) e (\mathbf{v}_\epsilon)\|^2 \leq C | e (\mathbf{v}_\epsilon)|_{0,\Omega}, \] which proves inequality \eqref{e4.23}. \end{proof} \section{Existence of Regularized Weak Solutions} In this section, we prove the existence of at least one weak solution of the problem \eqref{e4.2}-\eqref{e4.12}, by characterizing it as the limit of a sequence of approximated solutions defined on finite dimensional spaces. We state this result in what follows and the proof is done at the end of this section. \begin{theorem} \label{thm5.1} Under hypotheses {\rm (s1)--(S5)}, Problem \eqref{e4.2}--\eqref{e4.12} admits at least one weak solution. \end{theorem} If we set out $ \tilde \theta_\epsilon = \theta_\epsilon - \theta_\delta$, then to find a weak solution of the problem \eqref{e4.2}-\eqref{e4.12} becomes equivalent to: Find the triplet of functions $( c_\epsilon, \tilde \theta_\epsilon, \mathbf{v}_\epsilon ) \in H^1(\Omega)\times H_\theta(\Omega) \times H_\mathbf{v}(\Omega)$, such that \begin{gather} (\eta (c_\epsilon, \tilde \theta_\epsilon + \theta_\delta) \nabla c_\epsilon, \nabla \varphi)+(\mathbf{v}_\epsilon \cdot \nabla c_l (c_\epsilon, \tilde\theta_\epsilon + \theta_\delta), \varphi)= 0, \label{e5.1}\\ (\kappa (c_\epsilon, \tilde\theta_\epsilon + \theta_\delta) \nabla \tilde\theta_\epsilon, \nabla \psi)+ \rho C_p ( \mathbf{v}_\epsilon \cdot \nabla \tilde \theta_\epsilon, \psi )+ \rho C_p ( \mathbf{v}_\epsilon \cdot \nabla \theta_\delta, \psi ) \nonumber\\ = - (\kappa (c_\epsilon, \tilde\theta_\epsilon + \theta_\delta) \nabla \theta_\delta, \nabla \psi), \label{e5.2} \\ 2 ( \nu (c_\epsilon, \tilde\theta_\epsilon + \theta_\delta) e (\mathbf{v}_\epsilon), e (\mathbf{w}))+ \rho ((\mathbf{v}_\epsilon \cdot \nabla ) \mathbf{v}_\epsilon, \mathbf{w})+(F_i^\epsilon (c_\epsilon, \tilde \theta_\epsilon + \theta_\delta) \mathbf{v}_\epsilon, \mathbf{w}) \nonumber\\ = (\mathbf{F}_e (c_\epsilon, \tilde\theta_\epsilon + \theta_\delta), \mathbf{w}), \label{e5.3} \end{gather} for all $\varphi \in H^1(\Omega)$, $\psi \in H_\theta(\Omega)$, $\mathbf{w} \in H_\mathbf{v}(\Omega)$, and \begin{equation} \label{e5.4} \int_\Omega c_\epsilon (x)dx = c_g. \end{equation} Thus, if we find $( c_\epsilon, \tilde \theta_\epsilon, \mathbf{v}_\epsilon ) \in H^1(\Omega)\times H_\theta(\Omega) \times H_\mathbf{v}(\Omega)$ solution of \eqref{e5.1} - \eqref{e5.4}, we have that $( c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon )$ is a weak solution of \eqref{e4.2}-\eqref{e4.12}, where $ \theta_\epsilon = \tilde\theta_\epsilon + \theta_\delta$. To prove the existence of at least a weak solution to the stationary regularized problem \eqref{e4.2}-\eqref{e4.12}, we consider the Hilbert orthonormal basis $\{\varphi^i(x)\}^{\infty}_{i=1}$ of $H^1(\Omega)$, $\{\psi^i(x)\}^\infty_{i=1}$ of $H_\theta(\Omega)$ and $\{\mathbf{w}^i(x)\}^{\infty}_{i=1}$ of $H_\mathbf{v}(\Omega)$, whose elements could be chosen as solutions to the following spectral problems: \begin{gather*} (\nabla \varphi^i, \nabla c_\epsilon)= \alpha_i(\varphi^i, c_\epsilon) \quad \forall c_\epsilon \in H^1(\Omega), \\ (\nabla \psi^i,\nabla \tilde\theta_\epsilon)= \gamma_i(\psi^i, \tilde\theta_\epsilon) \quad \forall \tilde\theta_\epsilon \in H_\theta(\Omega), \\ (\nabla \mathbf{w}^i, \nabla \mathbf{v}_\epsilon)= \lambda_i (\mathbf{w}^i, \mathbf{v}_\epsilon) \quad \forall \mathbf{v}_\epsilon \in H_\mathbf{v}(\Omega). \end{gather*} Let $H^k$ be the subspace of $H^1(\Omega)$ spanned by $\{\varphi^1(x),\ldots,\varphi^k(x)\}$, $H_\theta^k$ be the subspace of $H_\theta(\Omega)$ spanned by $\{\psi^1(x),\ldots,\psi^k(x)\}$ and $H_\mathbf{v}^k$ be the subspace of $H_\mathbf{v}(\Omega)$ spanned by $\{\mathbf{w}^1(x),\ldots, \mathbf{w}^{3k}(x)\}$, respectively. For every $k \geq 1$, we define approximations $c^k_\epsilon,\tilde \theta^k_\epsilon$ and $\mathbf{v}^k_\epsilon$ of $c_\epsilon,\tilde \theta_\epsilon$ and $\mathbf{v}_\epsilon$ respectively, by means of the following finite expansions: \begin{equation} \label{e5.5} c^k_\epsilon(x) = \sum^k_{i=1} c_{ik}^\epsilon \varphi^i(x), \quad \tilde\theta^k_\epsilon(x) = \sum^k_{i=1} d_{ik}^\epsilon \psi^i(x), \quad \mathbf{v}^k_\epsilon(x) = \sum^{3k}_{i=1} e_{ik}^\epsilon \mathbf{w}^i(x), \end{equation} where the coefficients $c_{ik}^\epsilon$, $d_{ik}^\epsilon $ and $e_{ik}^\epsilon$ will be calculated in such way that $c^k_\epsilon, \tilde \theta^k_\epsilon$ and $\mathbf{v}^k_\epsilon$ solve the following equations, which are in fact an approximated system of \eqref{e5.1}-\eqref{e5.4}: \begin{gather} (\eta (c_\epsilon^k, \tilde \theta_\epsilon^k + \theta_\delta) \nabla c_\epsilon^k, \nabla \varphi^i)+(\mathbf{v}_\epsilon^k \cdot \nabla c_l (c_\epsilon^k, \tilde\theta_\epsilon^k + \theta_\delta), \varphi^i) = 0, \label{e5.6} \\ (\kappa (c_\epsilon^k, \tilde\theta_\epsilon^k + \theta_\delta) \nabla \tilde\theta_\epsilon^k, \nabla \psi^i)+ \rho C_p ( \mathbf{v}_\epsilon^k \cdot \nabla \tilde \theta_\epsilon^k, \psi^i )+ \rho C_p ( \mathbf{v}_\epsilon^k \cdot \nabla \theta_\delta, \psi^i ) \nonumber\\ = - (\kappa (c_\epsilon^k, \tilde\theta_\epsilon^k + \theta_\delta) \nabla \theta_\delta, \nabla \psi^i), \label{e5.7} \\ 2 ( \nu (c_\epsilon^k, \tilde\theta_\epsilon^k + \theta_\delta) e (\mathbf{v}_\epsilon^k), e (\mathbf{w}^i))+ \rho ((\mathbf{v}_\epsilon^k \cdot \nabla ) \mathbf{v}_\epsilon^k, \mathbf{w}^i)+(F_i^\epsilon (c_\epsilon^k, \tilde \theta_\epsilon^k + \theta_\delta) \mathbf{v}_\epsilon^k, \mathbf{w}^i) \nonumber\\ = (\mathbf{F}_e (c_\epsilon^k, \tilde\theta_\epsilon^k + \theta_\delta), \mathbf{w}^i), \label{e5.8} \end{gather} for all $\varphi^i \in H^k(\Omega)$, $\psi^i \in H_\theta^k$, $\mathbf{w}^i \in H_\mathbf{v}^k$, and \begin{equation} \label{e5.9} \int_\Omega c_\epsilon^k (x)dx = c_g. \end{equation} The proof of the solvability of the system \eqref{e5.6}-\eqref{e5.9} for any $k \in \mathbb{N}$, will be as done in \cite{du-fe-ro}, applying the Brouwer's Fixed Point theorem (see for instance \cite{eva}, \cite{gi-tr}, \cite{re-si}). In fact, we defined the operator \begin{align*} \Phi : H^k \times H^k_\theta \times H_\mathbf{v}^k & \to H^k \times H^k_\theta \times H_\mathbf{v}^k \\ (\xi, \phi, \mathbf{u}) & \mapsto \Phi(\xi, \phi, \mathbf{u})= (c, \tilde\theta, \mathbf{v}), \end{align*} where $(c, \tilde\theta, \mathbf{v})$ is the unique solution to the linearized problem \begin{gather} (\eta (\xi, \phi + \theta_\delta) \nabla c, \nabla \varphi^i)= - (\mathbf{v} \cdot \nabla c_l (\xi, \phi + \theta_\delta ), \varphi^i), \label{e5.10}\\ (\kappa (\xi, \phi + \theta_\delta ) \nabla \tilde \theta, \nabla \psi^i)+ \rho C_p ( \mathbf{u}\cdot \nabla \tilde \theta, \psi^i ) \nonumber\\ = -\rho C_p ( \mathbf{u} \cdot \nabla \theta_\delta, \psi^i ) - (\kappa (\xi, \phi + \theta_\delta ) \nabla \theta_\delta, \nabla \psi^i) , \label{e5.11}\\ 2 ( \nu (\xi, \phi + \theta_\delta) e (\mathbf{v}), e (\mathbf{w})) + \rho ((\mathbf{u} \cdot \nabla ) \mathbf{v}, \mathbf{w}^i) +(F_i (\xi, \phi + \theta_\delta ) \mathbf{v}, \mathbf{w}^i) \nonumber\\ = (\mathbf{F}_e (\xi, \tilde \theta + \theta_\delta), \mathbf{w}^i), \label{e5.12} \end{gather} for all $\varphi^i\in H^k$, $\psi^i\in H_\theta^k$, $\mathbf{w}^i\in H_\mathbf{v}^k$, and \begin{equation} \label{e5.13} \int_\Omega c (x) dx = c_g\,. \end{equation} From hypothesis (S4), we have $\mathbf{v}\cdot\nabla c_l (\xi,\phi+\theta_\delta)\in L^{3/2}(\Omega)$, so the right-hand side of \eqref{e5.10} defines a continuous linear operator in $H^k$ onto itself. As $\mathbf{u} \in L^6(\Omega), \nabla \theta_\delta \in L^2(\Omega)$ then $ \mathbf{u} \cdot\nabla \theta_\delta \in L^{3/2}(\Omega)$, and we conclude that the right-hand side of \eqref{e5.11} defines a continuous linear operator in $H^k_\theta$ onto itself. Since $\mathbf{F}_e (\xi, \tilde\theta +\theta_\delta)\in L^2(\Omega)$, we can conclude that the right-hand side of \eqref{e5.12} also defines a continuous linear operator in $H^k_\mathbf{v}$ onto itself (see \cite{bl-gas-ra}). The existence and uniqueness of $(c, \tilde\theta, \mathbf{v})$, the solution of \eqref{e5.10}-\eqref{e5.13}, follows directly from Lax-Milgram Lemma. In fact, given $(\xi, \phi, \mathbf{u})$ from \eqref{e5.11} we conclude the existence and uniqueness of function $\tilde\theta$ because the coercivity of the associated bilinear form. Next, from \eqref{e5.12} we deduce the existence and uniqueness of function $\mathbf{v}$, and finally, \eqref{e5.10} leads to the existence and uniqueness of function $c$. Thus the operator $\Phi$ is well defined and continuous from $H^k \times H^k_\theta \times H_\mathbf{v}^k$ onto itself. \begin{remark} \label{rmk5.2} \rm It is important to note that all concentration and temperature functions obtained as image of $\Phi$ operator satisfy the $\epsilon$-Maximal Principle property stated in Proposition \ref{prop4.3}. This implies $\|\tilde\theta + \theta_\delta\|_{0,3}\leq C$, where $C$ is a positive constant that does not depends on $k$ and $\epsilon$. \end{remark} \begin{proposition} \label{prop5.3} Assume the hypotheses of Theorem \ref{thm5.1}. Then the triplet $( c, \tilde\theta, \mathbf{v})$ defined by \eqref{e5.10}-\eqref{e5.13} satisfy the following estimates \begin{gather} \|\nabla c\| \leq R_c\,, \label{e5.14}\\ | e (\mathbf{v})|_{0,\Omega} \leq R_v\,,\label{e5.15}\\ \|\nabla \tilde \theta\| \leq R_\theta\,,\label{e5.16} \end{gather} where $R_c, R_v$ and $R_\theta$ are positive constants that do not depend on $k$ and $\epsilon$. \end{proposition} \begin{proof} Considering $\varphi^i=c$ and $\mathbf{w}^i=\mathbf{v}$ in \eqref{e5.10} and \eqref{e5.12} respectively, we have \begin{gather} \|\eta^{1/2} (\xi, \phi + \theta_\delta) \nabla c \|^2 = - (\mathbf{v} \cdot \nabla c_l (\xi, \phi + \theta_\delta ), c), \label{e5.17} \\ \| \nu^{1/2}(\xi, \phi + \theta_\delta) e (\mathbf{v})\|^2 +(F_i (\xi, \phi + \theta_\delta ) \mathbf{v}, \mathbf{v}) = (\mathbf{F}_e (\xi, \tilde \theta + \theta_\delta), \mathbf{v}). \label{e5.18} \end{gather} Since $F_i^\epsilon (\xi, \phi+\theta_\delta) \geq 0$ and considering (S5), from \eqref{e5.18} we have \begin{equation} \label{e5.19} \chi_0 | e (\mathbf{v})|^2_{0,\Omega} \leq |(\mathbf{F}_e (\xi, \tilde \theta + \theta_\delta), \mathbf{v})|. \end{equation} Similarly as done in \eqref{e4.26} and considering \eqref{e3.4}, from \eqref{e5.19} we can conclude \begin{equation} \label{e5.20} | e (\mathbf{v})|_{0,\Omega} \leq C (1+ \|\tilde\theta+\theta_\delta\|_{0,3} ). \end{equation} By Remark \ref{rmk5.2}, the function $\tilde\theta+ \theta_\delta$ satisfies the $\epsilon$ Maximum Principle property, then inequality \eqref{e5.20} implies that there exists a constant $R_v > 0$ such that \begin{equation} \label{e5.21} | e (\mathbf{v})|_{0,\Omega} \leq R_v. \end{equation} From \eqref{e5.17} and (S5), we get \begin{equation} \label{e5.22} \chi_0 \| \nabla c \|^2 \leq |(\mathbf{v} \cdot \nabla c_l (\xi, \phi + \theta_\delta ), c)| \leq C | e (\mathbf{v})|_{0,\Omega} \|\nabla c\| \|c_l (c , \theta)\|_{0,\infty}. \end{equation} Then, considering (S4), \eqref{e5.21} and \eqref{e5.22}, we obtain that there exists a constant $R_c > 0$ such that \begin{equation} \label{e5.23} \| \nabla c\| \leq R_c\,. \end{equation} Now, from \eqref{e5.11} just note that the function $\tilde\theta$ satisfies the following estimate \begin{equation} \label{e5.24} \chi_0 \|\nabla \tilde\theta\| \leq \rho C_p C_\Omega | e (\mathbf{u})|_{0,\Omega} \| \theta_\delta\|_{0,3}+ \chi_1 \|\nabla \theta_\delta\| \end{equation} with $\chi_0$ and $\chi_1$ given in (S5). Thus, if $| e (\mathbf{u})|_{0,\Omega}\leq R_v$, we deduce \begin{equation} \label{e5.25} \|\nabla \tilde \theta\| \leq R_\theta, \end{equation} where $R_\theta = \rho C_p C_\Omega \chi_0^{-1} R_v \| \theta_\delta\|_{0,3}+ \chi_0^{-1} \chi_1 \|\nabla \theta_\delta\|$. \end{proof} Now, taking into account Preposition 5.3, if we define the convex set \[ M= \{ \ ( \xi, \phi, \mathbf{u}) \in H^k \times H^k_\theta \times H^k_\mathbf{v}: \|\nabla \xi \|\leq R_c, \; \|\nabla \phi\| \leq R_\theta, \; | e (\mathbf{u})|_{0,\Omega} \leq R_v \}, \] it is straightforward to see that $\Phi(M) \subset M$ and thus we have proved the following results. \begin{lemma} \label{lem5.4} Under Hypotheses {\rm (S1)--(S5)} the operator $\Phi$ is well defined, continuous and possesses at least one fixed point in the set $M$. \end{lemma} \begin{corollary}[Existence of discrete regularized weak solutions] \label{coro5.5} Under assumptions (S1)--(S5), system \eqref{e5.6}-\eqref{e5.9} admits at least one solution $(c^k_\epsilon, \tilde\theta^k_\epsilon, \mathbf{v}_\epsilon^k)$, for all $\epsilon>0$, and all $k\geq 1$. \end{corollary} We conclude this section by proving its main theorem. \begin{proof}[Proof of Theorem \ref{thm5.1}] From Proposition \ref{prop5.3} and since the embedding $H^1(\Omega)\hookrightarrow L^2(\Omega)$ is compact, there exists $(c_\epsilon (x),\tilde\theta_\epsilon (x),\mathbf{v}_\epsilon (x))$ such that as $k\to \infty$ \begin{gather} c^k_\epsilon \to c_\epsilon \quad \mbox{weakly in }H^1(\Omega) \mbox{ and strongly in } L^2(\Omega), \label{e5.26}\\ \tilde\theta^k_\epsilon \to \tilde\theta_\epsilon \quad \mbox{weakly in }H_\theta(\Omega)\mbox{ and strongly in } L^2(\Omega), \label{e5.27}\\ \mathbf{v}^k_\epsilon \to \mathbf{v}_\epsilon \quad \mbox{ weakly in }H_\mathbf{v}(\Omega) \mbox{ and strongly in } H_\sigma(\Omega). \label{e5.28} \end{gather} We shall prove that the triplet $(c_\epsilon (x), \tilde \theta_\epsilon (x),\mathbf{v}_\epsilon (x))$ satisfies system \eqref{e5.1}-\eqref{e5.4}, which is in fact equivalent to showing that $(c_\epsilon (x),\theta_\epsilon (x),\mathbf{v}_\epsilon (x))$ is a weak solution of the regularized problem \eqref{e4.2}-\eqref{e4.12}, where $\theta_\epsilon = \tilde \theta_\epsilon + \theta_\delta$. We observe that if $\zeta^k_\epsilon(\cdot, \cdot), \zeta_\epsilon(\cdot, \cdot)$ satisfies (S5), \begin{align*} |(\zeta^k_\epsilon \nabla b_\epsilon^k - \zeta_\epsilon \nabla b_\epsilon, \nabla \xi^i)| &\leq |(\zeta^k_\epsilon ( \nabla b_\epsilon^k- \nabla b_\epsilon), \nabla \xi^i)| + |((\zeta^k_\epsilon - \zeta_\epsilon) \nabla b_\epsilon, \nabla \xi^i)|\,,\\ &\leq \chi_1 |(\nabla b_\epsilon^k - \nabla b_\epsilon, \nabla \xi^i)|+ \| \zeta_\epsilon^k - \zeta_\epsilon \|_{0,\infty} \|\nabla b_\epsilon \| \|\nabla \xi^i \|\,, \end{align*} then, by Proposition \ref{prop4.4}, we have \begin{equation} \label{e5.29} |(\zeta^k_\epsilon \nabla b_\epsilon^k - \zeta_\epsilon \nabla b_\epsilon, \nabla \xi^i)|\leq \chi_1 | (\nabla b_\epsilon^k - \nabla b_\epsilon, \nabla \xi^i)| + C \| \zeta_\epsilon^k - \zeta_\epsilon\|_{0,\infty}. \end{equation} Therefore, considering (S5) and \eqref{e5.26}-(5.27) by the weak convergence, from (5.29) we have as $k \to \infty$ \begin{gather} (\eta (c_\epsilon^k, \tilde \theta_\epsilon^k + \theta_\delta) \nabla c_\epsilon^k - \eta (c_\epsilon, \tilde \theta_\epsilon + \theta_\delta) \nabla c_\epsilon, \nabla \varphi^i) \to 0, \label{e5.30}\\ ( \kappa ( c_\epsilon^k, \tilde \theta_\epsilon^k + \theta_\delta) \nabla \tilde\theta^k_\epsilon - \kappa ( c_\epsilon, \tilde \theta_\epsilon + \theta_\delta)\nabla \tilde\theta_\epsilon, \nabla \psi^i) \to 0. \label{e5.31} \end{gather} Similarly as \eqref{e5.29}, considering (S5) and \eqref{e2.8}, as $k \to \infty$, we obtain \begin{equation} \label{e5.32} (\nu ( c_\epsilon^k, \tilde \theta_\epsilon^k + \theta_\delta) e (\mathbf{v}^k_\epsilon) - \nu ( c_\epsilon, \tilde \theta_\epsilon + \theta_\delta) e (\mathbf{v}_\epsilon), e (\mathbf{w}^i)) \to 0. \end{equation} In non linear terms, we have \begin{align*} |(\mathbf{v}^k_\epsilon \cdot \nabla b^k_\epsilon -\mathbf{v}_\epsilon \cdot \nabla b_\epsilon , \xi^i)| &=|-(\mathbf{v}^k_\epsilon \cdot \nabla \xi^i, b^k_\epsilon )+(\mathbf{v}_\epsilon \cdot \nabla \xi^i,b_\epsilon )|\\ &=|((\mathbf{v}_\epsilon -\mathbf{v}^k_\epsilon )\cdot \nabla \xi^i, b^k_\epsilon )-(\mathbf{v}_\epsilon \cdot \nabla \xi^i,(b^k_\epsilon -b_\epsilon ))|\\ &\leq \|\mathbf{v}^k_\epsilon -\mathbf{v}_\epsilon \| \|\nabla \xi^i\|_{0, 3} \|b^k_\epsilon \|_{0, 6} +\|\mathbf{v}_\epsilon \|_{0, 6}\|\nabla \xi^i\|_{0, 3} \|b^k_\epsilon - b_\epsilon \|\\ &\leq C (\|\mathbf{v}^k_\epsilon -\mathbf{v}_\epsilon \| \|\nabla b^k_\epsilon \| + \|\nabla \mathbf{v}_\epsilon \| \|b^k_\epsilon -b_\epsilon \|) \|\nabla \xi^i\|_{0, 3}\\ &\leq C \|\mathbf{v}^k_\epsilon -\mathbf{v}_\epsilon \| \|\nabla b^k_\epsilon \| +C \|\nabla \mathbf{v}_\epsilon \| \|b^k_\epsilon -b_\epsilon \|, \end{align*} which together \eqref{e5.26}-\eqref{e5.28}, implies \begin{gather} (\mathbf{v}^k_\epsilon \cdot \nabla c_l(c_\epsilon^k, \tilde \theta^k_\epsilon +\theta_\delta) -\mathbf{v}_\epsilon \cdot \nabla c_l (c_\epsilon, \tilde \theta_\epsilon+\theta_\delta), \varphi^i) \to 0, \label{e5.33}\\ (\mathbf{v}^k_\epsilon \cdot \nabla \tilde\theta^k_\epsilon - \mathbf{v}_\epsilon \cdot \nabla \tilde\theta_\epsilon\,, \psi^i) \to 0, \label{e5.34} \\ (\mathbf{v}^k_\epsilon \cdot \nabla \mathbf{v}^k_\epsilon - \mathbf{v}_\epsilon \cdot \nabla \mathbf{v}_\epsilon\,, \mathbf{w}^i) \to 0. \label{e5.35} \end{gather} For the external force term, we observe that \[ \mathbf{F}_e(c_\epsilon^k,\tilde\theta_ \epsilon^k+\theta_\delta)- \mathbf{F}_e(c_\epsilon, \tilde\theta_ \epsilon+\theta_\delta) = \rho \mathbf{g} [ \alpha (\tilde \theta_ \epsilon^k-\tilde \theta_\epsilon) + \beta (c_l(c_\epsilon^k,\tilde\theta_\epsilon^k + \theta_\delta) -c_l( c_\epsilon,\tilde\theta_\epsilon + \theta_\delta))], \] so \begin{align*} &\|\mathbf{F}_e(c_\epsilon^k, \tilde\theta_\epsilon^k+\theta_\delta)- \mathbf{F}_e(c_\epsilon, \tilde\theta_\epsilon+\theta_\delta)\| \\ &\leq C \|\mathbf{g}\|_\infty (\|\tilde \theta_ \epsilon^k- \tilde \theta_\epsilon \| + \|c_l(c_\epsilon^k,\tilde\theta_\epsilon^k + \theta_\delta) -c_l( c_\epsilon,\tilde\theta_\epsilon + \theta_\delta)\|), \end{align*} then, by the strong convergence in $L^2(\Omega)$ and the continuity of $c_l(\cdot, \cdot)$, we can conclude that as $k \to \infty$, \begin{equation} \label{e5.36} (\mathbf{F}_e(c_\epsilon^k, \tilde\theta_\epsilon^k+\theta_\delta)- \mathbf{F}_e(c_\epsilon, \tilde\theta_\epsilon+\theta_\delta), \mathbf{w}^i) \to 0. \end{equation} For the internal force term $F_i^\epsilon (\cdot,\cdot)$, we have \begin{align*} &F_i^\epsilon(c_\epsilon^k, \tilde\theta_ \epsilon^k+\theta_\delta) \mathbf{v}_\epsilon^k - F_i^\epsilon(c_\epsilon, \tilde\theta_\epsilon+ \theta_\delta)\mathbf{v}_\epsilon \\ &=(F_i^\epsilon(c_\epsilon^k, \tilde \theta_\epsilon^k +\theta_\delta)-F_i^\epsilon(c_\epsilon, \tilde \theta_\epsilon+\theta_\delta))\mathbf{v}_\epsilon^k + F_i^\epsilon(c_\epsilon, \tilde\theta_ \epsilon+\theta_\delta) (\mathbf{v}_\epsilon^k-\mathbf{v}_\epsilon). \end{align*} Then \begin{equation} \label{e5.37} \begin{aligned} &|(F_i^\epsilon(c_\epsilon^k, \tilde\theta_ \epsilon^k+\theta_\delta) \mathbf{v}_\epsilon^k -F_i^\epsilon (c_\epsilon, \tilde\theta_ \epsilon +\theta_\delta)\mathbf{v}_\epsilon, \mathbf{w}^i)| \\ & \leq C \|F_i^\epsilon(c_\epsilon^k, \tilde\theta_\epsilon^k +\theta_\delta)-F_i^\epsilon(c_\epsilon, \tilde\theta_\epsilon+\theta_\delta) \|_{0, \infty}\|\mathbf{v}_\epsilon^k \|\\ &\quad +C \|F_i^\epsilon(c_\epsilon, \tilde\theta_\epsilon+\theta_\delta) \|_{0, \infty} \|\mathbf{v}_\epsilon^k -\mathbf{v}_\epsilon \|\,. \end{aligned} \end{equation} By the regularity of $F_i^\epsilon(\cdot, \cdot)$, together \eqref{e5.28} and \eqref{e5.37} and the continuity of $c_l(\cdot, \cdot)$, we deduce as $k\to \infty$, \begin{equation} \label{e5.38} (F_i^\epsilon(c_\epsilon^k, \tilde\theta_ \epsilon^k+\theta_\delta) \mathbf{v}_\epsilon^k -F_i^\epsilon (c_\epsilon, \tilde\theta_ \epsilon +\theta_\delta)\mathbf{v}_\epsilon, \mathbf{w}^i) \to 0. \end{equation} Therefore, from \eqref{e5.30}-\eqref{e5.36} and \eqref{e5.38}, we have that $(c_\epsilon, \tilde\theta_ \epsilon, \mathbf{v}_\epsilon)$ satisfy \eqref{e5.1}-\eqref{e5.4} and consequently $(c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon)$ is a weak solution of \eqref{e4.2}-\eqref{e4.12}, with $\theta_\epsilon = \tilde\theta_\epsilon+\theta_\delta $. \end{proof} \section{Existence of a Weak Solution to the Stationary Problem} In this section, we prove the main result of this paper. \begin{proof}[Proof of Theorem \ref{thm3.3}] First we show the existence of a weak solution. From Proposition \ref{prop4.4}, for $\epsilon \in (0,1]$, any regularized weak solution $(c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon)$ of the problem \eqref{e2.2}-\eqref{e2.13} is uniformly bounded in $H^1(\Omega)^2 \times H_\mathbf{v}$ for a constant independent of $\epsilon$. Then, since the Sobolev embedding $H^1(\Omega)\hookrightarrow L^2(\Omega)$ is compact, we have that there exists a function $(c, \theta, \mathbf{v})\in H^1(\Omega)^2\times H_\mathbf{v}(\Omega)$ and a subsequence, still indexed by $\epsilon$, such that as $\epsilon\to 0$ \begin{gather} c_\epsilon \to c \quad \textrm{ weakly in } H^1(\Omega) \textrm{ and strongly in } L^2(\Omega), \label{e6.1}\\ \theta_\epsilon \to \theta \quad \textrm{ weakly in } H^1(\Omega) \textrm{ and strongly in } L^2(\Omega), \label{e6.2} \\ \mathbf{v}_\epsilon \to \mathbf{v} \quad \textrm{ weakly in } H^1(\Omega) \textrm{ and strongly in } L^2(\Omega). \label{e6.3} \end{gather} The proof that $(c, \theta, \mathbf{v})$ is a weak solution to the problem \eqref{e2.2}-\eqref{e2.13} is similar as we prove that $(c_\epsilon, \theta_\epsilon, \mathbf{v}_\epsilon)$ is a weak solution to the regularized problem. We must only take care in the limit process with respect to internal force term. We note that regions $ \Omega_s, \Omega_m$ and $ \Omega_l$ are defined by limit functions $c $ and $\theta$, as given in section 2. Let $\mathbf{w} \in C_0^\infty$ be with compact support $K \subset\Omega_{ml}$. Then, for an certain $\delta >0$ and for every $x \in K$, \begin{equation} \label{e6.4} f_s(c(x),\theta(x))<1-\delta, \end{equation} which implies \begin{equation} \label{e6.5} \|F_i(c, \theta)\|_{0,\infty,K} 0$ such that for all $\epsilon \in (0,\epsilon_\delta)$ and $\forall x \in K$, \begin{equation} |f_s(c_\epsilon(x),\theta_\epsilon(x))-f_s(c(x), \theta(x))|<\frac{\delta}{2}\,. \end{equation} Thus, from the last inequality and \eqref{e6.4}, we have \begin{equation} \label{e6.7} f_s(c_\epsilon(x), \theta_\epsilon(x)) < \frac{\delta}{2}+f_s(c(x), \theta(x)) < \frac{\delta}{2}+ ( 1 - \delta) < 1-\frac{\delta}{2} \quad \forall x \in K. \end{equation} Consequently, when $\epsilon \to 0$, \begin{equation} \label{e6.6} F_i^\epsilon(c_\epsilon (x), \theta_\epsilon(x)) \to F_i(c(x),\theta (x)) \quad \textrm{ in } C^0(K). \end{equation} Similarly as in \eqref{e5.37}, using the H\"older inequality and \eqref{e6.5} we obtain \begin{equation} |(F_i^\epsilon(c_\epsilon, \theta_\epsilon ) \mathbf{v}_\epsilon - F_i (c, \theta) \mathbf{v}, \mathbf{w} ) | \leq C \|F_i^\epsilon (c_\epsilon, \theta_\epsilon )-F_i(c,\theta)\|_{0, \infty,K} +C(\delta) \|\mathbf{v}_\epsilon - \mathbf{v}\|_{0,2,K}. \end{equation} Therefore, using \eqref{e6.3} and \eqref{e6.6} in \eqref{e6.7} we deduce \begin{equation} \label{e6.8} (F_i^\epsilon(c_\epsilon, \theta_\epsilon ) \mathbf{v}_\epsilon - F_i (c, \theta) \mathbf{v}, \mathbf{w} ) \to 0 \quad \mbox{ as } \epsilon \to 0, \end{equation} for all $\mathbf{w} \in C_0^\infty$ with compact support $K \subset \Omega_{ml}$. It follows by density arguments and the definition of $H_\mathbf{v} (\Omega_{ml})$, that equation \eqref{e3.9} holds for any $\mathbf{w} \in H_\mathbf{v}(\Omega)$ with $\mathop{\rm supp}\mathbf{w} \subset\Omega_{ml}$. Now, we must prove that $\mathbf{v}= {\bf 0}$ \for all $x \in\Omega_s$. We fix a compact set $K \subset \Omega_s$, and since the solid domain $\Omega_s$ is open, there exists $\epsilon_K>0$ such that $f_s(c_\epsilon(x),\theta_\epsilon (x))=1$ for $x \in K$, whenever $\epsilon \in (0,\epsilon_K)$. This implies \begin{equation} \label{e6.9} F_i^\epsilon (c_\epsilon(x), \theta_\epsilon (x))=\frac{C_0}{\epsilon^3} \quad\textrm{for}\ x \in K. \end{equation} Choosing $\mathbf{w} = \mathbf{v}_\epsilon $ in \eqref{e4.15}, we have \[ \| \nu^{1/2}(c_\epsilon, \theta_\epsilon ) e (\mathbf{v}_\epsilon) \|^2 +(F_i^\epsilon(c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon, \mathbf{v}_\epsilon) = (\mathbf{F}_e(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon) \] and consequently from \eqref{e4.26} and Proposition \ref{prop4.4}, we have \[ (F_i^\epsilon(c_\epsilon, \theta_\epsilon) \mathbf{v}_\epsilon, \mathbf{v}_\epsilon) \leq (\mathbf{F}_e(c_\epsilon, \theta_\epsilon), \mathbf{v}_\epsilon) \leq C ( 1+\|\theta_\epsilon\|_{0,3} ) | e (\mathbf{v}_\epsilon)|_{0,\Omega} \leq C, \] where $C$ does not depends on $\epsilon$. Therefore, the last inequality and \eqref{e6.9}, implies \[ \frac{C_0}{\epsilon^3} | e (\mathbf{v}_\epsilon)|^2_{0,K} \leq C. \] So, as $\epsilon \to 0$ the term ${\frac{C_0}{\epsilon^3}}\to \infty$, this compels the term $| e (\mathbf{v}_\epsilon)|^2_{0,K}$ to converge to $0$ and consequently $| e (\mathbf{v}_\epsilon)|_{0,K}$ converge to $0$. Then $\mathbf{v} =\mathbf{0}$ on $K$ and the arbitrary choice of compact set $K\subset \Omega_s$ mean that $\mathbf{v}=\mathbf{0}$ in $\Omega_s$. \end{proof} \begin{thebibliography}{99} \bibitem{ad} R. Adams, \textit{Sobolev Spaces}, Academic Press New York (1975). \bibitem{ah} N. M. Ahmad, \textit{Numerical Simulation of Transport Proceses in Multicomponent System realted to Solidification Problems}, Ph. D. Thesis DMS-EPFL 1349 (1995). \bibitem{ba-du} A. Badillo \& M. Dur\'an, \textit{New Trends in Modelling of Solidification Processes for Cooper Alloys}, Copper 2003- Cobre 2003, Vol. I-Plenary Lectures, Economics and Applications of Cooper (2003). \bibitem{ba-du-or} A. Badillo, M. Dur\'an \& E. Ortega-Torres, \textit{C\'alculo de Inestabilidades de un Proceso de Solidificaci\'on en Dominios a Simetr\'{\i}a Cil\'{\i}ndrica, Parte II: Resultados Num\'ericos}, Revista Internacional de M\'etodos Num\'ericos para C\'alculo y Dise\~no en Ingenier\'{\i}a \textbf{21}, 4 (2005), 307-325. \bibitem{be-sa} A. Berm\'udez \& C. Saguez, \textit{Mathematical formulation and numerical solution of an solidification problem}, Free boundary problems: theory and applications, A. Fasano \& M. Primicerio (Eds.), \textbf{1} Pitmann, (1983), 237-247. \bibitem{bl-gas} Ph. Blanc \& L. Gasser, \textit{Existence of a weak solution of a binary alloy problem}, Technical Report 12.92 DMA-EPFL (1992). \bibitem{bl-gas-ra} Ph. Blanc, L. Gasser and J. Rappaz, \textit{Existence for a stationary model of binary alloy solidification}, Math. Model. and Num. Anal. \textbf{29} (1995), 687-699. \bibitem{co-le} H. Combeau \& G. Lesoult, \textit{Simulation of freckles formation and related segregation during directional solidification of metallic alloys}, Modelling of casting, welding and advanced solidification processes 6 (Voller V., Piwonka T.S. \& Katgerman L. eds.) (1993), 201-208. \bibitem{don} D. P. Donnelly, \textit{A model for non-equilibrium thermodynamic processes involving phase changes}, J. Inst. Math. Appl. \textbf{24} (1979), 425-438. \bibitem{du-fe-ro} M. Dur\'an, J. Ferreira \& M. Rojas, \textit{Reproductive weak solutions of magneto-micropolar fluid equations in exterior domains}, Math. Comput. Modelling \textbf{35}, 7-8 (2002), 779-791. \bibitem{du-or-ra} M. Dur\'an, E. Ortega-Torres \& J. Rappaz, \textit{Weak Solutions of a Stationary Convection-Diffusion Model Describing Binary Allow Solidification Processes}, Math. Comput. Modelling \textbf{42}, 11-12 (2005), 1269-1286. \bibitem{du-or-st} M. Dur\'an, E. Ortega-Torres \& R. Stein, \textit{C\'alculo de Inestabilidades de un Proceso de Solidificaci\'on en Dominios a Simetr\'{\i}a Cil\'{\i}ndrica, Parte I: Formulaci\'on del Modelo}, Revista Internacional de M\'etodos Num\'ericos para C\'alculo y Dise\~no en Ingenier\'{\i}a \textbf{17}, 2 (2001), 127-148. \bibitem{eva} L. C. Evans, \textit{Partial Differential Equations}. Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, (1998). \bibitem{gai} F. Gaillard, \textit{Mod\'elisation et analyse num\'erique d'instabilit\'es d\'evelopp\'ess lors de la solidification d'alliages binaires}, Ph. D. Thesis DMS-EPFL (1998). \bibitem{gai-ra} F. Gaillard \& J. Rappaz, \textit{Analysis and numerical simulation for models of binary alloy solification}. Periaux Contributed Book, Jonh Wiley \& Sons Ltd. (1997). \bibitem{gas} L. Gasser, \textit{Existence analysis and numerical schemes for models of binary alloy solidification}, Ph. D. Thesis DMS-EPFL 1421 (1995). \bibitem{gi-tr} D. Gilbarg \& N. Trudinger, \textit{Elliptic Partial Differential Equations of Second Order}. Springer, (1977). \bibitem{hi-lo-ro} N. Hills, D. E. Loper \& P.H. Roberts, \textit{A thermodynamically consistent model of a mushy zone}, Q. J. Mech. Appl. Math. \textbf{36} (1983), 505-539. \bibitem{Lions} J.-L. Lions, \textit{Quelques M\'ethodes de R\'esolution des Probl\`emes aux Limites Non-Lin\'eaires}, Dunod Ed. Paris (1969). \bibitem{luc-vi} S. Luckhaus \& A. Visintin, \textit{Phase transition in multicomponent systems}, Manuscripta math. \textbf{43} (1983), 261-288. \bibitem{mi-th} L. M. Milne-Thomson, \textit{Theoretical Hydrodinamics}, Dover publications, Inc. New York (1996). \bibitem{pr-vo} C. Prakrash \& V. R. Voller, \textit{A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems}, Int. J. Heat Mass Transfer. \textbf{30} (1987), 1709-1719. \bibitem{ra-vo} M. Rappaz \& V.R. Voller, \textit{Modelling of micro-macrosegragation in solidification processes}, Metall. Trans. \textbf{A 21 A} (1990), 749--753. \bibitem{re-si} M. Reed \& B. Simon, \textit{Methods of Modern Mathematical Physics, I: Functional Analysis}. Academic Press Inc., Orlando, (1980). \bibitem{sh} A. H. Shapiro, \textit{The Dynamics and Thermodynamics of Compressible Flow}, Ronald Press, New York (1953). \end{thebibliography} \end{document}