一类非齐次变型Bessel方程边值问题的构造解法及应用

The Constructive Solution and Application of BVP for a Class of Non-Homogeneous Bessel Equations

摘要: 本文给出了一类二阶线性常微分方程边值问题的构造解法，运用此方法进一步求得了在相同边界条件下的非齐次零阶变型Bessel方程边值问题的解，并将该理论应用于具有源汇影响的变流率平面径向渗流力学问题中，结合Laplace变换法求得了该渗流力学模型在三种外边界条件下的连分式解式。研究表明，该方法可以方便、快捷地构造出此类边界条件下的统一解式，对同类边界条件下的扩散(热传导、渗流)方程求解具有重要意义。

Abstract: This paper gives a constructive solution of a class of second order linear ordinary differential equation boundary value problem, by using this method, a non-homogeneous zero order modified Bessel equation BVP in the same boundary conditions is solved. Then the theory is also applied to the plane radial seepage mechanics with variable flow rate and source sink effect, combined with Laplace transform method, the continued fraction solutions of the seepage mechanics model under three kinds of external boundary conditions are obtained. The results show that this method can easily and quickly construct the unified solution under the same boundary conditions, which is of great significance for solving the diffusion (heat conduction, seepage) equation under the same boundary conditions.

1. 引言

在解决数学物理问题如热的传导、浓度扩散以及工程中的流体渗流力学等问题中常常会遇到扩散方程，它们在数学本质上都是二阶抛物型微分方程，因此研究二阶微分方程的求解方法具有重要意义。

对于二阶线性齐次偏微分方程的解法已经非常成熟，常见的有分离变量法，是将偏微分方程分离为多个常微分方程分别求解，再根据线性方程叠加原理求得形式解；积分变换法是通过构造正交函数系求得积分变换对，将偏微分方程化为常微分方程，将常微分方程化为代数方程来求解 [1]。具有汇源影响的方程中通常会出现非齐次方程或非齐次边界条件，求解此类方程比较常用的是Green函数法，但需要借助分离变量法求得对应的齐次问题的解，从而构造相应的Green函数，然后将此Green函数与方程的非齐次项或非齐次边界作卷积运算求得解 [2]，且运算过程较为繁琐。对非齐次方程的求解更为一般的有常数变易法、微分算子法、Laplace变换法、待定系数法等，虽解法各有优劣，但在很大程度上都需要求解常微分方程。一般的二阶常系数线性微分方程的求解归结为求取方程的特征值和特征函数 [3]。对特殊的变系数二阶常微分方程，如Bessel方程 [4]、Euler方程、Legendre方程等，其通解是由各方程对应的特殊函数所构成 [5]。在非稳态扩散(热传导、渗流)方程中，Bessel方程最为常见，求解其所构成的初边值问题是研究此类问题的关键。李顺初 [6] [7] [8] [9] 在求解二阶齐次线性微分方程(组)初边值问题的研究中发现该类问题的解式都具有相似的连分式结构，并在研究油气藏试井解释模型中 [10] [11]，同样得到具有相似性的连分式解。进一步地，该研究团队又提出了求解二阶常微分方程边值问题的相似构造法 [12] [13] [14]，并获得了丰硕的成果。

本文将进一步把该理论应用于求解二阶非齐次零阶变型Bessel方程的边值问题，探索该非齐次方程边值问题的构造解法，并将其应用于具有源汇影响的变流率渗流问题中，结合Laplace变换求得了模型的半解析解，应用数值反演技术可以满足工程需要。该研究表明此方法快捷有效、应用范围广，是对微分方程求解方法的拓展，对此类扩散方程初边值问题求解具有重要意义。

2. 二阶线性微分方程的构造解

引理1 如果二阶线性微分方程边值问题

$y^{″}+p\left(x\right)y^{\prime }+q\left(x\right)y=0$ (1)

${\left[Ey+\left(1+EF\right)y^{\prime }\right]}_{x=a}=D$ (2)

${\left[Gy+H{y}^{\prime }\right]}_{x=b}=0$ (3)

(其中 $D,E,F,G,H,a,b$ 均为已知的实常数，且 $D\ne 0,G^2+H^2\ne 0,a<b$ ；已知函数 $p\left(x\right)\in C^1\left[a,b\right],q\left(x\right)\in C^2\left[a,b\right]$ )有唯一解，则其解可以表示为如下的连分式形式 [12]

$y\left(x\right)=D\cdot \frac{1}{E+\frac{1}{F+\Phi \left(a\right)}}\cdot \frac{1}{F+\Phi \left(a\right)}\cdot \Phi \left(x\right)$ (4)

其中

$\Phi \left(x\right)=\frac{G{\phi }_{0,0}\left(x,b\right)+H{\phi }_{0,1}\left(x,b\right)}{G{\phi }_{1,0}\left(a,b\right)+H{\phi }_{1,1}\left(a,b\right)}$ (5)

称为解的相似核函数。

假设 $y_1\left(x\right)$， $y_2\left(x\right)$ 是(1)式方程的两个线性无关的解，构造以下二元函数

${\phi }_{0,0}\left(\xi ,x\right)=y_1\left(\xi \right)y_2\left(x\right)-y_2\left(\xi \right)y_1\left(x\right)$ (6)

并记

${\phi }_{1,0}\left(\xi ,x\right)=\frac{\partial }{\partial \xi }{\phi }_{0,0}\left(\xi ,x\right)={y^{\prime }}_1\left(\xi \right)y_2\left(x\right)-{y^{\prime }}_2\left(\xi \right)y_1\left(x\right)$ (7)

${\phi }_{0,1}\left(\xi ,x\right)=\frac{\partial }{\partial x}{\phi }_{0,0}\left(\xi ,x\right)=y_1\left(\xi \right){y^{\prime }}_2\left(x\right)-y_2\left(\xi \right){y^{\prime }}_1\left(x\right)$ (8)

${\phi }_{1,1}\left(\xi ,x\right)=\frac{{\partial }^2}{\partial \xi \partial x}{\phi }_{0,0}\left(\xi ,x\right)=\frac{{\partial }^2}{\partial x\partial \xi }{\phi }_{0,0}\left(x,\xi \right)={y^{\prime }}_1\left(\xi \right){y^{\prime }}_2\left(x\right)-{y^{\prime }}_2\left(\xi \right){y^{\prime }}_1\left(x\right)$ (9)

称为相似核函数的生成函数。

引理2 在(3)式外边界条件 ${\left[Gy+H{y}^{\prime }\right]}_{x=b}=0$ 中，若 $G\ne 0,H=0$，则外边界条件

变为第一类边界条件 $y\left(b\right)=0$，相应的核函数为

$\Phi \left(x\right)=\frac{{\phi }_{0,0}\left(x,b\right)}{{\phi }_{1,0}\left(a,b\right)}$ (10)

在(3)式外边界条件 ${\left[Gy+H{y}^{\prime }\right]}_{x=b}=0⤢$ 中，若 $G=0,H\ne 0$，则外边界条件变为第二类边界条件 $y^{\prime }\left(b\right)=0$，相应的核函数为

$\Phi \left(x\right)=\frac{{\phi }_{0,1}\left(x,b\right)}{{\phi }_{1,1}\left(a,b\right)}$ (11)

在边值问题(1)~(3)式的解(4)式中连分式结构中的第一个连分式有如下性质

${\left[y\left(x\right)+F{y}^{\prime }\left(x\right)\right]}_{x=a}=\frac{D}{E+\frac{1}{F+\Phi \left(a\right)}}$ (12)

此式反映了解在内边界处的本质性的特征，在渗流力学理论中起着十分重要的作用。

定理1 如果Bessel函数 $I_0\left(\sqrt{m}x\right),K_0\left(\sqrt{m}x\right)$ 是零阶齐次变型Bessel方程

$x^2{y}^{″}+x{y}^{\prime }-m{x}^{2}y=0$ (13)

的两个线性无关的解，则其通解为

$y\left(x\right)=C_1{I}_0\left(\sqrt{m}x\right)+C_2{K}_0\left(\sqrt{m}x\right)$ (14)

则方程(13)式与边界条件(2) (3)式所构成的边值问题仍然具有形如(4)式的解式 [6]，

其中的相似核函数形式也与(5)式完全相同，区别只在于二元生成函数

$\begin{array}{c}{\phi }^{*}{}_{0,0}\left(\xi ,x\right)=y_1\left(\xi \right)y_2\left(x\right)-y_2\left(\xi \right)y_1\left(x\right)\\ =I_0\left(\sqrt{m}\xi \right)K_0\left(\sqrt{m}x\right)-K_0\left(\sqrt{m}\xi \right)I_0\left(\sqrt{m}x\right)=-{\psi }_{0,0}^{*}\left(\xi ,x,\sqrt{m}\right)\end{array}$ (15)

根据Bessel函数的微分性质有

${\phi }^{*}{}_{1,0}\left(\xi ,x\right)={I^{\prime }}_0\left(\sqrt{m}\xi \right)K_0\left(\sqrt{m}x\right)-{K^{\prime }}_0\left(\sqrt{m}\xi \right)I_0\left(\sqrt{m}x\right)=\sqrt{m}{\psi }_{1,0}^{*}\left(\xi ,x,\sqrt{m}\right)$ (16)

${\phi }^{*}{}_{0,1}\left(\xi ,x\right)=I_0\left(\sqrt{m}\xi \right){K^{\prime }}_0\left(\sqrt{m}x\right)-K_0\left(\sqrt{m}\xi \right){I^{\prime }}_0\left(\sqrt{m}x\right)=-\sqrt{m}{\psi }_{0,1}^{*}\left(\xi ,x,\sqrt{m}\right)$ (17)

${\phi }^{*}{}_{1,1}\left(\xi ,x\right)={I^{\prime }}_0\left(\sqrt{m}\xi \right){K^{\prime }}_0\left(\sqrt{m}x\right)-{K^{\prime }}_0\left(\sqrt{m}\xi \right){I^{\prime }}_0\left(\sqrt{m}x\right)=m{\psi }_{1,1}^{*}\left(\xi ,x,\sqrt{m}\right)$ (18)

这里

${\psi }_{m,n}^{*}\left(x,y,t\right)=K_m\left(xt\right)I_n\left(yt\right)+{\left(-1\right)}^{m-n+1}{I}_m\left(xt\right)K_n\left(yt\right)$ (19)

其中 $K_v\left(x\right),I_v\left(x\right)$ 是 $v$ 阶变型Bessel函数。

利用上式得相似核函数

${\Phi }^{\text{*}}\left(x\right)=\frac{G{\phi }_{0,0}^{\text{*}}\left(x,b\right)+H{\phi }_{0,1}^{\text{*}}\left(x,b\right)}{G{\phi }_{1,0}^{\text{*}}\left(a,b\right)+H{\phi }_{1,1}^{\text{*}}\left(a,b\right)}$

其中上标“*”表示由变型Bessel函数生成的函数。

定理2 对于非齐次零阶变型Bessel方程

$x^2{y}^{″}+x{y}^{\prime }-m{x}^{2}y=-\frac{1}{n}$ (20)

采用常数变易法不难求得特解为(证明见附录)

$y^{*}\left(x\right)=\frac{1}{mn}$ (21)

则非齐次方程(20)式的通解为

$y\left(x\right)=C_1{I}_0\left(\sqrt{m}x\right)+C_2{K}_0\left(\sqrt{m}x\right)+\frac{1}{mn}$ (22)

则由非齐次零阶变型Bessel方程(20)式与边界条件(2) (3)式联立构成的边值问题的解为对应齐次方程边值问题的解与非齐次方程特解之和

$y\left(x\right)=D\cdot \frac{1}{E+\frac{1}{F+{\Phi }^{*}\left(a\right)}}\cdot \frac{1}{F+{\Phi }^{*}\left(a\right)}\cdot {\Phi }^{*}\left(x\right)+\frac{1}{mn}$ (23)

3. 非稳态变流率渗流模型的构造解

将上述边值问题解的构造理论应用于带有源汇影响的非稳态变流率平面径向渗流模型可以直接构造得出模型的解。

无因次渗流方程

$\frac{{\text{d}}^2{p}_D}{\text{d}{r}_D^2}+\frac{1}{r_D}\frac{\text{d}{p}_D}{\text{d}{r}_D}+\frac{1}{\lambda }{q}_D\left(t_D\right)=\frac{\text{d}{p}_D}{\text{d}{t}_D}\left(r_D>0\right)$ (24)

考虑井筒储集效应和表皮影响的内边界条件

$P_{wD}\left(t_D\right)={\left[P_D-S{r}_D\frac{\text{d}{p}_D}{\text{d}{r}_D}\right]}_{r_{D=1}}$ (25)

$C_D\frac{\text{d}{p}_{wD}}{\text{d}{t}_D}-\frac{\partial p_D}{\partial r_D}|{}_{r_{D=1}}=f\left(t_D\right)$ (26)

初始条件 $P\left(r_D,0\right)=0$ (27)

无因次外边界条件

$P_D\left(\infty ,t_D\right)=0\left(无穷大外边界\right)$ (28)

$\frac{\partial p_D}{\partial r_D}|{}_{r_{D=R_D}}=0\left(不渗透封闭外边界\right)$ (29)

$P\left(R_D,t_D\right)=0\left(定压外边界\right)$ (30)

其中： $C_D$ ——无因次井筒储集系数；

$S$ ——表皮因子；

$q_D\left(t_D\right)$ ——无因次井源汇强度；

$f\left(t_D\right)$ ——无因次变流率；

无因次变量定义如下

$p_D=\frac{2\pi kh\left(p_0-p\right)}{Bq\mu },t_D=\frac{kt}{\mu \phi C_t{r}_w^2},r_D=\frac{r}{r_w}，R_D=\frac{R}{R_w}$

将渗流问题(24)~(30)式作关于时间 $t_D$ 的Laplace变换得以下二阶常微分非齐次方程边值问题

$\frac{{\text{d}}^2{\bar{p}}_D}{\text{d}{r}_D^2}+\frac{1}{r}\frac{\text{d}{\bar{p}}_D}{\text{d}{r}_D}+\frac{1}{\lambda }{\bar{q}}_D\left(z\right)=z{\bar{p}}_D\left(r_D>0\right)$ (31)

${\bar{P}}_{wD}\left(z\right)={\left[{\bar{p}}_D-S{r}_D\frac{\text{d}{\bar{p}}_D}{\text{d}{r}_D}\right]}_{r_{D=1}}$ (32)

$C_{D}z{\bar{p}}_{wD}-\frac{\text{d}{\bar{p}}_D}{\text{d}{r}_D}|{}_{r_{D=1}}=\bar{f}\left(z\right)$ (33)

$\begin{array}{l}\bar{p}\left(\infty ,z\right)=0\\ \frac{\text{d}{\bar{p}}_D}{\text{d}{r}_D}|{}_{r_{D=R_D}}=0\\ {\bar{p}}_D\left(R_D,z\right)=0\end{array}$ (34)

由线性方程叠加原理不难证明，方程(31)的解

${\bar{p}}_D\left(r_D,z\right)=q\left(z\right){\bar{p}}_{1D}\left(r_D,z\right)$ (35)

其中 ${\bar{p}}_{1D}\left(r_D,z\right)$ 是方程

$\frac{{\text{d}}^2{\bar{p}}_{1D}}{\text{d}{r}_D^2}+\frac{1}{r}\frac{\text{d}{\bar{p}}_{1D}}{\text{d}{r}_D}+\frac{1}{\lambda }=z{\bar{p}}_{1D}$ (36)

的解，因此渗流问题(31)~(34)式的求解可以转化为以下边值问题的求解，同时合并(32) (33)式有

$\frac{{\text{d}}^2{\bar{p}}_{1D}}{\text{d}{r}_D^2}+\frac{1}{r}\frac{\text{d}{\bar{p}}_{1D}}{\text{d}{r}_D}+\frac{1}{\lambda }=z{\bar{p}}_{1D}$ (37)

$C_{D}z{\bar{p}}_{1D}-\left(1+C_{D}zS{r}_D\right)\frac{\text{d}{\bar{p}}_{1D}}{\text{d}{r}_D}|{}_{r_{D=1}}=\frac{\bar{f}\left(z\right)}{{\bar{q}}_D\left(z\right)}$ (38)

$\begin{array}{l}{\bar{p}}_{1D}\left(\infty ,z\right)=0\\ \frac{\text{d}{\bar{p}}_{1D}}{\text{d}{r}_D}|{}_{r_{D=R_D}}=0\\ {\bar{p}}_{1D}\left(R_D,z\right)=0\end{array}$ (39)

只要令，，，，

则边值问题(37)~(39)式与(20)、(2)、(3)式所构成的边值问题完全等价，根据(23)式可得其解为

(40)

根据(35)式，可得渗流问题(31)~(34)式的Laplace空间解

(41)

其中的相似核函数为

(42)

由引理2中(10)式知当时，可得对应渗流模型定压外边界条件下的相似核函数

(43)

由(11)式知，当时，可得对应渗流模型不渗透封闭外边界条件下的相似核函数

(44)

当渗流模型外边界无穷大时相似核函数变为

(45)

只要将渗流模型不同外边界条件下对应的相似核函数(43)~(45)式代入(41)式即可得到该油藏模型在三种外边界条件下的统一连分式解。

由于对(41)式作Laplace逆变换较为困难，运用Stehfest数值反演技术结合计算机软件，可以获得实空间压力动态典型曲线，可以满足工程需要。

4. 结论

1) 本文提出的一类二阶线性常微分方程边值问题在三类外边界条件下具有统一的连分式形式解，并且内边界条件系数决定解的结构，外边界条件决定解的相似核函数；

2) 根据线性微分方程解的结构及叠加原理，求得了非齐次零阶变型Bessel方程边值问题的解式，其解可以根据相应的齐次变型Bessel方程边值问题的构造求得；

3) 多孔介质内流体的非稳态渗流方程、热传导方程以及扩散方程在本质上都可以归结为本文提出的非齐次零阶变型Bessel方程，他们在数学本质上是一致的。该研究是对二阶非齐次扩散方程求解方法的一种扩展，对求解此类方程具有重要意义。

基金项目

四川省教育厅项目(No. 15ZB0450、18ZB0439)，教育部产学研合作办学项目(No. 201702028006)，西华大学人才引进项目(No. Z201076)，宜宾职业技术学院基金项目(No. YBZYSC20BK05)。

附录

关于定理2中非齐次零阶变型Bessel方程特解的求得证明如下

定理2中非齐次零阶变型Bessel方程

(A1)

对应的齐次方程为

(A2)

由定理(1)知方程通解为

(A3)

下面采用常数变易法求特解，假设非齐次方程(A1)式的通解为

(A4)

则由以下方程组确定

(A5)

根据Bessel函数的性质

可得Wronski行列式

根据线性方程组的Cramer法则可得

(A6)

再根据Bessel函数的积分性质

得

若只求特解，则有

(A7)

(A8)

代入(A4)式可得方程(A1)的特解

(A9)