Using boundary-layer coordinates to compute unsteady laminar two-dimensional viscous and convective flows

In computational fluid dynamics, it is often necessary to implement a variable grid in order to adequately resolve the flow field. There are many techniques available to achieve this. For example, in one dimension, we can discretize the interval 0 ≤ Z ≤ H into M unequal subintervals by positioning the grid points, Z_i, such that

with

This discretization scheme will space the grid points closer together near the bottom (Z = 0) when the parameter β < 0, whereas for β > 0, the grid points will be clustered near the top (Z = H). This is a computational approach since the governing equations are not altered in any way. In some situations, this approach may be sufficient. However, for unsteady boundary-layer flows, this scheme will be inadequate because the variable grid spacing needed must be able to vary with time, t.

In this study, a coordinate transformation that is ideal for computing flows possessing two unsteady boundary layers in different directions is presented. The transformation leads to a new mathematical formulation of the problem that accounts for the inherent evolving boundary-layer structure of the flow and, thus, has the effect of implementing a variable grid spacing that changes with time. The method proposed is an alternative to the formula listed above. Essentially, the approach outlined here represents an analytical approach to performing grid reduction rather than a computational approach.

The proposed method has already been successfully applied to compute various flows around cylinders.^1–7 In these flows, a single boundary layer forms along the surface; hence, introducing a boundary-layer coordinate orthogonal to the surface makes a lot of sense. In this study, we will explore flows possessing two interacting boundary layers, and to properly describe these flows, we need to implement two boundary-layer coordinates. Thus, this work will focus on extending the concepts used in the previous studies^1–7 to flows, which naturally require the use of boundary-layer coordinates in two dimensions. Examples of such flows include flow around a corner, free convective flow from a heated corner, and mixed convective vortex flow. In addition, a numerical solution procedure to solve the transformed equations is also proposed.

Corner flow occurs when fluid is trapped between two intersecting planes. There are several different problems related to corner flow; for example, for two-dimensional steady creeping flow near a corner, there is the so-called paint-scraper problem solved by Taylor⁸ where one plane moves parallel to itself with constant velocity, and there is also the hinged-plate problem solved by Moffatt⁹ where the planes are rotating toward or away from each other about their line of intersection with constant angular velocity. Furthermore, if the planes are stationary and the motion is driven by far-field stirring, then Moffatt showed that there exists an infinite sequence of reversing eddies of decreasing size provided the corner angle is less than 73°. In addition, there is also the three-dimensional problem of flow along a right-angle corner. This problem was originally solved within the framework of laminar boundary-layer theory by Carrier¹⁰ and later by Rubin and Grossman.¹¹ In these steady-state investigations, the following coordinate transformation was used:

Here, the flow is along the x direction with U denoting the uniform far-field velocity, ν is the kinematic viscosity, and ξ and η are boundary-layer coordinates. In this work, we will propose an alternate choice of boundary-layer coordinates, which are suitable for unsteady two-dimensional flows. A recent study involving corner flow was carried out by Hinton, Hogg, and Huppert.¹² They investigated the steady gravity-driven free-surface viscous flow down an inclined plane around a corner using lubrication theory. They showed that far downslope after detachment, the motion is accurately modeled by a similar solution.

The corner flow problem considered in this study regards the unsteady two-dimensional flow around a right-angle corner, and here, we will solve the fully nonlinear Navier–Stokes equations instead of the simplified creeping flow or boundary-layer equations. To our knowledge, this problem has not yet been solved. However, a well-known closely related problem is the steady two-dimensional stagnation flow first studied by Hiemenz¹³ (also referred to as Hiemenz flow¹⁴) and later modified by Howarth.¹⁵ Here, a uniform flow along the y axis impinges on a flat wall located at y = 0. The flow takes place in the upper x–y plane and the stagnation point, located at the origin, divides the flow into two streams along the wall, which then depart in opposite directions: one in the positive x direction and the other in the negative x direction. The case of unsteady flow was later investigated by Rott,¹⁶ and a class of similar solutions pertaining to stagnation point flows having a far-field accelerating or decelerating incoming flow was discovered by Kolomenskiy and Moffatt.¹⁷ 

Free convection from heated plates (vertical and horizontal) is well documented in Refs. 14 and 18. The specific problem considered in this paper is the unsteady two-dimensional free convective laminar flow from an unbounded heated right-angle corner. The flow starts from rest and is set into motion by a buoyancy force, which is triggered by a constant heat flux imposed along the walls. A problem similar to that discussed in this study involves thermally driven laminar flows in bounded cavities, and some previous studies include it.^19–23 Although the coordinate transformation proposed in this study is best suited for unbounded domains, it can be applied to a bounded domain by using the method of domain decomposition.²⁴ This involves decomposing the domain into two or more regions to properly resolve the flow. The boundary-layer coordinate transformation can then be implemented in the regions adjacent to the boundaries.

The last example discussed is that of mixed convection enhanced by vortex flow. Vortex flows can be described as swirling motions or whirlwinds, and because of their association with geophysical phenomena, such as dust devils and tornadoes, numerous studies have been devoted to understanding them including two review articles.^25,26 Early investigations^27,28 have approached vortex flows through the lens of boundary-layer theory. Some numerical simulations of laminar rotating flows include the research reported in Refs. 29 and 30, while experimental studies are reported in Refs. 31–33. The problem solved here involves an unsteady axisymmetric laminar thermal flow formulated in a rotating reference frame and heated from below and is motivated by the studies.^34,35 We will see that even though there is only one wall located at z = 0, the boundary-layer coordinates still apply since boundary layers will form along the axis of rotation as well as along the wall for rotating flows. A related problem is a flow near a rotating disk, which is described in Refs. 14 and 36.

This paper is structured as follows: the next three sections focus on computational examples of the proposed method. In Sec. II, the flow around a two-dimensional corner is discussed. Following that, in Sec. III, two-dimensional free convective flow from a heated corner is presented. Then in Sec. IV, mixed convective vortex flow is addressed under the assumption of azimuthal symmetry. A brief summary is provided in Sec. V. Finally, three appendices are included, which present comparisons with previous studies as a means of validating the proposed formulation.

We begin with a straightforward example to illustrate the idea and mathematical formulation. This involves the flow past a 90° corner shown in Fig. 1. Here, we consider the unsteady two-dimensional laminar flow of a viscous incompressible Newtonian fluid around the corner in the first quadrant. We will assume that the flow remains laminar and two dimensional for all time. The governing Navier–Stokes and continuity equations cast in the dimensionless form are given by

In the above, u and v denote the velocities in the x, y directions, respectively, while P is the pressure. The parameter Re refers to the familiar Reynolds number defined by Re = UL/ν, where U and L are velocity and length scales, respectively, and ν is the kinematic viscosity. Equations (1)–(3) must be solved for x, y > 0 under the no-slip boundary conditions of u = v = 0 along the walls of x = 0 and y = 0, and the far-field conditions illustrated in Fig. 1.

To understand the origin of the proposed coordinate transformation, we examine the dominant terms in Eq. (1) for large x. In this limit, we can discard derivatives with respect to x and set v ≈ 0. Then, Eq. (1) can be approximated by

and the solution satisfying u = 0 at y = 0 and u → 1 as y → ∞ is easily found to be

Applying the same argument for large y, we obtain a similar result,

The parameter λ is a measure of the boundary-layer thickness formed along the axes and suggests the introduction of boundary-layer coordinates ξ and η, which are defined by

An interpretation of working in terms of boundary-layer coordinates ξ and η is that the physical coordinates x and y become moving coordinates; that is, lines of constant ξ and η expand in time when plotted in a Cartesian coordinate system. This is ideal from a numerical point of view since the grid lines are alive and allowed to expand with the growing boundary layers to ensure adequate resolution with the passage of time; for example, lines of constant ξ when plotted in the Cartesian plane represent vertical lines that are initially clustered near x = 0 and spread apart from each other with time at a rate that is consistent with the growth of the boundary layer. This is particularly important in the early development of the flow when the boundary layer grows more rapidly. At later times when the boundary layer thickens appreciably, it may make sense to switch back to the original Cartesian coordinates. An optimal time to do so would be at time t = Re/4 since λ = 1 at this instant, and thus, x = ξ and y = η. To illustrate the evolving grid, Fig. 2 shows the case with Re = 100 at times t = 1 and t = 25. Here, the closest five vertical and horizontal grid lines to the origin are plotted in physical coordinates at the two times. The uniform grid spacing h = Δξ = Δη = 0.2 was used, and note that when t = 25, λ = 1 for Re = 100.

In terms of ξ and η and noting that

Eqs. (1)–(3) transform to

Since the flow is two dimensional, it is worth working in terms of a stream function (ψ) - vorticity (ω) formulation. This is accomplished by defining

which automatically satisfies the continuity equation, and

System (5)–(7) then simplifies to

Before proceeding to solve this system, we rescale ψ and ω according to

Then,

and Eqs. (8) and (9) become

To solve systems (10) and (11), we first note that the initial conditions for Ψ and Ω, denoted by Ψ₀ and Ω₀, can be obtained by setting t = 0 in (10) and (11) and thus satisfy

Since Eqs. (12) and (13) form a linear coupled system, they can easily be solved numerically using finite differences. This involves replacing all derivatives with second-order, central-difference approximations and then numerically solving a large linear system of algebraic equations. The only complication that arises stems from the fact that the surface vorticity is not known. That is, the vorticity is underspecified while the stream function is overspecified. We can overcome this difficulty by deriving a second-order accurate formula for the surface vorticity as follows: We present the details for the boundary η = 0 with the understanding that the boundary ξ = 0 can be obtained using a similar procedure. Along η = 0, the no-slip conditions in terms of the stream function become

In addition, from (10), we have that

Expanding Ψ in a Taylor series about η = 0 yields

where h is the uniform grid spacing in both the ξ and η directions. Using the known boundary conditions for Ψ together with the expressions for the higher derivatives yields the following formula for the surface vorticity:

In arriving at this result, we have made use of the second-order finite difference approximation,

The far-field boundary conditions for Ψ and Ω along ξ_∞ can be derived by making use of u → erf(η), v → 0, from which we obtain

for the vorticity, while the stream function will satisfy

Similar conditions apply along the boundary η_∞.

Equations (10) and (11) were solved on the square domain [0, ξ_∞] × [0, η_∞] with ξ_∞ = η_∞ = 5 using the boundary conditions discussed above along the edges. The domain was discretized into N = 100 equally spaced subintervals in both the ξ and η directions having a grid spacing of h = 0.05. Adopting the notation $fi,jk=f(ξi,ηj,tk)$ and discretizing Eq. (10) using second-order central differences lead to

To solve Eq. (11), we begin by rewriting it in the form

where

Assuming the solution at time t is known, we can advance the solution to time t + Δt by integrating (16) and using integration by parts as follows:

where Δt is the time increment. We now approximate the integrals using

where ϖ is a weight factor and χ is a generic function. In general, 0 ≤ ϖ ≤ 1 and we treat ϖ as a computational parameter. Here, we have used ϖ = 1/2, which yields the well-known Crank–Nicolson scheme. With this approximation in place, we obtain

Upon substituting the expression for R(ξ, η, t + Δt) and replacing all derivatives using central-difference approximations, Eq. (17) becomes a nonlinear system of algebraic equations, which was solved in conjunction with (15) using the Gauss–Seidel iterative procedure to determine both Ψ and Ω at time t + Δt at all the grid points. In our computations, Δt = 0.005 was used, and no convergence problems were encountered in our numerical scheme with the parameters listed.

We first present the initial solution that is independent of the Reynolds number. This solution was used as an initial condition for obtaining unsteady solutions for a given Reynolds number. In Fig. 3, the surface vorticity distribution along the walls ξ = 0 and η = 0 is plotted. The distributions along the walls are identical, suggesting that the flow is initially symmetrical. This is confirmed in the streamline plot shown in Fig. 4, which is also plotted in boundary-layer coordinates. We notice a region of clockwise circulating flow near the corner, and based on the magnitude of the values of the streamlines, the flow is nearly stagnant as it is rotating slowly. This explains why the surface vorticity profile increases slowly with distance from the origin before attaining its maximum value. Afterward, it decreases monotonically to its asymptotic value. The flow near the walls is largely tangential to the walls and reverses direction; within the circulating eddy region, the tangential flow is toward the corner, while outside the eddy, the flow is directed away from the corner. The maximum value in vorticity occurs where the transverse gradient in the tangential flow is the largest, and the location of the maximum indicates that this occurs near the middle of the eddy.

To ensure that the flow remains laminar and two dimensional, we limit our results to Re = 100, although solutions for larger Re can easily be obtained. Simulations were carried out for 0 ≤ t ≤ 25. As time increased, a larger computational domain was needed to fully capture the growing eddy. In Figs. 5 and 6, surface vorticity distributions and streamlines for Re = 100 at t = 25 in physical coordinates x and y are shown. These figures clearly illustrate that the flow no longer remains symmetrical with the passage of time. The eddy size increases significantly with time as does the strength of the circulation, and this is reflected in the surface vorticity distribution. The vorticity profile plotted in Fig. 5 shows a more rapid increase, leading to a much larger maximum value than that shown in Fig. 3. The maximum value is larger along the wall x = 0 because the shape of the eddy is such that it stretches further in the y direction than it does in the x direction, and so the gradient in tangential flow is larger along the wall x = 0 than it is along the wall y = 0. This also explains why the maxima occur at different locations. The profiles also take on other extremes before assuming their asymptotic values. Again, this is because of the variation in the gradient of the tangential flow along the walls.

Since stagnation flow is closely related to our problem, it was used to validate our formulation. The comparison between our findings against known results for stagnation flow is presented in  Appendix A. It is shown that the agreement is good.

We next consider the unsteady laminar two-dimensional free convective flow induced by a heat flux. The fluid is initially at rest at a constant temperature, and the flow is confined to the first quadrant. At t = 0 an impulsive heat flux is applied along the walls x = 0 and y = 0 with gravity acting in the vertical direction. The equations governing the motion and heat transfer process of a viscous incompressible fluid are the continuity, Navier–Stokes, and energy equations. The fluid is characterized by the following properties: the kinematic viscosity, ν, the thermal diffusivity, κ, the thermal expansion coefficient, α, and the thermal conductivity, k. While these fluid properties are assumed to remain constant, the fluid density, ρ, is allowed to vary with temperature in the usual linear fashion given by

where T denotes the temperature while T₀ is the far-field constant ambient temperature and ρ₀ is the corresponding density.

In Cartesian coordinates, the equations in dimensionless form become

Here, the Boussinesq approximation was used to model the driving buoyancy force, and viscous dissipation was omitted. The parameters Gr and Pr are known as the Grashof and Prandtl numbers, respectively, and are defined by

where ΔT = LQ₀/k with Q₀ representing the heat flux scale, and L is the length scale. The velocity scale used in rendering the equations in dimensionless form was $U=αgLΔT$⁠, where g is the acceleration due to gravity, and the dimensionless temperature, ϕ, is related to the dimensional temperature, T, through the relation ϕ = (T − T₀)/ΔT.

Since the flow is assumed to remain two dimensional for all time, it makes sense to work in terms of a stream function and vorticity as previously defined. Then, Eqs. (18)–(21) reduce to

Using a similar argument as in Sec. II, the boundary-layer coordinates in this case become

Here, μ is a measure of the thermal boundary-layer thickness. The parameter μ is very similar to the parameter λ from Sec. II with $Gr$⁠, playing the role of the Reynolds number, Re. Even though the fluid is initially at rest, significant temperature gradients will form along the walls as a result of the applied heat flux and, hence, the need for scaled coordinates. Thus, the methodology discussed in Sec. II for velocity gradients can easily be extended to include temperature gradients. In terms of these boundary-layer coordinates, Eqs. (22)–(24) transform to

and need to be solved subject to the quiescent far-field conditions ψ, ϕ, ω → 0 as ξ² + η² → ∞ with ξ, η > 0. Along the wall η = 0, we enforce no-slip, and the applied heat flux is given by

while along ξ = 0, we have the similar conditions

Here, Q is the applied heat flux, which is taken to vary with position along the walls. In our simulations, we have set $Q(s)=e−s2$⁠. Thus, the heat flux is concentrated near the origin where the walls meet. Since the fluid is initially at rest with uniform temperature T₀, the initial conditions are simply ψ₀ = ω₀ = ϕ₀ = 0. To determine the surface vorticity, we follow the approach discussed in Sec. II. Along the wall η = 0, the following expression is easily obtained:

where again, h denotes the uniform grid spacing. A similar expression results along the wall ξ = 0. From the heat flux condition, the following second-order accurate formula was derived to compute the surface temperature:

along

and a similar expression was derived to obtain the temperature along the wall ξ = 0.

Since system (26)–(28) bears a close resemblance to systems (10) and (11), the same numerical scheme was used to integrate (26)–(28), including the same grid size, h, and time step, Δt. Here, results for the case Gr = 100 and Pr = 0.7 (for air) are presented. Figures 7–10 illustrate that the surface vorticity, surface temperature, streamline, and isotherm plots, respectively, at t = 2.5 and are shown in physical coordinates. Both the surface vorticity and surface temperature profiles are nearly symmetrical, which is indicative of the fact that during the early stages, the dominant mode of heat transfer is conduction. This is further supported by the isotherm plot, which reveals almost circular isotherms as a result of conduction in the radial direction. The streamline plot shows a slowly rotating pattern that is circulating in the clockwise direction.

In Figs. 11–14, the surface vorticity, surface temperature, streamline, and isotherm plots, respectively, at a later time t = 25 are shown again in physical coordinates. Here, we see noticeable asymmetry as a result of convection, which is intensifying due to the rising motion of fluid caused by buoyancy. The isotherm plot reveals a mushroom pattern, which is a signature of convective flows and has been reported in previous studies such as Refs. 5 and 35. The clockwise circulating flow is much stronger as noted by the magnitude of the values of the streamlines. The surface vorticity and surface temperature profiles are significantly larger along the wall x = 0 than they are along y = 0. This can be explained as follows: as the fluid adjacent to the wall x = 0 is heated, it will rise and, therefore, increase the fluid temperature near the wall above as it drifts. The rising fluid motion creates larger velocity gradients along the wall x = 0 than it does along the wall y = 0. Hence, the vorticity will be larger along x = 0, especially in the vicinity of the corner where the heat flux is the largest. The flatter section in the temperature profile shown in Fig. 12 from y ≈ 3 to 6 is consistent with the wider spacing between isotherms shown in Fig. 14 over that range of y values.

Finally, to validate our formulation, we made a comparison with the closely related problem of free convection from a vertical hot wall. This comparison is carried out in  Appendix B, and good agreement was found.

We now consider the axisymmetric unsteady rotating laminar flow of a viscous incompressible fluid with a variable density that is heated from below. Expressed in cylindrical coordinates (r, θ, z), the nondimensional Navier–Stokes equations in a rotating reference frame, energy, and continuity equations are given by

In the above, u, v, and w are the radial, azimuthal, and axial velocity components, respectively, P is the pressure, f is the rate of rotation, which, in geophysical flows, is known as the Coriolis parameter, and Pr and Re are the Prandtl and Reynolds numbers, respectively. The parameter Ra is the Rayleigh number, which is defined by

with T_w > T₀ denoting the constant temperature maintained along the bottom wall z = 0, T₀ is the initial surrounding constant ambient temperature, and the remaining quantities have been previously defined. The dimensionless temperature is again denoted by ϕ, but this time it is related to the dimensional temperature T through the relation ϕ = (T − T₀)/(T_w − T₀). The density was allowed to vary linearly with temperature as in Sec. III.

Before introducing the boundary-layer coordinates, we first define new velocity components,

and also define the stream function, ψ, in the r–z plane and the azimuthal component of vorticity, ω, as follows:

Then, Eqs. (31)–(35) reduce to

where ζ = rω.

For rotating flows, thin boundary layers will develop along the bottom wall z = 0 and along the axis of rotation r = 0. Thus, it makes sense to once again introduce the boundary-layer coordinates,

and to rescale the stream function and vorticity according to

Then, systems (36)–(39) transform to

System (40)–(43) is solved on the domain ξ, η > 0. The initial solutions, Ψ₀, ω₀, V₀, and ϕ₀, can be obtained by setting t = 0 and thus satisfy

The exact solution to (45) satisfying the boundary conditions V₀(ξ = 0, η) = V₀(ξ, η = 0) = 0 with Re = 0 is given by

where Γ denotes the circulation at infinity and will be treated as an additional parameter controlling the flow. It is also referred to as the swirl ratio. Also, the exact solution to (47) satisfying the following conditions:

with Re = 0 is found to be

and indicates that the isotherms are initially horizontal lines. The initial solutions Ψ₀ and Ω₀ with Re = 0 can be determined by numerically solving (44) and (46) using the known solution for V₀. These equations were solved on the square domain [0, ξ_∞] × [0, η_∞] with ξ_∞ = η_∞ by finite differences using the scheme outlined in Sec. II to solve (12) and (13). Along the bottom wall η = 0, the no-slip conditions were applied,

while along the axis of rotation ξ = 0, the following conditions were imposed:

We note that the above conditions also apply for all t > 0. To determine the surface vorticity along the bottom wall η = 0, a second-order accurate formula was derived using the approach outlined in Sec. II. Along the outer boundaries ξ_∞ and η_∞, conditions Ψ = Ω = 0 were used for all t ≥ 0. The validity of using these conditions will be addressed shortly.

The flow is completely characterized by the parameters Re, Ra, f, Γ, and Pr. The aim here is not to conduct an extensive parameter investigation but rather to thoroughly present some key results for selected values. With this in mind, we have fixed the values of Pr, Ra, and Re to be Pr = 0.7, Ra = 1, and Re = 25 while f and Γ were allowed to take on the values f = 0, 0.05 and Γ = 0, 0.5 to demonstrate their influence on the flow. This choice of values ensures the flow remains laminar and two-dimensional for all time. Unsteady calculations were carried out by numerically integrating systems (40)–(43) using the method from Sec. II with Δt = 0.005 and h = 0.05. The simulations were performed over the time interval 0 < t ≤ 25, and no convergence problems were encountered with these parameter values.

We begin with the case f = Γ = 0 (Re = 25, Ra = 1, Pr = 0.7), which we refer to as the base case since it amounts to free convective flow with the Reynolds number, assuming the role of the Grashof number according to the relation $Re∼2Gr14$⁠. As a result of f = Γ = 0, the azimuthal velocity remains zero for all time. The plots to be presented are all at t = 25. In Fig. 15, the surface vorticity profile is shown, while Fig. 16 illustrates the temperature variation along the axis r = 0. In Fig. 17, the local surface heat transfer coefficient is displayed, which is defined by

also known as the Nusselt number. Figures 18–20 illustrate contour plots for streamlines, isotherms, and vorticity, respectively. The vorticity contour plot reveals that vorticity is generated along the bottom wall and is confined to a thin region adjacent to the wall. The streamline plot, on the other hand, illustrates the formation of a large central rotating toroidal cell flanked by two smaller counter-rotating toroidal cells. The counter-rotating cells are consistent with the sign reversals in the surface vorticity. Also, the maximum value of the surface vorticity aligns well with the convergence of streamlines shown near the bottom wall of Fig. 18. The isotherms are no longer horizontal, and the spacing between them is in accordance with the heat transfer coefficient shown in Fig. 17. That is, the maximum in the heat transfer coefficient occurs where the isotherms are closer together, and the minimum happens when the isotherms are furthest apart. The axial temperature profile shows a rapid decrease along the axis.

We next consider the case f = 0 and Γ = 0.5. Since Γ ≠ 0, there is a noticeable change in the azimuthal velocity as seen in the contour plots of Figs. 21 and 22. Figure 21 is the contour plot corresponding to Eq. (48) plotted in boundary-layer coordinates while Fig. 22 is obtained from our numerical simulation at t = 25. The waviness observed in the contour lines of Fig. 22 is due to convection. Both of these plots reveal rapid variations in the azimuthal velocity along the bottom wall and the axis of rotation, while in the interior of the domain, there is little variation. The contour plots for streamlines, isotherms, and vorticity displayed very similar features as those illustrated in Figs. 18–20.

In the last case considered, we set f = 0.05 and Γ = 0.5. Here, we see significant changes in the contour plot for V shown in Fig. 23. Due to the rotating reference frame, there is a significant change in the azimuthal velocity, which contributes to the tightly packed contour lines seen in the diagram. We also notice more variation in V throughout the entire domain when compared to Fig. 22. The isotherm, streamline, and vorticity contour plots showed the same features as those already observed and discussed in Figs. 18–20. Finally, in Figs. 24–26, the surface vorticity, axial temperature, and heat transfer profiles for the three cases are contrasted. While there is little change in the surface vorticity distributions, there are significant differences in the axial temperature and heat transfer coefficient near the axis (r < 5). The impact of rotation and vortex flow is evident in the broadening of the axial temperature profile, and the net effect is to slightly enhance convection as indicated by the increase in the average value of the Nusselt number, denoted by $Nū$⁠. For the base case, $Nū=2.6$⁠, while for the other two cases, $Nū=2.7$⁠.

We next validate the far-field conditions used along the boundary ξ_∞. An asymptotic solution can be derived by arguing that along ξ_∞, Eq. (42) is well approximated by