Two-Dimensional Steady Boussinesq Convection: Existence, Computation and Scaling

Abstract

This research investigates laser-induced convection through a stream function-vorticity formulation. Specifically, this paper considers a solution to the steady Boussinesq Navier–Stokes equations in two dimensions with a slip boundary condition on a finite box. A fixed-point algorithm is introduced in stream function-vorticity variables, followed by a proof of the existence of steady solutions for small laser amplitudes. From this analysis, an asymptotic relationship is demonstrated between the nondimensional fluid parameters and least upper bounds for laser amplitudes that guarantee existence, which accords with numerical results implementing the algorithm in a finite difference scheme. The findings indicate that the upper bound for laser amplitude scales by $\mathcal{O}\left({\mathrm{Re}}^{-2}{\mathrm{Pe}}^{-1}{\mathrm{Ri}}^{-1}\right)$ when $\mathrm{Re}\gg \mathrm{Pe}$, and by $\mathcal{O}\left({\mathrm{Re}}^{-1}{\mathrm{Pe}}^{-2}{\mathrm{Ri}}^{-1}\right)$ when $\mathrm{Pe}\gg \mathrm{Re}$. These results suggest that the existence of steady solutions is heavily dependent on the size of the Reynolds (Re) and Peclet (Pe) numbers, as noted in previous studies. The simulations of steady solutions indicate the presence of symmetric vortex rings, which agrees with experimental results described in the literature. From these results, relevant implications to thermal blooming in laser propagation simulations are discussed.

1. Introduction

In this article, we investigate steady Boussinesq fluid convection driven by a laser source. More precisely, we present and interrogate a solution to the steady-state case of a buoyancy-driven fluid with a forcing term due to heat transfer from a laser. The Boussinesq approximation is highly accurate for such buoyancy-driven flows and is commonplace in laser propagation studies, especially those involving thermal blooming [1]. This phenomenon of thermal blooming plays an important role in dictating the behavior of laser-fluid interaction within the context of laser propagation, and especially in domains such as manufacturing, remote sensing, directed energy, imaging, and several others [2,3,4]. Motivated by the applications in these respective fields, we will consider steady-state laser convection within a 2D setting. Several authors have elucidated the means by which laser-induced convection, both in steady and unsteady regimes, can be applied to engineering and applied science problems. For example, Vela [5] describes a method to manipulate micro-beads at mesoscale velocities within a convection chamber. In a similar analysis, Flores–Flores [6] considers microparticle transport driven by laser-induced convective currents in gas at microscale velocities. In nanofluids, Mourad et al. [7] and Jamshed et al. [8] consider the thermal characteristics of fluid flow undergoing natural convection. At larger fluid velocities, Masoomi [9] considered steady-state convective heat transfer from lasers in the context of additive manufacturing, and Rennie [10] has discussed the impact laser-induced convection can have on material heating times in directed energy applications.

In treating the steady-state laser convection problem, we introduce an iterative numerical method to simulate steady solutions over a 2D finite grid. We provide a proof of existence of steady solutions for small enough laser amplitudes and we give a theorem to describe the scaling between nondimensional fluid parameters and a least upper bound for laser amplitude which permits computable solutions. There exist several papers in the relevant literature that prove convergence of iterations in weak and strong formulations for the Navier–Stokes equations. Work has already been done [11,12] to prove the existence of weak solutions in the 2D and 3D steady Navier–Stokes equations with a Boussinesq approximation and slip boundary conditions, but all in the primitive variable setting and without an implementable iteration to find the solutions. Other studies each considered an unforced problem or a problem with forcing along the boundary of the considered domain [13,14]. By choosing an appropriate function space, convergent iterations exist in at least two dimensions for steady and unsteady cases [15,16,17,18]. The iteration presented in this work, however, is unique in that it provides solutions to an internally forced convection problem. In our formulation, we establish a fixed-point iteration with existence in Sobolev spaces in two dimensions. From this iteration, we implement a finite-difference numerical algorithm to compute solutions over a grid. Previous work surrounding the simulation of laser-induced convection has included discussions about free versus forced convection [19], difficulties with nondimensional scaling [20], difficulties with deriving proper boundary conditions [21], and several other associated problems. Similar finite difference algorithms do exist [22,23] along with methods that implement finite elements [24] and radial basis functions [25]; however, these algorithms do not consider a convection problem with internal forcing. Furthermore, in our analysis of the induced finite-difference scheme, we give particular consideration to the asymptotic scaling of the nondimensional fluid parameters.

Since the Navier–Stokes equations may be represented in either primitive variables (temperature, velocity, pressure), or in stream function-vorticity form, the chosen formulation informs on the techniques used to analyze potential solution methods. Both representations have their benefits. The stream function-vorticity formulation, however, is advantageous in that it implicitly guarantees divergence-free flow without the requirement of a pressure-correction term. In numerical implementations of the primitive Navier–Stokes equations, this pressure correction term often leads to solution inaccuracies arising from inappropriate choices in pressure boundary conditions [26]. With the stream function-vorticity formulation, the choice of boundary conditions for each fluid variable is thus simplified. We consider here a no-friction slip condition in stream function-vorticity variables with homogeneous Dirichlet conditions for temperature. The slip condition is easily manifested in stream-function vorticity form as a homogeneous Dirichlet condition on the fluid stream function and vorticity. This condition allows for free tangential movement along the boundary, allowing for the development of vortex rings within the domain.

The novel findings from this study will include the existence of solutions arising from a numerically implementable iteration, as well as the asymptotic relationship between nondimensional fluid parameters and a maximum laser amplitude. In the following chapters, we discuss the relationship between the physical problem and the stream function-vorticity representation, we introduce and prove convergence of an iteration to solve the Boussinesq Navier–Stokes equations, and we conclude with results and a discussion on the asymptotic scaling of nondimensional fluid parameters with respect to convergence. Following this analysis, we provide numerical solutions to the steady-state and we examine the applications to laser beam propagation.

2. Formulation

We now introduce our approach to mathematically describing the laser-fluid interaction of interest. We consider a model for laser-fluid interaction on a finite, 2D domain such that this domain represents a rectangular cross section along the propagation path (z coordinate) of the laser beam. That is, the domain will be a rectangular subsection of the $xy$ plane where we assume the Gaussian laser intensity profile is centered on the origin of the domain. This assumption would be most valid when considering thermal blooming at the laser aperture. We denote the domain as $\Omega$ and the boundary of the domain as $\partial \Omega$, where we begin by considering a square domain $\Omega =[-2,2]\times [-2,2]$ as depicted in Figure 1.

Several assumptions on the fluid and laser properties must be imposed in the pursuit of governing equations which dictate fluid behavior. We assume an incompressible, Newtonian fluid which is governed by the Boussinesq approximation for buoyancy-driven flows. This Boussinesq model assumes that variations in fluid density have a linear relationship with temperature. Further, we assume a representation of the normalized laser irradiance as a function defined over the domain $\Omega$ such that the laser heats the fluid and induces temperature variations. We first introduce the governing equations for fluid flow in primitive variables where we seek the functions

$$T:\Omega \to \mathbb{R},p:\Omega \to \mathbb{R}\mathbf{u}:\Omega \to {\mathbb{R}}^2$$

representing the nondimensional temperature fluctuation, pressure, and velocity field of the fluid which obey the incompressible, Boussinesq Navier–Stokes equations undergoing forcing [27]:

$$\begin{array}{ccc}\hfill {\mathbf{u}}_t+(\mathbf{u}·\nabla )\mathbf{u}& =& \nabla p+\frac{1}{\mathrm{Re}}\Delta \mathbf{u}+\mathrm{Ri}T{\mathbf{e}}_{\mathbf{2}},\hfill \end{array}$$

(1a)

$$\begin{array}{ccc}\hfill \nabla ·\mathbf{u}& =& 0,\hfill \end{array}$$

(1b)

$$\begin{array}{ccc}\hfill T_t+(\mathbf{u}·\nabla )T& =& \frac{1}{\mathrm{Pe}}\Delta T+F,\hfill \end{array}$$

(1c)

where $F:\Omega \to \mathbb{R}$ represents the laser forcing and ${\mathbf{e}}_{\mathbf{2}}$ represents the unit vector in the vertical y direction. For a characteristic length scale L, characteristic velocity scale U, force due to gravity g, kinematic viscosity $\nu$, and thermal diffusivity $\mu$, the Reynolds, Peclet, and Richardson numbers are respectively defined as

$$\mathrm{Re}=\frac{UL}{\nu },\mathrm{Pe}=\frac{UL}{\mu },\mathrm{Ri}=\frac{gL}{U^2}.$$

(2)

We seek steady-state solutions such that all temporal variation is negligible and each ${\mathbf{u}}_t=T_t=p_t=0$. In two dimensions, it is often easier to consider the stream function-vorticity form of the Boussinesq equations. This formulation reduces the total number of fluid variables by one, and often allows for easier treatment of boundary conditions. The transformation from the primitive form to the stream function-vorticity form is a well established procedure, where the two representations are equivalent for valid boundary data [28]. We state the equivalent stream function-vorticity form of the steady Boussinesq equations below:

$$\begin{array}{}\mathrm{(3a)}& \hfill \Delta T& =& \mathrm{Pe}\left({\psi }_y{T}_x-{\psi }_x{T}_y-F\right),\hfill \mathrm{(3b)}& \hfill \Delta \omega & =& \mathrm{Re}\left({\psi }_y{\omega }_x-{\psi }_x{\omega }_y\right)+\mathrm{Ri}\mathrm{Re}{T}_x,\hfill \mathrm{(3c)}& \hfill \Delta \psi & =& \omega ,\hfill \end{array}$$

where $\omega$ is the vorticity and $\psi$ is the stream function of the fluid. Like the primitive formulation, we may interpret the system as a type of vector-valued, nonlinear Poisson equation over the domain $\Omega$.

In this study, we consider a slip boundary condition on a finite box corresponding to zero shear stress along the boundary. This slip condition enforces boundary impermeability and assumes that the coefficient of friction between the fluid and the boundary is zero. That is, the normal component of velocity at the boundary vanishes but the tangential component is unrestricted. This boundary condition is useful in investigating laser-fluid interaction as it localizes the fluid within a finite domain without having to consider boundary layer effects. In primitive variables, where $\mathbf{u}=(u,v)$, this is satisfied through the condition

$$\left\{\begin{array}{cc}v=0,\hfill & (x,y)\in \partial {\Omega }_{\mathrm{top}\&\mathrm{bottom}},\hfill \\ \frac{\partial u}{\partial y}=0,\hfill & (x,y)\in \partial {\Omega }_{\mathrm{top}\&\mathrm{bottom}},\hfill \\ u=0,\hfill & (x,y)\in \partial {\Omega }_{\mathrm{left}\&\mathrm{right}},\hfill \\ \frac{\partial v}{\partial x}=0,\hfill & (x,y)\in \partial {\Omega }_{\mathrm{left}\&\mathrm{right}}.\hfill \end{array}\right.$$

(4)

Based on the definitions for the stream function

$$\frac{\partial \psi }{\partial y}=u,\frac{\partial \psi }{\partial x}=-v,$$

(5)

and the definition for vorticity

$$\omega =\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y},$$

(6)

we may transform the primitive slip condition into an equivalent condition on the stream function and vorticity. First, if $v=0$ on the top and bottom and $u=0$ on the left and right, then

$${\psi }_x=0,(x,y)\in \partial {\Omega }_{\mathrm{top}\&\mathrm{bottom}}$$

and

$${\psi }_y=0,(x,y)\in \partial {\Omega }_{\mathrm{left}\&\mathrm{right}}$$

which is satisfied by taking $\psi =0$ on all of $\partial \Omega$. For the vorticity, if $u_y=0$ on the top and bottom, then $\omega =v_x=0$ on the top and bottom. Similarly, if $v_x=0$ on the left and right, then $\omega =u_y=0$ on the left and right. It thus follows that the slip condition in primitive variables is equivalent to a homogeneous Dirichlet condition on the stream function and vorticity. For simplicity, we impose a homogeneous Dirichlet condition on the temperature, but a Neumann condition can also yield useful results. The slip condition in stream function-vorticity form is thus stated as

$$\left\{\begin{array}{cc}\omega (x,y)=0,\hfill & (x,y)\in \partial \Omega ,\hfill \\ \psi (x,y)=0,\hfill & (x,y)\in \partial \Omega ,\hfill \\ T(x,y)=0,\hfill & (x,y)\in \partial \Omega .\hfill \end{array}\right.$$

(7)

With the stream function-vorticity formulation posed as a system of coupled Poisson boundary value problems, we now define $\widehat{T}$, $\widehat{\omega }$, and $\widehat{\psi }$ as the exact temperature fluctuation, vorticity, and stream function, respectively. We restate Equation (3) for the steady Navier–Stokes equations in the Boussinseq approximation as

$$\begin{array}{}\mathrm{(8a)}& \Delta \widehat{T}\hfill & =& \mathrm{Pe}\left({\widehat{\psi }}_y{\widehat{T}}_x-{\widehat{\psi }}_x{\widehat{T}}_y-F\right),\hfill \mathrm{(8b)}& \Delta \widehat{\omega }\hfill & =& \mathrm{Re}\left({\widehat{\psi }}_y{\widehat{\omega }}_x-{\widehat{\psi }}_x{\widehat{\omega }}_y\right)+\mathrm{Ri}\mathrm{Re}{\widehat{T}}_x,\hfill \mathrm{(8c)}& \Delta \widehat{\psi }\hfill & =& \widehat{\omega },\hfill \end{array}$$

We then introduce the fixed-point iteration for $n>1$ as

$$\begin{array}{}\mathrm{(9a)}& \hfill \Delta {\widehat{T}}_n& =& \mathrm{Pe}\left(\frac{\partial {\widehat{\psi }}_{n-1}}{\partial y}\frac{\partial {\widehat{T}}_{n-1}}{\partial x}-\frac{\partial {\widehat{\psi }}_{n-1}}{\partial x}\frac{\partial {\widehat{T}}_{n-1}}{\partial y}-\epsilon f\right),\hfill \mathrm{(9b)}& \hfill \Delta {\widehat{\omega }}_n& =& \mathrm{Re}\left(\frac{\partial {\widehat{\psi }}_{n-1}}{\partial y}\frac{\partial {\widehat{\omega }}_{n-1}}{\partial x}-\frac{\partial {\widehat{\psi }}_{n-1}}{\partial x}\frac{\partial {\widehat{\omega }}_{n-1}}{\partial y}\right)+\mathrm{Ri}\mathrm{Re}\frac{\partial {\widehat{T}}_n}{\partial x},\hfill \mathrm{(9c)}& \hfill \Delta {\widehat{\psi }}_n& =& {\widehat{\omega }}_n.\hfill \end{array}$$

where the laser irradiance takes the form of an amplitude $\epsilon$ scaled by a normalized Gaussian f

$$F(x,y)=\epsilon e^{-c(x^2+y^2)}=\epsilon f(x,y),$$

(10)

and we initialize the solution as

$$\begin{array}{}\mathrm{(11a)}& \hfill \Delta {\widehat{T}}_0& =& -\mathrm{Pe}\left(\epsilon f\right),\hfill \mathrm{(11b)}& \hfill \Delta {\widehat{\omega }}_0& =& \mathrm{Ri}\mathrm{Re}\frac{\partial {\widehat{T}}_0}{\partial x},\hfill \mathrm{(11c)}& \hfill \Delta {\widehat{\psi }}_0& =& {\widehat{\omega }}_0.\hfill \end{array}$$

This initialization is the linearization of the system of equations and thus should serve as a strong initial guess for small amplitude solutions. We show in the following section that unique, convergent solutions exist for this iteration when laser amplitude $\epsilon$ is small.

In the numerical implementation of this iteration, we employ a finite difference scheme to successively solve the Poisson equations in each fluid variable at each step n. We proceed by uniformly discretizing the domain $\Omega$ with step size h, and we utilize a second order, centered difference approximation to discretize the Laplacian

$${(\Delta u)}_{ij}=\frac{u_{i+1,j}+u_{i-1,j}-4u_{i,j}+u_{i,j+1}+u_{i,j-1}}{h^2}+\mathcal{O}\left(h^2\right).$$

(12)

For an $N\times N$ discrete grid, the solution to the discrete Poisson equation solution can be found in $\mathcal{O}(N^{2}log{N}^2)$ flops [29]. Since this represents the most computationally expensive step in the iteration, we thus conclude that the scheme (9) has an asymptotic cost of $\mathcal{O}(N^{2}log{N}^2)$ when it converges to a unique solution. In a MATLAB implementation, we establish stopping criteria based on the Cauchy error of successive iterations. That is, the algorithm terminates when the numerical solutions converge based on the criterion

$$\parallel T_n-T_{n-1}{\parallel }_2+\parallel {\omega }_n-{\omega }_{n-1}{\parallel }_2+{\parallel {\psi }_n-{\psi }_{n-1}\parallel }_2<p$$

(13)

using vector 2-norms for some appropriately small p. The following section is dedicated to demonstrating that these convergent solutions exist for small enough laser amplitudes.

3. Existence Proof

Recall that the iteration (9) was constructed so that its fixed points are steady solutions of (8), so this convergence of the scheme gives existence of solutions. We state this result as a theorem:

Theorem 1

(Convergence of Steady Solutions to the Stream Function-Vorticity Form of the Boussinesq Navier–Stokes equations). The iteration defined in (9) and (11) converges to unique solutions for small forcing and Dirichlet boundary conditions.

Proof. 

We begin by stating a definition from Brezis [30] and two lemmas which assist in our analysis:

Definition 1 

(Sobolev Spaces). Let $\Omega \subset {\mathbb{R}}^N$ and let $p\in \mathbb{R}$ with $1\le p\le \infty$. If $m\ge 2$ is an integer, then the Sobolev space $W^{m,p}(\Omega )$ is defined for all α with $\left|\alpha \right|\le m$,

$$W^{m,p}(\Omega )=\{\left.u\in L^p(\Omega )\right|\exists g_{\alpha }\in L^p(\Omega ):{\int }_{\Omega }u{D}^{\alpha }\phi ={(-1)}^{\left|\alpha \right|}{\int }_{\Omega }{g}_{\alpha }\phi \forall \phi \in C_c^{\infty }(\Omega )\},$$

(14)

where we employ the multi-index notation $\alpha =({\alpha }_1,{\alpha }_2,...,{\alpha }_N)$ such that

$$\left|\alpha \right|=\sum _{i=1}^N{\alpha }_i,D^{\alpha }\phi =\frac{{\partial }^{\left|\alpha \right|}\phi }{\partial x_1^{{\alpha }_1}\partial x_2^{{\alpha }_2}...\partial x_N^{{\alpha }_N}}.$$

(15)

We assign the $g_{\alpha }=D^{\alpha }u$. The space $W^{m,p}(\Omega )$ equipped with the norm

$${\parallel u\parallel }_{W^{m,p}}=\sum _{\left|\alpha \right|\le m}{\parallel D^{\alpha }u\parallel }_{L^p}$$

(16)

is a Banach space and we define the space

$$H^m(\Omega )=W^{m,2}(\Omega )$$

(17)

originally defined by Sobolev in [31,32].

Lemma 1 

($H^s$ is an algebra). Let $\Omega \subset {\mathbb{R}}^d$. Suppose $s_1,s_2\ge s$ and $s_1+s_2>s+d/2$. Then $u\in H^{s_1}(\Omega )$ and $v\in H^{s_2}(\Omega )$ gives that $\left(uv\right)\in H^s(\Omega )$, and that

$${\parallel uv\parallel }_{H^s}\le {M\parallel u\parallel }_{H^{s_1}}{\parallel v\parallel }_{H^{s_2}}$$

where M depends only on $s_1,s_2$ and d.

See [33] for proof.

Lemma 2 

(Elliptic Estimate). If $F\in L^p(\Omega )$ then there exists a unique solution to

$$\Delta u=Fx\in \Omega =[a,b]\times [c,d]$$

with

$${u|}_{\partial \Omega }=0,$$

And that solution satisfies ${\parallel u\parallel }_{W^{2,p}(\Omega )}\le C{\parallel F\parallel }_{L^p(\Omega )}$ where $C=C(a,b,c,d,p)$ is a positive constant.

(see [34] for proof. A similar result also exists for inhomogeneous and Neumann boundary conditions for boundary data which satisfy compatibility conditions as described in [35,36].)

We now introduce the small amplitude assumption on the dynamics, where each $\widehat{T},\widehat{\omega },\widehat{\psi }$ is represented by a perturbation from the resting fluid case:

$$\begin{array}{}\mathrm{(18a)}& \hfill \widehat{T}& =\epsilon T,\mathrm{(18b)}& \hfill \widehat{\omega }& =\epsilon \omega ,\mathrm{(18c)}& \hfill \widehat{\psi }& =\epsilon \psi .\end{array}$$

The key result of the proof involves constructing the equation for the difference between two iterations $T_n$ and $T_m$ (and similarly ${\omega }_n$ and ${\omega }_m$, ${\psi }_n$ and ${\psi }_m$) by applying the representation in (18) and subtracting (9) at the two indices m and n:

$$\begin{array}{}\mathrm{(19a)}& \hfill \begin{array}{cc}\hfill & \Delta (T_n-T_m)\hfill \\ \hfill & =\epsilon \mathrm{Pe}\left[\left(\frac{\partial {\psi }_{n-1}}{\partial y}\frac{\partial T_{n-1}}{\partial x}-\frac{\partial {\psi }_{n-1}}{\partial x}\frac{\partial T_{n-1}}{\partial y}\right)-\left(\frac{\partial {\psi }_{m-1}}{\partial y}\frac{\partial T_{m-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\frac{\partial T_{m-1}}{\partial y}\right)\right],\hfill \end{array}\mathrm{(19b)}& \hfill \begin{array}{cc}\hfill & \Delta ({\omega }_n-{\omega }_m)\hfill \\ & =\epsilon \mathrm{Re}\left[\left(\frac{\partial {\psi }_{n-1}}{\partial y}\frac{\partial {\omega }_{n-1}}{\partial x}-\frac{\partial {\psi }_{n-1}}{\partial x}\frac{\partial {\omega }_{n-1}}{\partial y}\right)-\left(\frac{\partial {\psi }_{m-1}}{\partial y}\frac{\partial {\omega }_{m-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\frac{\partial {\omega }_{m-1}}{\partial y}\right)\right]\hfill \\ \hfill & +\mathrm{Ri}\mathrm{Re}\left(\frac{\partial T_n}{\partial x}-\frac{\partial T_m}{\partial x}\right),\hfill \end{array}\mathrm{(19c)}& \begin{array}{cc}\hfill & \Delta ({\psi }_n-{\psi }_m)={\omega }_n-{\omega }_m.\hfill \end{array}\hfill \end{array}$$

We proceed by investigating the form of Equation (19a) first. Define

$$F_T(n,m)=\epsilon \mathrm{Pe}\left[\left(\frac{\partial {\psi }_{n-1}}{\partial y}\frac{\partial T_{n-1}}{\partial x}-\frac{\partial {\psi }_{n-1}}{\partial x}\frac{\partial T_{n-1}}{\partial y}\right)-\left(\frac{\partial {\psi }_{m-1}}{\partial y}\frac{\partial T_{m-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\frac{\partial T_{m-1}}{\partial y}\right)\right]$$

(20)

where we manipulate the form of the equation to arrive at a more convenient expression:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & F_T(n,m)=\epsilon \mathrm{Pe}\left[\left(\frac{\partial {\psi }_{n-1}}{\partial y}\frac{\partial T_{n-1}}{\partial x}-\frac{\partial {\psi }_{n-1}}{\partial x}\frac{\partial T_{n-1}}{\partial y}\right)-\left(\frac{\partial {\psi }_{m-1}}{\partial y}\frac{\partial T_{m-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\frac{\partial T_{m-1}}{\partial y}\right)\right]\hfill \\ \hfill & =\frac{\epsilon \mathrm{Pe}}{2}\left[\left(\frac{\partial {\psi }_{n-1}}{\partial y}-\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}+\frac{\partial T_{m-1}}{\partial x}\right)+\left(\frac{\partial {\psi }_{n-1}}{\partial y}+\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}-\frac{\partial T_{m-1}}{\partial x}\right)\right.\hfill \\ \hfill & -\left.\left(\frac{\partial {\psi }_{n-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}+\frac{\partial T_{m-1}}{\partial y}\right)-\left(\frac{\partial {\psi }_{n-1}}{\partial x}+\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}-\frac{\partial T_{m-1}}{\partial y}\right)\right].\hfill \end{array}\end{array}$$

(21)

We now prove that the function $F_T(n,m)\in L^2(\Omega )$ such that the norm $\parallel F_T{\parallel }_{L^2(\Omega )}<\infty$ by showing that each $T_n,{\omega }_n,{\psi }_n\in H^2(\Omega )$ for $n,m\ge 0$.

Lemma 3 

(Sobolev Inclusion). Let $f(x,y)$ be a laser forcing such that $f\in L^2(\Omega )$. Then for all $n,m\ge 0$, each $T_n,{\omega }_n,{\psi }_n\in H^2(\Omega )$ and $F_T(n,m)\in L^2(\Omega )$ as defined in (20).

Proof. 

Suppose that $f(x,y)\in L^2(\Omega )$. By recalling the initialization (11), we apply Lemma 2 to see that $T_0\in H^2(\Omega )$. By the definition of the Sobolev space (14), it follows that $D^1{T}_0\in L^2(\Omega )$ and thus by applying Lemma 2 again it follows that ${\omega }_0\in H^2(\Omega )$ and similarly ${\psi }_0\in H^2(\Omega )$. We now hypothesize for induction that each $T_n,{\omega }_n,{\psi }_n\in H^2(\Omega )$. The base case was just shown. Based on the iteration (9) and by applying the small amplitude assumption (18), we see that

$$\Delta T_{n+1}=\epsilon \mathrm{Pe}\left[\left(\frac{\partial {\psi }_n}{\partial y}\frac{\partial T_n}{\partial x}-\frac{\partial {\psi }_n}{\partial x}\frac{\partial T_n}{\partial y}\right)-f\right].$$

Since ${\psi }_n$ and $T_n\in H^2(\Omega )$ by assumption, we apply Lemma 1 with $s=0$ and $s_1=1$ on the $T_n$ derivatives and $s_2=1$ on the ${\psi }_n$ derivatives to see that the right hand side is an $L^2(\Omega )$ function such that $T_{n+1}\in H^2(\Omega )$ from Lemma 2. A similar argument holds for each ${\omega }_{n+1},{\psi }_{n+1}$ so we conclude that $T_n,{\omega }_n,{\psi }_n\in H^2(\Omega )$. Based on this Sobolev inclusion, we may again apply Lemma 1 to the terms in the expression for $F_T(n,m)$ to see that $F_T(n,m)\in L^2(\Omega )$.   □

We now apply the triangle inequality to the norm of the expression (21) which yields

$$\begin{array}{cc}\hfill \parallel F_T{\parallel }_{L^2(\Omega )}& \le \frac{\epsilon \mathrm{Pe}}{2}\left[\parallel \left(\frac{\partial {\psi }_{n-1}}{\partial y}-\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}+\frac{\partial T_{m-1}}{\partial x}\right){\parallel }_{L^2}\right.\hfill \\ & +\parallel \left(\frac{\partial {\psi }_{n-1}}{\partial y}+\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}-\frac{\partial T_{m-1}}{\partial x}\right){\parallel }_{L^2}\hfill \\ \hfill & +\parallel \left(\frac{\partial {\psi }_{n-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}+\frac{\partial T_{m-1}}{\partial y}\right){\parallel }_{L^2}\hfill \\ \hfill & +\left.\parallel \left(\frac{\partial {\psi }_{n-1}}{\partial x}+\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}-\frac{\partial T_{m-1}}{\partial y})\right){\parallel }_{L^2}\right].\hfill \end{array}$$

(22)

Recall that we can define the Sobolev space

$$H^2(\Omega )=\{f\in L^2(\Omega ):D^{\alpha }f\in L^2(\Omega )\forall \left|\alpha \right|\le 2\}.$$

The norm for $f\in H^2(\Omega )$ is then

$${\parallel f\parallel }_{H^2(\Omega )}=\sum _{\left|\alpha \right|\le 2}{\parallel D^{\alpha }f\parallel }_{L^2(\Omega )}.$$

By this definition, it then follows that

$$\frac{\partial {\psi }_{n-1}}{\partial x},\frac{\partial {\psi }_{m-1}}{\partial x},\frac{\partial {\psi }_{n-1}}{\partial y},\frac{\partial {\psi }_{m-1}}{\partial y},\frac{\partial T_{n-1}}{\partial x},\frac{\partial T_{m-1}}{\partial x},\frac{\partial T_{n-1}}{\partial y},\frac{\partial T_{m-1}}{\partial y}\in L^2(\Omega ).$$

Hence, by applying Lemma 1 with $s=0$ and $s_1=s_2=1$, we can establish the following inequalities on the norms:

$$\begin{array}{cc}\hfill & \parallel \left(\frac{\partial {\psi }_{n-1}}{\partial y}-\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}+\frac{\partial T_{m-1}}{\partial x}\right){\parallel }_{L^2(\Omega )}\hfill \\ \hfill \le & M\parallel \frac{\partial {\psi }_{n-1}}{\partial y}-\frac{\partial {\psi }_{m-1}}{\partial y}{\parallel }_{H^2(\Omega )}{\parallel \frac{\partial T_{n-1}}{\partial x}+\frac{\partial T_{m-1}}{\partial x}\parallel }_{L^2(\Omega )},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \parallel \left(\frac{\partial {\psi }_{n-1}}{\partial y}+\frac{\partial {\psi }_{m-1}}{\partial y}\right)\left(\frac{\partial T_{n-1}}{\partial x}-\frac{\partial T_{m-1}}{\partial x}\right){\parallel }_{L^2(\Omega )}\hfill \\ \hfill \le & M\parallel \frac{\partial {\psi }_{n-1}}{\partial y}+\frac{\partial {\psi }_{m-1}}{\partial y}{\parallel }_{H^2(\Omega )}{\parallel \frac{\partial T_{n-1}}{\partial x}-\frac{\partial T_{m-1}}{\partial x}\parallel }_{L^2(\Omega )},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \parallel \left(\frac{\partial {\psi }_{n-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}+\frac{\partial T_{m-1}}{\partial y}\right){\parallel }_{L^2(\Omega )}\hfill \\ \hfill \le & M\parallel \frac{\partial {\psi }_{n-1}}{\partial x}-\frac{\partial {\psi }_{m-1}}{\partial x}{\parallel }_{H^2(\Omega )}{\parallel \frac{\partial T_{n-1}}{\partial y}+\frac{\partial T_{m-1}}{\partial y}\parallel }_{L^2(\Omega )},\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \parallel \left(\frac{\partial {\psi }_{n-1}}{\partial x}+\frac{\partial {\psi }_{m-1}}{\partial x}\right)\left(\frac{\partial T_{n-1}}{\partial y}-\frac{\partial T_{m-1}}{\partial y}\right){\parallel }_{L^2(\Omega )}\hfill \\ \hfill \le & M\parallel \frac{\partial {\psi }_{n-1}}{\partial x}+\frac{\partial {\psi }_{m-1}}{\partial x}{\parallel }_{H^2(\Omega )}{\parallel \frac{\partial T_{n-1}}{\partial y}-\frac{\partial T_{m-1}}{\partial y}\parallel }_{L^2(\Omega )}.\hfill \end{array}$$

Note that by definition of the norm (16) in the Sobolev space $H^2(\Omega )$ for some $\phi \in H^2(\Omega )$, we can establish the inequality

$$\parallel \frac{\partial \phi }{\partial x}{\parallel }_{L^2(\Omega )},\parallel \frac{\partial \phi }{\partial y}{\parallel }_{L^2(\Omega )}\le {\parallel \phi \parallel }_{H^2(\Omega )}.$$

Finally, by applying Lemma 2 to the Equation (19a), we see that for $\parallel F_T{(n,m)\parallel }_{L^2(\Omega )}<\infty$, $T_n-T_m\in H_0^2(\Omega )$ with $\parallel T_n-T_m{\parallel }_{H^2(\Omega )}\le C{\parallel F_T\parallel }_{L^2(\Omega )}$. We now establish the inequality on $\parallel T_n-T_m{\parallel }_{H^2(\Omega )}$ by:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel T_n-T_m{\parallel }_{H^2(\Omega )}& \le \epsilon \mathrm{Pe}CM\left(\parallel ({\psi }_{n-1}-{\psi }_{m-1}){\parallel }_{H^2(\Omega )}{\parallel (T_{n-1}+T_{m-1})\parallel }_{H^2(\Omega )}\right.\hfill \\ \hfill & +\left.\parallel ({\psi }_{n-1}+{\psi }_{m-1}){\parallel }_{H^2(\Omega )}{\parallel (T_{n-1}-T_{m-1})\parallel }_{H^2(\Omega )}\right).\hfill \end{array}\end{array}$$

(23)

With a similar approach, we provide the equivalent inequalities for $\omega$ and $\psi$:

$$\begin{array}{cc}\hfill \parallel {\omega }_n-{\omega }_m{\parallel }_{H^2(\Omega )}& \le \epsilon \mathrm{Re}CM\left(\parallel ({\psi }_{n-1}-{\psi }_{m-1}){\parallel }_{H^2(\Omega )}{\parallel ({\omega }_{n-1}+{\omega }_{m-1})\parallel }_{H^2(\Omega )}\right.\hfill \\ \hfill & +\left.\parallel ({\psi }_{n-1}+{\psi }_{m-1}){\parallel }_{H^2(\Omega )}{\parallel ({\omega }_{n-1}-{\omega }_{m-1})\parallel }_{H^2(\Omega )}\right)\hfill \\ \hfill & +\mathrm{Ri}\mathrm{Re}C\parallel T_n-T_m{\parallel }_{H^2(\Omega )}\hfill \end{array}$$

(24)

$$\parallel {\psi }_n-{\psi }_m{\parallel }_{H^2(\Omega )}\le C{\parallel {\omega }_n-{\omega }_m\parallel }_{H^2(\Omega )}$$

(25)

With a bound established on the difference between two terms in the iteration, we now show that for small enough $\epsilon$ the entire sequence remains in a small ball.

Lemma 4 

(Bounded Sequence). There exists ${\epsilon }_1>0$ and $R\left({\epsilon }_1\right)<\infty$ for which the sequence defined by the iteration (9) with initialization (11) is bounded in a ball for all n.

Proof. 

We show by strong induction that the sequence (9) is bounded by the ball ${\mathcal{B}}_0\left(R\right)$ centered at zero with radius R. Hence, let $R\left({\epsilon }_1\right)$ be defined such that

$$\parallel T_0{\parallel }_{H^2},\parallel T_1{\parallel }_{H^2},\parallel {\omega }_0{\parallel }_{H^2},\parallel {\omega }_1{\parallel }_{H^2},\parallel {\psi }_0{\parallel }_{H^2},{\parallel {\psi }_1\parallel }_{H^2}\le \frac{R}{2},$$

where we say

$$R=2max\{\parallel T_0{\parallel }_{H^2},\parallel T_1{\parallel }_{H^2},\parallel {\omega }_0{\parallel }_{H^2},\parallel {\omega }_1{\parallel }_{H^2},\parallel {\psi }_0{\parallel }_{H^2},\parallel {\psi }_1{\parallel }_{H^2},1\}$$

to keep R sufficiently large where each $H^2$ norm is measured over the domain $\Omega$. This establishes the base case for induction. Now, we hypothesize by induction that each

$$\parallel T_{n-1}{\parallel }_{H^2},\parallel T_{n-2}{\parallel }_{H^2},...,\parallel {\omega }_{n-1}{\parallel }_{H^2},\parallel {\omega }_{n-2}{\parallel }_{H^2},....,\parallel {\psi }_{n-1}{\parallel }_{H^2},{\parallel {\psi }_{n-2}\parallel }_{H^2},....<R$$

where we seek to show that $\parallel T_n{\parallel }_{H^2},\parallel {\omega }_n{\parallel }_{H^2},{\parallel {\psi }_n\parallel }_{H^2}<R$. Thus, by the inequality (23), it follows that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel T_n-T_1{\parallel }_{H^2}& \le 2\epsilon \mathrm{Pe}CM\left(\parallel ({\psi }_{n-1}-{\psi }_0){\parallel }_{H^2}{\parallel (T_{n-1}+T_0)\parallel }_{H^2}\right.\hfill \\ \hfill & +\left.\parallel ({\psi }_{n-1}+{\psi }_0){\parallel }_{H^2}{\parallel (T_{n-1}-T_0)\parallel }_{H^2}\right).\hfill \end{array}\end{array}$$

(26)

By applying the definition of R and the inductive hypothesis it follows that $\parallel {\psi }_{n-1}-{\psi }_0{\parallel }_{H^2}\le \frac{3R}{2}$, $\parallel {\psi }_{n-1}+{\psi }_0{\parallel }_{H^2}\le \frac{3R}{2}$, $\parallel T_{n-1}-T_0{\parallel }_{H^2}\le \frac{3R}{2}$, and $\parallel T_{n-1}+T_0{\parallel }_{H^2}\le \frac{3R}{2}$. Hence, we can apply the inequality

$$\parallel T_n-T_1{\parallel }_{H^2}\le 2\epsilon \mathrm{Pe}CM\left(\left(\frac{3R}{2}\right)\left(\frac{3R}{2}\right)\right)$$

where we select ${\epsilon }_T$ such that for $\epsilon <{\epsilon }_T$, $\parallel T_n-T_1{\parallel }_{H^2}<\frac{R}{2}$. Hence, let

$${\epsilon }_T=\frac{1}{9\mathrm{Pe}CMR}$$

and it follows that if $\parallel T_n-T_1{\parallel }_{H^2}<\frac{R}{2}$ then by the reverse triangle inequality

$$|\parallel T_n{\parallel }_{H^2}-\parallel T_1{\parallel }_{H^2}|\le \parallel T_n-T_1{\parallel }_{H^2}<\frac{R}{2}$$

so

$$\parallel T_n{\parallel }_{H^2}<\frac{R}{2}+{\parallel T_1\parallel }_{H^2}<R.$$

We now proceed by finding the similar terms ${\epsilon }_{\omega }$ and ${\epsilon }_{\psi }$ such that each respective term stays within the ball of radius R. From the inequality (24):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel {\omega }_n-{\omega }_1{\parallel }_{H^2}& \le \epsilon \mathrm{Re}CM\left(\parallel ({\psi }_{n-1}-{\psi }_0){\parallel }_{H^2}{\parallel ({\omega }_{n-1}+{\omega }_0)\parallel }_{H^2}\right.\hfill \\ \hfill & +\left.\parallel ({\psi }_{n-1}+{\psi }_0){\parallel }_{H^2}{\parallel ({\omega }_{n-1}-{\omega }_0)\parallel }_{H^2}\right)\hfill \\ & +\mathrm{Ri}\mathrm{Re}C\parallel T_n-T_1{\parallel }_{H^2}\hfill \end{array}\end{array}$$

(27)

We recall that $\parallel T_n-T_1{\parallel }_{H^2}\le \frac{9\epsilon \mathrm{Pe}CM{R}^2}{2}$ and we then apply a similar analysis to obtain:

$$\parallel {\omega }_n-{\omega }_1{\parallel }_{H^2}\le \frac{9}{2}\mathrm{Re}CM\epsilon R^2+C\mathrm{Ri}\mathrm{Re}\frac{9\epsilon \mathrm{Pe}CM{R}^2}{2}.$$

To force $\parallel {\omega }_n-{\omega }_1{\parallel }_{H^2}<\frac{R}{2}$, we then pick

$${\epsilon }_{\omega }=\frac{1}{(9\mathrm{Re}CM+9C^2\mathrm{Ri}\mathrm{Re}\mathrm{Pe}M)R}.$$

From the inequality (25), we have:

$$\parallel {\psi }_n-{\psi }_1{\parallel }_{H^2}\le C{\parallel {\omega }_n-{\omega }_1\parallel }_{H^2}$$

(28)

So from applying the result for ${\epsilon }_{\omega }$, we obtain

$${\epsilon }_{\psi }=\frac{{\epsilon }_{\omega }}{C}=\frac{1}{(9\mathrm{Re}{C}^{2}M+9C^3\mathrm{Ri}\mathrm{Re}\mathrm{Pe}M)R}.$$

We now let ${\epsilon }_1=inf\{{\epsilon }_T,{\epsilon }_{\omega },{\epsilon }_{\psi }\}$ where for all $\epsilon <{\epsilon }_1$,

$$\parallel T_n{\parallel }_{H^2},\parallel {\omega }_n{\parallel }_{H^2},{\parallel {\psi }_n\parallel }_{H^2}<R$$

so we conclude by induction that there exists ${\epsilon }_1$ such that the sequence remains bounded in a ball of radius R.   □

Theorem 2 

(Contraction Mapping). The iteration defined by (9) and (11) is a contraction mapping.

Proof. 

From the above lemma, we establish by the triangle inequality that each

$$\parallel T_n+T_m{\parallel }_{H^2},\parallel {\omega }_n+{\omega }_m{\parallel }_{H^2},{\parallel {\psi }_n+{\psi }_m\parallel }_{H^2}\le 2R.$$

Hence, by the inequality (23) for $\parallel T_n-T_m{\parallel }_{H^2}$, we have:

$$\parallel T_n-T_m{\parallel }_{H^2}\le 2\epsilon \mathrm{Pe}CMR\left(\parallel ({\psi }_{n-1}-{\psi }_{m-1}){\parallel }_{H^2}+{\parallel (T_{n-1}-T_{m-1})\parallel }_{H^2}\right).$$

Similarly, the other variables satisfy

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel {\omega }_n-{\omega }_m{\parallel }_{H^2}& \le \epsilon R\left(2ReCM+2PeRiRe{C}^{2}M\right){\parallel {\psi }_{n-1}-{\psi }_{m-1}\parallel }_{H^2}\hfill \\ \hfill & +2\epsilon \mathrm{Re}CMR\parallel {\omega }_{n-1}-{\omega }_{m-1}{\parallel }_{H^2}\hfill \\ \hfill & +2\epsilon \mathrm{Pe}\mathrm{Ri}\mathrm{Re}{C}^{2}MR{\parallel T_{n-1}-T_{m-1}\parallel }_{H^2},\hfill \end{array}\end{array}$$

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel {\psi }_n-{\psi }_m{\parallel }_{H^2}& \le \epsilon CR\left(2ReCM+2PeRiRe{C}^{2}M\right){\parallel {\psi }_{n-1}-{\psi }_{m-1}\parallel }_{H^2}\hfill \\ \hfill & +2\epsilon \mathrm{Re}{C}^{2}MR{\parallel {\omega }_{n-1}-{\omega }_{m-1}\parallel }_{H^2}\hfill \\ \hfill & +2\epsilon \mathrm{Pe}\mathrm{Ri}\mathrm{Re}{C}^{3}MR{\parallel T_{n-1}-T_{m-1}\parallel }_{H^2}.\hfill \end{array}\end{array}$$

We now consider the space $S=H^2(\Omega )\times H^2(\Omega )\times H^2(\Omega )$ consisting of triples $X=(X_1,X_2,X_3)$ of $H^2(\Omega )$ functions with metric defined by:

$$d(X,Y)=\parallel X_1-Y_1{\parallel }_{H^2}+\parallel X_2-Y_2{\parallel }_{H^2}+{\parallel X_3-Y_3\parallel }_{H^2}\forall X,Y\in S.$$

Thus, for our set of fluid variables, we establish a new inequality for the sum of differences of index:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \parallel T_n-T_m{\parallel }_{H^2}+\parallel {\omega }_n-{\omega }_m{\parallel }_{H^2}+{\parallel {\psi }_n-{\psi }_m\parallel }_{H^2}\le \hfill \\ \hfill & 2\epsilon CMR\left[(\mathrm{Pe}+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}C+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}{C}^2){\parallel T_{n-1}-T_{m-1}\parallel }_{H^2}\right.\hfill \\ \hfill & +(\mathrm{Re}+\mathrm{Re}C){\parallel {\omega }_{n-1}-{\omega }_{m-1}\parallel }_{H^2}\hfill \\ \hfill & \left.(\mathrm{Pe}+\mathrm{Re}+\mathrm{Re}C+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}C+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}{C}^2){\parallel {\psi }_{n-1}-{\psi }_{m-1}\parallel }_{H^2}\right].\hfill \end{array}\end{array}$$

Now, define

$$A=\mathrm{Pe}+\mathrm{Re}+\mathrm{Re}C+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}C+\mathrm{Pe}\mathrm{Ri}\mathrm{Re}{C}^2$$

such that

$$\begin{array}{cc}\hfill & \parallel T_n-T_m{\parallel }_{H^2}+\parallel {\omega }_n-{\omega }_m{\parallel }_{H^2}+{\parallel {\psi }_n-{\psi }_m\parallel }_{H^2}\le \hfill \\ & 6A\epsilon CMR(\parallel T_{n-1}-T_{m-1}{\parallel }_{H^2}+\parallel {\omega }_{n-1}-{\omega }_{m-1}{\parallel }_{H^2}+\parallel {\psi }_{n-1}-{\psi }_{m-1}{\parallel }_{H^2})\hfill \end{array}$$

(29)

where we note that $A\lesssim \mathrm{Pe}\mathrm{Ri}\mathrm{Re}$. Setting $\epsilon <\frac{1}{6ACMR}$ and defining $F:S\to S$ by $F\left(X_{(n-1)}\right)=X_{\left(n\right)}$ for functions $X\in S$ then ensures the Lipschitz condition such that for $X_{\left(n\right)}=(T_n,{\omega }_n,{\psi }_n)$:

$$d(F\left(X_{(n-1)}\right),F\left(Y_{(n-1)}\right))=d(X_{\left(n\right)},Y_{\left(n\right)})\le d(X_{(n-1)},Y_{(n-1)}).$$

(30)

Hence, the iteration is a contraction mapping.   □

Since $F:S\to S$ is a contraction mapping, we thus conclude by the Banach fixed-point theorem that F converges to a unique fixed-point $X^{*}$ such that the sequence $X_{\left(0\right)},X_{\left(1\right)},...,X_{\left(n\right)}$$\to X^{*}$. It then follows that this $X^{*}$ is the steady solution to the stream function-vorticity formulation of the Boussinesq Navier–Stokes equations for small amplitude forcing.   □

Since the sequence $X_{\left(0\right)},X_{\left(1\right)},...,X_{\left(n\right)}$ is a convergent, Cauchy sequence, it follows that the inequality of the form (30) is equivalent to the inequality

$$d(X_{(n+1)},X^{*})\le d(X_{\left(n\right)},X^{*})$$

(31)

which implies linear convergence to the unique solution $X^{*}$ for $\epsilon$ which satisfies the Lipschitz condition. From this result, we conclude that the iteration defined by (9) guarantees linear convergence for a sufficiently small amplitude forcing. The following section considers the relationship between the laser amplitudes for convergent solutions and the Re, Ri, and Pe numbers.

4. Numerical Results

4.1. Maximum Laser Amplitude

Based on the condition that the laser amplitude must be small enough for the iteration (9) to converge, we are motivated to find the relationship between the maximum allowable laser amplitude and the nondimensional fluid parameters Re, Ri, and Pe. We first consider a lemma which describes the asymptotic size of solutions as a function of these numbers.

Lemma 5

(Asymptotic Size of Solutions). Let $\parallel ·\parallel$ be the $H^2(\Omega )$ norm and let f be normalized such that ${\parallel f\parallel }_{L^2}=1$. If $Re\gg Pe$, then for $n\ge 1$,

$$\begin{array}{c}\hfill \parallel {\omega }_n\parallel =\mathcal{O}\left({Re}^{2^{n+1}-1}{Pe}^{2^n}{Ri}^{2^n}{\epsilon }^{2^n-1}\right),\end{array}$$

(32a)

$$\begin{array}{c}\hfill \parallel T_n\parallel =\mathcal{O}\left({Re}^{2(2^n-1)-n}{Pe}^{2^n+n}{Ri}^{2^n-1}{\epsilon }^{2^n-1}\right).\end{array}$$

(32b)

Similarly, if $Pe\gg Re$, then

$$\begin{array}{c}\hfill \parallel {\omega }_n\parallel =\mathcal{O}\left({Re}^{2^n}{Pe}^{2^{n+1}-1}{Ri}^{2^n}{\epsilon }^{2^n-1}\right),\end{array}$$

(33a)

$$\begin{array}{c}\hfill \parallel T_n\parallel =\mathcal{O}\left({Re}^{2^n-1}{Pe}^{2^{n+1}-1}{Ri}^{2^n-1}{\epsilon }^{2^n-1}\right).\end{array}$$

(33b)

Proof. 

The proof proceeds by a simple inductive argument on the iteration (9).   □

We now introduce a theorem which identifies an upper bound for laser amplitude as an asymptotic of the Re, Ri, Pe numbers.

Theorem 3

(Upper Bound for Laser Amplitude). Let E be the set of all amplitudes ε such that the iteration (9) converges. If $Re\gg Pe$, then

$${\epsilon }_{\sup }=supE=\mathcal{O}\left({Re}^{-2}{pe}^{-1}{Ri}^{-1}\right).$$

(34)

Similarly, if $Pe\gg Re$, then

$${\epsilon }_{\sup }=supE=\mathcal{O}\left({Re}^{-1}{pe}^{-2}{Ri}^{-1}\right).$$

(35)

Proof. 

Recall from Theorem 2 that for the fixed-point iteration to converge, we require

$$\epsilon <\mathcal{O}\left(A^{-1}{R}^{-1}\right)$$

for an $R\left(\epsilon \right)$ such that the entire sequence of solutions stay bounded in a ball of size R. A least upper bound for the allowable laser amplitude, then, may be determined by finding a greatest lower bound for R such that the sequence of solutions stays bounded. This $R_{\inf }$ will thus be a least upper bound for solutions in the $H^2(\Omega )$ norm for a given $\epsilon$. Suppose that $\mathrm{Re}\gg \mathrm{Pe}$. Then from the above lemma, we observe that the size of the nth iteration may be bounded by the asymptotic $\mathcal{O}\left({\mathrm{Re}}^{2^{n+1}-1}{\mathrm{Pe}}^{2^n}{\mathrm{Ri}}^{2^n}{\epsilon }^{2^n-1}\right)$. Then for the supremum of the amplitude,

$$\begin{array}{cc}\hfill {\epsilon }_{\sup }& =\mathcal{O}\left(A^{-1}{R}_{\inf }^{-1}\right)\hfill \\ \hfill & =\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{-1}{\mathrm{Pe}}^{-1}{\mathrm{Ri}}^{-1}{\left({\mathrm{Re}}^{2^{n+1}-1}{\mathrm{Pe}}^{2^n}{\mathrm{Ri}}^{2^n}{\epsilon }_{\sup }^{2^n-1}\right)}^{-1}\right)\hfill \\ \hfill & =\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{2^{n+1}}{\mathrm{Pe}}^{2^n+1}{\mathrm{Ri}}^{2^n+1}{\epsilon }_{\sup }^{2^n-1}\right){,}^{-1}\mathrm{so}\hfill \end{array}$$

$${\epsilon }_{\sup }^{2^n}=\mathcal{O}\left({\mathrm{Re}}^{-2^{n+1}}{\mathrm{Pe}}^{-2^n-1}{\mathrm{Ri}}^{-2^n-1}\right)\mathrm{and}\mathrm{thus}$$

$${\epsilon }_{\sup }=\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{-2}{\mathrm{Pe}}^{-\frac{2^n-1}{2^n}}{\mathrm{Ri}}^{-\frac{2^n-1}{2^n}}\right)=\mathcal{O}\left({\mathrm{Re}}^{-2}{\mathrm{Pe}}^{-1}{\mathrm{Ri}}^{-1}\right).$$

Now, suppose that $\mathrm{Pe}\gg \mathrm{Re}$. By again applying Lemma 5, we observe that the size of the nth iteration may now be bounded by the asymptotic $\mathcal{O}\left({\mathrm{Re}}^{2^n}{\mathrm{Pe}}^{2^{n+1}-1}{\mathrm{Ri}}^{2^n}{\epsilon }^{2^n-1}\right)$. Then for the supremum of the amplitude,

$$\begin{array}{cc}\hfill {\epsilon }_{\sup }& =\mathcal{O}\left(A^{-1}{R}_{\inf }^{-1}\right)\hfill \\ \hfill & =\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{-1}{\mathrm{Pe}}^{-1}{\mathrm{Ri}}^{-1}{\left({\mathrm{Re}}^{2^n}{\mathrm{Pe}}^{2^{n+1}-1}{\mathrm{Ri}}^{2^n}{\epsilon }_{\sup }^{2^n-1}\right)}^{-1}\right)\hfill \\ \hfill & =\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{2^n+1}{\mathrm{Pe}}^{2^{n+1}}{\mathrm{Ri}}^{2^n+1}{\epsilon }_{\sup }^{2^n-1}\right){,}^{-1}\mathrm{so}\hfill \end{array}$$

$${\epsilon }_{\sup }^{2^n}=\mathcal{O}\left({\mathrm{Re}}^{-2^n-1}{\mathrm{Pe}}^{-2^{n+1}-1}{\mathrm{Ri}}^{-2^n-1}\right)\mathrm{and}\mathrm{thus}$$

$${\epsilon }_{\sup }=\underset{n\in \mathbb{N}}{inf}\mathcal{O}\left({\mathrm{Re}}^{-\frac{2^n-1}{2^n}}{\mathrm{Pe}}^{-2}{\mathrm{Ri}}^{-\frac{2^n-1}{2^n}}\right)=\mathcal{O}\left({\mathrm{Re}}^{-1}{\mathrm{Pe}}^{-2}{\mathrm{Ri}}^{-1}\right).$$

  □

Based on the relationships defined in Theorem 3, we expect the maximum computable amplitude in a numerical implementation of the iteration to adhere to an asymptotic relationship in these dimensionless parameters. By considering the homogeneous Dirichlet boundary condition in the stream function-vorticity formulation, we run the algorithm for increasing laser amplitudes until the iteration diverges to determine the max amplitude within an error of $1\%$. Figure 2 and Figure 3 provide these numerical results with a comparison to the expected asymptotic relationship.

From these results, we observe that the convergence analysis gives accurate predictions for the maximum laser amplitude as computed in the numerical algorithm. This is significant, since we see that the maximum allowable laser amplitude scales by the inverse square of the Reynolds number or Peclet number depending on their relative size. For physical laser-fluid setups, this result suggests that the computation of steady solutions is particularly sensitive to fluid viscosity and fluid thermal diffusivity. Physical intuition suggests that this is true, since less viscous flows with lower rates of heat transfer will have a more elusive steady state. Hence, finding convergent solutions for both large amplitude laser forcing and large Reynolds and Peclet numbers proves difficult.

4.2. Numerical Simulations

We now consider the convergent solutions found through a numerical implementation of the iteration (9) in the slip boundary condition. First, we specify the Re, Ri, Pe numbers to appropriately characterize fluid flow. Akers and Reeger [1] provide a discussion on the relationship between these parameters and the fluid behavior in thermal blooming. For the normalized laser forcing $f(x,y)=e^{-c(x^2+y^2)}$, we pick $c=5$ to represent rapid Gaussian decay from the origin towards the boundary of the domain. This parameter, however, can also be changed based on the exact laser specifications. We assign $\mathrm{Re}=\mathrm{Pe}=10$ and $\mathrm{Ri}=100$ with a grid spacing of $h=\frac{1}{50}$ and a laser amplitude of $\epsilon ={10}^{-3}$. Figure 4 provides a streamline plot to the left and a plot of the temperature fluctuation to the right for the steady-state solution to the slip boundary condition.

We can clearly see here the result of the buoyancy driven convection within the domain. Fluid heated by the laser in the center of the domain rises until it contacts the top, then travels along the boundary to satisfy the slip condition. Steady vortex rings are formed to the left and right of center, which are centered horizontally at the halfway point between the middle and the edge of the domain. The vertical center of the vortex rings shift with the peak temperature fluctuation, which occurs just above the domain center due to the heated fluid rising in a buoyancy driven flow. Figure 5 provides a plot of the temperature fluctuation along the $x=0$ line as a function of the vertical position.

Due to the influence of gravity in the steady-state, the temperature fluctuation exhibits asymmetry in the vertical direction. The peak temperature is shifted slightly upwards in the domain before asymmetric decay in temperature forces temperature to zero along the boundary.

5. Conclusions

This study considered the existence of solutions to a 2D steady, Boussinesq fluid flow problem in a stream function-vorticity formulation. Motivated by the fact that thermal blooming often plays an immense role in the simulation of laser propagation, the algorithm introduced in this paper provides a rapid numerical method for computing the steady-state of a buoyancy driven fluid while heated by a laser. For a small enough laser amplitude, we showed with a fixed-point argument that steady solutions indeed exist in an appropriate function space. The threshold which determines how small this amplitude must be is related to the nondimensional fluid parameters Re, Ri, Pe through an established asymptotic relationship. The theoretical derivation of this relationship was corroborated by numerical simulations, with a key result being that high amplitude convergent solutions are difficult to find for large Re and Pe numbers. We showed that the slip boundary condition for fluid velocity is equivalent to a homogeneous Dirichlet condition on the stream function and vorticity based on their respective definitions. The numerical simulations of these steady solutions described the existence of symmetric vortex rings in the left and right sides of the domain, implying that the fluid heated by the laser rises in a circulatory fashion before being recirculated within the domain. The peak temperature fluctuation was shifted upwards in the domain due to the same convective behavior, which prognosticates asymmetric laser propagation when slip boundary conditions are imposed in a finite box. In building upon this work, the authors are seeking to couple the presented fluid iteration with a laser propagation model to determine how buoyancy driven thermal blooming impacts on-target laser irradiance. Included in this pursuit is the investigation of more apt boundary conditions around the domain for a free-space laser propagation setup.