**3. Integral Forms of Equations**

The solution of the stated problem will be performed using the integral method. The essence of the integral method is the reduction in the original system of differential transport equations in partial derivatives into a system of ordinary differential equations, in which the marching coordinate *x* remains the only independent variable [29,30]. To do this, one should assume proper approximating functions for the profiles of the longitudinal velocity component *u* and the local temperature *T* in the boundary layer, and integrate transport Equations (10) and (11) within the limits from the wall to the outer boundary δ of the boundary layer. In this case, the continuity equation is incorporated into the integral equation of momentum transfer (to exclude the transverse velocity component v). The unknown quantities in the so-called integral equations of the boundary layer (which, in fact, after integration across the boundary layer, become ordinary differential equations) are: (a) the boundary layer thickness δ, which for gases is the same for the velocity and temperature profiles; and (b) the scaling (maximum) velocity in the boundary layer *U*, which increases with the development of the boundary layer and is itself unknown, and which must be found from the solution of the problem. Therefore, the only boundary conditions along the marching coordinate *x*, which are needed for solving the integral equations of the boundary layer, are the equality to zero of the thickness of the boundary layer δ and the scaling velocity *U* at the leading edge of the plate at *x* = 0. These boundary conditions will be implicitly incorporated into the mathematical form of the functions δ(*x*) and *U*(*x*) in the subsequent solution of the problem.

Integral equations for the boundary layer in the case of the free convection of an ideal gas are given in the classical works [30,31]. Let us derive integral equations for the case of a real gas considered in this paper and described by the van der Waals equation. To obtain the integral equation of motion, we multiply the continuity Equation (9) by the velocity component and add it to the left side of Equation (19). As a result, we have:

$$\frac{\partial u^2}{\partial x} + \frac{\partial u\mathbf{v}}{\partial y} = \mathbf{v}\frac{\partial^2 u}{\partial y^2} + g\left(1 - \frac{\mathbf{Z}}{1 + \beta\Delta T\theta} \left(1 - \frac{\mathbf{W}\mathbf{a}\_b}{1 + \beta\Delta T\theta} + \frac{\mathbf{W}\mathbf{a}\_d}{\left(1 + \beta\Delta T\theta\right)^2}\right)\right),\tag{22}$$

Let us integrate Equation (22) over the thickness of the boundary layer, taking into account the boundary conditions (12) and (13). This gives:

$$\frac{d}{dx}\int\_{0}^{\delta}u^{2}dy=-\text{v}\left(\frac{\partial u}{\partial y}\right)\_{y=0}+\text{g}\int\_{0}^{\delta}\left(1-\frac{Z}{1+\mathsf{f}\Delta T\mathsf{\theta}}\left(1-\frac{\text{Wa}\_{\mathsf{b}}}{1+\mathsf{f}\Delta T\mathsf{\theta}}+\frac{\text{Wa}\_{\mathsf{d}}}{\left(1+\mathsf{f}\Delta T\mathsf{\theta}\right)^{2}}\right)\right)dy\,\tag{23}$$

Let us further rewrite this equation using dimensionless variables and functions:

$$\frac{d(\delta \mathcal{U}^2)}{dx} \int\_0^1 w^2 d\eta = -\mathbf{v} \left(\frac{\partial w}{\partial \eta}\right)\_{\eta=0} \frac{\mathcal{U}}{\delta} + g\delta \int\_0^1 \left(1 - \frac{Z}{1 + \beta \Delta T \theta} \left(1 - \frac{\text{Wa}\_b}{1 + \beta \Delta T \theta} + \frac{\text{Wa}\_d}{\left(1 + \beta \Delta T \theta\right)^2}\right)\right) d\eta \tag{24}$$

where *U* is a scaling velocity that depends on the streamwise coordinate (see below), whereas:

*<sup>w</sup>* <sup>=</sup> *<sup>u</sup> <sup>U</sup>* , *<sup>η</sup>* <sup>=</sup> *<sup>y</sup> δ* , (25)

In a similar way, one can integrate the energy Equation (11). This brings:

$$\frac{d(\mathcal{U}\mathfrak{E})}{d\mathfrak{x}}\int\_{0}^{1} w\mathfrak{H}d\mathfrak{\eta} = -\frac{\mathfrak{a}}{\mathfrak{d}} \left(\frac{\partial\mathfrak{G}}{\partial\mathfrak{\eta}}\right)\_{\mathfrak{\eta}=\mathfrak{0}} \tag{26}$$

To solve Equations (24) and (26), it is necessary to specify the velocity and temperature profiles in the boundary layer. Authors [30,31] proposed the following equations for these profiles:

$$w = \eta (1 - \eta)^2,\tag{27}$$

$$\theta = \left(1 - \eta\right)^{3},\tag{28}$$

The substitution of profiles (27) and (28) into the integral Equation (24) gives:

$$\frac{1}{105}\frac{d\left(\delta l I^2\right)}{d\chi} = -\mathbf{v}\frac{\mathcal{U}}{\delta} + \mathbf{g}\delta F\_{\prime} \tag{29}$$

where

$$\begin{pmatrix} F(\mathbf{\tilde{B}}\Delta T, \mathbf{W}\mathbf{a}\_d, \mathbf{W}\mathbf{a}\_b) \\ \frac{1}{54} \begin{pmatrix} 54 - \frac{3(-6\mathbf{W}\mathbf{a}\_b(1+\mathbf{\tilde{B}}\Delta T) + \mathbf{W}\mathbf{a}\_d(8+5\mathbf{\tilde{B}}\Delta T))}{(1+\mathbf{W}\mathbf{a}\_d - \mathbf{W}\mathbf{a}\_b)(1+\mathbf{\tilde{B}}\Delta T)^2} + \\ + \frac{(9+5\mathbf{W}\mathbf{a} - 6\mathbf{W}\mathbf{a}\_b)\left(6\arctan\left(\frac{1-2\frac{3}{\sqrt{5}\Delta T}}{\sqrt{5}}\right) + \sqrt{3}\ln\left(1-\frac{3\frac{3}{\sqrt{5}\Delta T}}{\left(1+\frac{3}{\sqrt{5}\Delta T}\right)^2}\right) - \pi\right) \end{pmatrix} \end{pmatrix} \tag{30}$$

In the limiting case of βΔ*T* → 0, Equation (29) transforms into:

$$\frac{1}{105}\frac{d\left(\\$lI^2\right)}{d\mathbf{x}} = -\mathbf{v}\frac{U}{\delta} + \frac{1}{4}g\delta\beta\Delta T\frac{1+3\mathbf{W}\mathbf{a}\_d - 2\mathbf{W}\mathbf{a}\_b}{1+\mathbf{W}\mathbf{a}\_d - \mathbf{W}\mathbf{a}\_b}.\tag{31}$$

For an ideal gas (Wa*<sup>a</sup>* = Wa*<sup>b</sup>* = 0), Equation (31) reduces to the respective equation obtained in the classical works [30,31]. Energy Equation (26) takes the following form:

$$\frac{1}{42}\frac{d(\delta \mathcal{U})}{dx} = 3\frac{\alpha}{\delta},\tag{32}$$

If Equation (28) for the temperature profile is replaced by another profile used in our study [28]:

$$\boldsymbol{\Theta} = \left(1 - \boldsymbol{\eta}\right)^{2},\tag{33}$$

then the form of the function *F* in Equation (30) will change as follows:

$$\begin{array}{l} F(\beta \Delta T, \text{Wa}\_d, \text{Wa}\_b) = \\ 1 - \frac{-4 \text{Wa}\_b(1 + \beta \Delta T) + \text{Wa}\_d(5 + 3 \beta \Delta T)}{8(1 + \text{Wa}\_d - \text{Wa}\_b)(1 + \beta \Delta T)^2} - \frac{(8 + 3 \text{Wa}\_d - 4 \text{Wa}\_b) \text{arctan} \left(\sqrt{\beta \Delta T}\right)}{8\sqrt{\beta \Delta T}(1 + \text{Wa}\_d - \text{Wa}\_b)} \end{array} \tag{34}$$

Then, in the limiting case of βΔ*T* → 0 used in the Boussinesq approximation [30,31], the integral Equation (29) looks as follows:

$$\frac{1}{105}\frac{d\left(\\$lI^2\right)}{d\mathbf{x}} = -\mathbf{v}\frac{U}{\delta} + \frac{1}{3}g\delta\beta\Delta T \frac{1+3\text{Wa}\_d-2\text{Wa}\_b}{1+\text{Wa}\_d-\text{Wa}\_b},\tag{35}$$

and the energy equation takes the form:

$$\frac{1}{30}\frac{d(\delta \mathcal{U})}{dx} = 2\frac{\alpha}{\delta},\tag{36}$$

Summing up, we can write the system of integral equations in a general form:

$$\mathrm{s\,s\,}\frac{d\left(\delta \mathcal{U}^2\right)}{d\mathbf{x}} = -\mathbf{v}\frac{\mathcal{U}}{\delta} + \mathbf{g}\delta \mathcal{F}\_{\prime} \tag{37}$$

$$t\frac{d(\delta \mathcal{U})}{d\mathfrak{x}} = z\frac{\mathfrak{a}}{\mathfrak{s}'}\tag{38}$$

where the function *F* is defined either by Equation (30) or Equation (34). Here *s* and *t* are numerical coefficients. Substituting Equation (27) into Equation (24), one can obtain *s* = 1/105; a substitution of Equation (27) into Equation (24) yields *t* = 1/30.
