Doklady of the Academy of Sciences of the USSR
1962, Volume 146, No. 5

MATHEMATICAL PHYSICS

L. I. RUBINSTEIN

ON THE ASYMPTOTICS OF THE SOLUTION OF A CONTACT AXISYMMETRIC THERMOCONVECTIVE PROBLEM FOR LARGE VALUES OF THE CONVECTIVE PARAMETER*

(Presented by Academician S. L. Sobolev, 21 V 1962)

Below we consider the asymptotic behavior, with respect to the parameter \(n \to \infty\), of the solution \(u(r,z,t)\) of the boundary-value problem:

\[ \frac{\partial^2 u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r} +\frac{\partial^2 u}{\partial z^2} = \frac{1}{a^2}\frac{\partial u}{\partial t}; \qquad 0<r<\infty;\quad 0<z<\infty;\quad t>0. \]

\[ \frac{\partial^2 u}{\partial r^2} +\frac{1-n}{r}\frac{\partial u}{\partial r} +\alpha\frac{\partial u}{\partial z} = \frac{\partial u}{\partial t}; \qquad 0<r<\infty;\quad z=0;\quad t>0; \tag{1} \]

\[ u\big|_{t=0}=0;\qquad u\big|_{r^2+z^2\to\infty}=0;\qquad \lim_{r\to0}\lim_{z\to0;\,t>0} u=1. \]

This problem was posed and solved by us earlier \((^{1,2})\) by means of the double transform

\[ \int_0^\infty w(p,s,z)J_0(sr)\,s\,ds = p\int_0^\infty e^{-pt}u(r,z,t)\,dt. \tag{2} \]

It was shown, in particular, that for \(a^2=1\) and integer \(n\)

\[ w(p,s,z)= n e^{-z\sqrt{s^2+p}}(\alpha+\sqrt{p})^n (\alpha+\sqrt{p+s^2})^{-(n+1)}(p+s^2)^{-1/2}. \tag{3} \]

The inversion of (2) proposed in \((^1)\), using the direct application of the Bromwich integral and Hankel’s integral theorem, gives a representation of the solution in a form unsuitable for numerical calculations for large values of the convective parameter \(n\), owing to the strongly oscillatory character of the integrand in the integral representation of the solution. Meanwhile, for the applications that were intended when formulating problem (1), precisely this case is of interest, which creates the need to obtain an asymptotic representation of the solution as \(n\to\infty\).

Let first \(a^2=1\) and \(n\) be an integer. Let \(\div\) denote the correspondence sign under the Laplace–Carson transform between the image and the original, let \(\sim\) denote the correspondence sign between Hankel transforms, and let \(\overset{\sim}{\div}\), \(\underset{\sim}{\div}\) denote the correspondence signs under the double Laplace–Carson transform in \(t\) and Hankel transform in \(r\), the upper bar corresponding to the inner transform and the lower to the outer one. Let \(w \div W\). Using the operational rules listed in Table B of the handbook \((^3)\) under Nos. 0.1; 0.3; 0.7; 0.25; 0.35; 0.50, expanding \((\alpha+\sqrt p)^n\) by the binomial formula, and using integration by parts and Leibniz’ formula, we readily find

\[ W= \frac{(-1)^n}{\Gamma(n)} \int_z^\infty e^{\alpha z} \frac{\partial^n}{\partial \xi^n} \left( e^{-\alpha \xi}\operatorname{erfc}\frac{\xi}{2\sqrt{t}} \right)d\xi \int_z^\xi J_0\!\left(s\sqrt{\xi^2-\eta^2}\right) e^{-\alpha(\eta-z)}(\eta-z)^n\,d\eta. \tag{4} \]

* The work was carried out on assignment from the association “Bashneft.”

From (2) it follows that \(w \simeq u\). Assuming that the change in the order of integration is legitimate, we find that \(W \sim u\). Hence, from (4) and the Hankel integral theorem, after simple transformations it follows that

\[ u=\frac{(-1)^n}{\Gamma(n)}\int_0^\infty e^{aR}\frac{\partial^n}{\partial R^n} \left(e^{-\alpha R}\operatorname{erfc}\frac{R}{2\sqrt{t}}\right)e^{-\alpha s}s^nR^{-1}ds . \tag{5} \]

Here \(R^2=(s+z)^2+r^2\).

Direct substitution of (5) into (1) shows that (5) is indeed a solution of the problem for \(a^2=1\) and integral \(n\).

Let \(C\) be the contour in the plane of the complex variable \(\zeta\), formed by the rays \(\arg \zeta=\pm\left(\frac{\pi}{4}-\varepsilon\right)\), \(\varepsilon>0\). Then (5) can be represented in the form

\[ u=\frac{(-1)^n n}{2\pi i}\int_0^\infty \frac{e^{a(R-s)}}{R}s^n ds \int_C e^{-\alpha\zeta}\operatorname{erfc}\frac{\zeta}{2\sqrt{t}}\, \frac{d\zeta}{(\zeta-R)^{n+1}}, \tag{6} \]

where the integral, considered as a double integral, converges absolutely. Put

\[ x=\frac{az}{n};\qquad y=\frac{ar}{n};\qquad \tau=\frac{a^2\alpha^2t}{n^2};\qquad \omega=\frac{a\zeta}{n};\qquad \lambda=\frac{as}{n}. \tag{7} \]

Remembering that, by assumption, \(a^2=1\), we obtain

\[ u=\frac{n}{2\pi i}\int_0^\infty \exp[n(1-\lambda+\ln\lambda)]\frac{d\lambda}{R_\lambda} \int_C \exp\{n[R_\lambda-\omega-1-\ln(R_\lambda-\omega)]\}\times \]

\[ {}\times \operatorname{erfc}\frac{\omega}{2\sqrt{\tau}}\, \frac{d\omega}{\omega-R_\lambda}. \tag{8} \]

Here \(R_\lambda^2=(\lambda+x)^2+y^2\). We shall use the saddle-point method twice, in the same way as was done by Szegő in deriving Perron’s formula for Laguerre polynomials \((^4)\). Let

\[ R^2=(1+x)^2+y^2;\qquad f_{s,m}=\frac{\partial^s}{\partial x^s} \left(\frac{1}{R}\frac{\partial^m}{\partial R^m} \operatorname{erfc}\frac{R-1}{2\sqrt{\tau}}\right); \]

\[ \alpha_{k,s}=\int_{-\infty}^{+\infty} e^{-\rho^2/2}(\rho i)^s A_{2k-s}(\rho)\,d\rho;\qquad \beta_{k,s}=\int_{-\infty}^{+\infty} e^{-\rho^2/2}\rho^s B_{2k-s}(\rho)\,d\rho; \]

\[ \exp\{-n[\rho i n^{-1/2}-\ln(1-\rho i n^{-1/2})]\} \simeq e^{-\rho^2/2}\sum_{k=0}^{\infty}A_k(\rho)n^{-k/2}, \]

\[ \exp\{n[-\rho n^{-1/2}+\ln(1+\rho n^{-1/2})]\} \simeq e^{-\rho^2/2}\sum_{k=0}^{\infty}B_k(\rho)n^{-k/2}. \]

Here \(\simeq\) denotes the sign of asymptotic correspondence. It is not difficult to verify that

\[ u\simeq \sum_{j=0}^{\infty}u_j n^{-j}, \tag{9} \]

where

\[ u_j=\frac{1}{2\pi}\sum_{p=0}^{j}\sum_{k=0}^{2p}\sum_{m=0}^{k}\sum_{s=0}^{2(j-p)} \frac{1}{m!s!}\,\alpha_{p,k}f_{s,m}\beta_{j-p,s}. \tag{10} \]

The calculation gives, in particular,

\[ u_0=\frac{1}{R}\operatorname{erfc}\frac{R-1}{2\sqrt{\tau}}, \qquad u_1=\left[\frac{\partial}{\partial x} -\frac{1}{2}\left(\frac{\partial^2}{\partial y^2} +\frac{1}{y}\frac{\partial}{\partial y}\right)\right]u_0 . \tag{11} \]

It is easy to see that \(u_0\) and \(u_1\) are solutions of the boundary-value problems

\[ \left(\frac{\partial^2}{\partial x^2} +\frac{\partial^2}{\partial y^2} +\frac{1}{y}\frac{\partial}{\partial y} -\frac{\partial}{\partial \tau}\right)u_i=0; \qquad x>0;\quad y>0;\quad \tau>0; \]

\[ \left(\frac{\partial}{\partial x} -\frac{1}{y}\frac{\partial}{\partial y}\right)u_i+\psi_i=0; \qquad x=0;\quad y>0;\quad \tau>0; \tag{12} \]

\[ u_i\big|_{\tau=0}=0;\qquad u_i\big|_{x^2+y^2\to\infty}=0;\qquad u_i\big|_{x=y=0;\,\tau>0}=\Psi_i;\qquad i=0,1,2,\ldots; \tag{12*} \]

\[ \psi_0\equiv 0;\qquad \Psi_0=1;\qquad \psi_i=\left(\frac{\partial^2}{\partial y^2} +\frac{1}{y}\frac{\partial}{\partial y} -a^2\frac{\partial}{\partial \tau}\right)u_{i-1}; \qquad \Psi_i=0,\quad i=1,2,\ldots, \]

if in them one sets \(a^2=1\).

We now abandon the assumption “\(a^2=1\) and \(n\) is an integer.” Let us introduce the variables (7) directly into (1) and seek \(u\) in the form of the series

\[ u=\sum_{j=0}^{\infty} u_j n^{-j}. \tag{13} \]

Assuming that the series (13) may be differentiated twice term by term, we find that the \(u_j\) must be solutions of the boundary-value problems (12). The convergence of the series (13) and the legitimacy of its twice term-by-term differentiation can hardly be proved. However, the coincidence of the first two terms of this series (for \(a^2=1\)) with the first two terms of the asymptotic series (9), obtained by the saddle-point method, permits one to regard the series (13) only as asymptotic. The terms of the series (13) can be found recursively. Namely, put \(\mathfrak{w}_m \doteq u_m;\ \mathfrak{w}_m \doteq W_m\), so that formally, under the assumption that the order of integration may be changed, \(W_m\sim u_m\). Repeating the arguments carried out in (1), we show that the \(\mathfrak{w}_m\) are determined recursively by the equalities

\[ \mathfrak{w}_0=(p+s^2)^{-1/2} \exp\left[\sqrt{p}-(1+x)\sqrt{p+s^2}\right]; \tag{14} \]

\[ \mathfrak{w}_{m+1} =\mathfrak{w}_0(p,s) \int_0^s \frac{\lambda^2+a^2p}{\lambda^2+p} \frac{\mathfrak{w}_m(p,\lambda)}{\mathfrak{w}_0(p,\lambda)} \,\lambda\,d\lambda -\frac{s^2+a^2p}{\sqrt{p+s^2}}\,\mathfrak{w}_m(p,s). \tag{15} \]

Inverting the transformation \(\mathfrak{w}_m \doteq u_m\) will give \(u_m\). Let us carry out the computation of \(u_0\) and \(u_1\) to the end. Since \(u_0\) does not depend on \(a^2\), it is evident that it coincides with the previously found \(u_0\) (see (11)). Substituting (14) into (15) and integrating, we find

\[ \mathfrak{w}_1 =\sum_{m=1}^{4}\mathfrak{w}_{1m} =\frac{s^2}{2}\mathfrak{w}_0 -\sqrt{s^2+p}\,\mathfrak{w}_0 -\frac{a^2-1}{\sqrt{s^2+p}}\,p\mathfrak{w}_0 +\frac{a^2-1}{2}\,p\ln\left(1+\frac{s^2}{p}\right)\mathfrak{w}_0. \tag{16} \]

Let \(\mathfrak{w}_{1i}\doteq W_{1i}\sim u_{1i}\), so that \(\mathfrak{w}_{1i}\doteq u_{1i}\), or formally \(\mathfrak{w}_1\doteq u_{1i}\); hence,

\[ u_1=\sum_{i=1}^{4}u_{1i}. \]

As is known (3),

\[ p\ln\left(1+\frac{s^2}{p}\right) \doteq \tau^{-1}\left[1-\exp(-s^2\tau)\right]. \]

Hence, and from the convolution theorem, it follows that

\[ W_{14}=\frac{a^2-1}{2}\frac{\partial}{\partial \tau} \int_0^\tau \frac{1-\exp(-s^2\lambda)}{\lambda}\, W_0(x,s,\tau-\lambda)\,d\lambda. \]

Therefore, formally,

\[ u_{14} =\frac{a^2-1}{2}\frac{\partial}{\partial \tau} \int_0^\tau d\lambda \int_0^\infty \frac{1-\exp(-s^2\lambda)}{\lambda}\, W_0(x,s,\tau-\lambda)J_0(sy)\,s\,ds. \tag{17} \]

But \(W_0\) is the Hankel transform of order zero of \(u_0\). Using this and Weber’s second exponential integral (5), we obtain

\[ u_{14}=\frac{a^2-1}{2}\frac{\partial}{\partial \tau} \int_0^\tau \frac{d\lambda}{\lambda}\left[ u_0(x,y,\tau-\lambda)- \right. \]

\[ \left. -\int_0^\infty u_0(x,\eta,\tau-\lambda) \frac{\exp[-(y^2+\eta^2)/4\lambda]}{2\lambda} I_0\!\left(\frac{y\eta}{2\lambda}\right)\eta\,d\eta \right]. \tag{18} \]

It can further be shown that \(W_0(x,s,0)=0\). Consequently, multiplication of \(w_0\) by \(p\) corresponds to differentiation of \(u_0\) with respect to \(\tau\). From (14) we see further that multiplication of \(W_0\) by \(-\sqrt{p+s^2}\) corresponds to differentiation of \(u_0\) with respect to \(x\), while division by \(\sqrt{p+s^2}\) corresponds to integration of \(u_0\) with respect to \(x\) from \(x\) to infinity. Finally, multiplication of \(w_0\) by \(s^2\) means applying the operator
\[ -\left(\frac{\partial^2}{\partial y^2}+\frac{1}{y}\frac{\partial}{\partial y}\right) \]
to \(u_0\). Comparing these remarks with (16) and using (18), we obtain

\[ u_1=\left[ \frac{\partial x}{\partial x} -\frac12\left(\frac{\partial^2}{\partial y^2} +\frac1y\frac{\partial}{\partial y}\right) \right]u_0 -(a^2-1)\frac{\partial}{\partial\tau}\int_x^\infty u_0(\xi,y,\tau)\,d\xi+ \]

\[ +\frac{a^2-1}{2}\frac{\partial}{\partial\tau} \int_0^\tau \frac{d\lambda}{\lambda}\left[ u_0(x,y,\tau-\lambda)- \right. \]

\[ \left. -\int_0^\infty u_0(x,\eta,\tau-\lambda) \frac{\exp[-(y^2+\eta^2)/4\lambda]}{2\lambda} I_0\!\left(\frac{y\eta}{2\lambda}\right)\eta\,d\eta \right], \tag{19} \]

which for \(a^2=1\) coincides with (11).

A more delicate analysis shows that all the operations that led to the representation of \(u_1\) in the form (19) are legitimate, so that (19) is indeed a solution of problem (12). A direct verification of this is difficult. We do not have estimates of the remainder terms of the asymptotic expansions obtained. However, in the case \(a^2=1\) and integer \(n\leq 20\), the computation of \(u\) can be carried out by the exact formula (5). Comparison of the computational results with the computation by the asymptotic formula \(u \cong u_0+\frac1n u_1\) shows that, in the range of values of the arguments for which the asymptotics gives a value of \(u\) not less than 0.2, the error from replacing the exact solution by the asymptotic one does not exceed 0.005 already for \(n=10\), if (for \(n=10\)) \(\tau\geq 0.1\).* Deterioration of the result occurs in the region of large \(x^2+y^2\) and, at the same time, small \(\tau\), where the solution becomes close to zero. In this region, to refine the result one has to take into account the next term of the asymptotic expansion, which we do not write down for lack of space. From the applied point of view, taking two terms of the asymptotic expansion into account proves to be quite sufficient.

All numerical computations connected with the present work, and their programming on the BESM-2, were carried out by Kh. Geiman and N. A. Avdonin, to whom we are pleased to express our gratitude.

Latvian State University
named after P. Stučka

Received
3 V 1962

REFERENCES

1. L. I. Rubinshtein, DAN, 135, No. 4 (1960).
2. L. I. Rubinshtein, Uch. zap. Kazan. gos. univ., 121, book 5, 129 (1961).
3. V. A. Ditkin, P. I. Kuznetsov, Handbook on Operational Calculus, M.—L., 1951.
4. G. Szegő, Orthogonal Polynomials, N. Y., 1939.
5. G. N. Watson, Theory of Bessel Functions, part 1, IL, 1949.

* As \(n\) increases, the lower bound of the admissible \(\tau\) decreases.