本发明涉及煤矿区矿井的动力学技术领域，具体而言，涉及一种不稳定承压－无压井流模型及承压含水层参数反演方法。
背景技术：
不稳定承压－无压抽水井流是煤矿区矿井疏干过程中的经典地下水动力学问题。近60年来，国内外众多学者就承压－无压井流的问题进行了针对性的研究，提出了多种基于泰斯井流假设或者裘布依假设的承压－无压完整井流的数值模型和解析模型。最先开始利用数值模型研究承压－无压井流渗流场特征的是一些国外学者。rushton和wedderburn(1971)首次建立了基于二维渗流基本微分方程的稳定承压－无压完整井流有限差分数值模型。该模型忽略了在承压流转换为无压流的过程中导水系数(transmissivity)的变化，在无压流区采用单位产水量(specificyield)取代承压区的储水系数(storativity)，提出利用基于电阻－电容网络模拟原理的有限差分法分析承压流转变为无压流的问题。elango和swaminathan(1980)提出基于四边混合曲线等参数单元的稳定承压－无压井流有限元数值模型。wang和zhan(2009)则提出一种更为完善的不稳定承压－无压完整井流的半解析数值模型。该模型基于二维渗流基本微分方程和潜水流动的布西涅斯克微分方程(boussinesqequation)，考虑了在承压流到无压流的转换过程中导水系数和给水度(storagecoefficient)的变化，采用四阶龙格－库塔法(runge–kuttamethod)求解非线性无压流问题。
解析模型相较于数值模型，具有易操作的优点。1961年，我国学者陈崇希等首次提出了单孔抽水形成的稳定承压－无压完整井流解析模型。moench和prickett(1972)认为当无压流区水位降深远远小于承压含水层本身的厚度时，无压流区的导水系数可视为唯一常量，并提出采用热传导中导管周围冻结问题的解来分析二维不稳定承压－无压完整井流的渗流场特征，即mp模型。该模型可用于研究渗流场中承压－无压转换界面随时间变化的时间动态分布特征，因此在解决不稳定承压－无压井流的工程实际问题中得到了广泛的应用(duffield，2007)。chenetal(2006)基于吉林斯基势函数(girinskii’spotential)，在考虑了河流补给边界条件作用下，提出了井群抽水形成的稳定承压－无压完整井流的解析模型。陈崇希(1961，1996)以及hu和chen(2008)则认为在承压－无压含水层中，导水系数和给水度分时段可被视为空间上的平均值，而不是在转换界面发生突变，并在此假设基础上提出了基于泰斯井流假设和吉林斯基势函数(girinskii’spotential)的二维不稳定承压－无压完整井流解析模型，即chen模型。
上述承压－无压井流模型为解决工程问题提供了一定的理论基础和模型方法。然而，在人类开采地下水的活动加剧和对非饱和流作用的认识不断深入的情况下，传统的不稳定承压－无压井流模型在实际应用中面临着挑战。传统的二维不稳定承压－无压完整井流解析模型均未能充分考虑承压流转换到无压流过程中导水系数、释水系数以及扩散系数(diffusivity)等水力参数变化对渗流场的影响。例如，在mp模型(moenchandprickett，1972)中将导水系数视为常量；chen模型(huandchen，2008)则假设不稳定承压－无压井流渗流场的扩散系数为分时段空间上的平均值。值得注意的是wang和zhan(2009)提出半解析数值模型虽然考虑了不稳定承压－无压完整井流渗流场中导水系数、给水度以及扩散系数等水力参数变化，但是该模型理论推导过程过于复杂，限制了其在实际工程中的推广应用。
技术实现要素：
本发明的目的在于提供一种不稳定承压－无压井流模型的计算方法，方法简单，有利于推广。
本发明的另一目的在于提供一种承压含水层参数反演的计算方法，方法简单，有利于推广。
本发明是采用以下技术方案实现的：
一种不稳定承压－无压井流模型的计算方法，包括如下步骤：
(1)、根据承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水速率、抽水时间、水平导水率、抽水井流函数和抽水井与观测井之间的距离计算出承压－无压转换界与抽水井之间的距离和无压区平均水头；
(2)、根据承压－无压转换界与抽水井之间的距离和无压区平均水头计算出承压－无压井流水位降深。
进一步地，本发明较佳的实施例中，上述计算承压－无压转换界与抽水井之间的距离和无压区平均水头的具体方法为：
通过求解第一方程组得到，第一方程组为：
其中，q为抽水速率，m3/s；t为抽水时间；sy为单位产水量；r为承压－无压转换界距离抽水井距离，m；r为抽水井和观测井之间的距离，m；b为含水层厚度，m；kr为水平导水率m/s；h为初始水头；s为给水度；t为导水系数，m2/s；w()为抽水井流函数；
进一步地，本发明较佳的实施例中，上述计算承压－无压井流水位降深的具体方法为：判定承压－无压转换界与抽水井之间的距离和抽水井与观测井之间的距离的大小关系，计算出承压－无压井流水位降深。
进一步地，本发明较佳的实施例中，上述当承压－无压转换界与抽水井之间的距离小于抽水井与观测井之间的距离时，通过求解第一方程得到承压－无压井流水位降深，第一方程为：
其中：h2(r，t)为承压态井流水头，m，即承压－无压井流水位降深。
进一步地，本发明较佳的实施例中，上述当承压－无压转换界与抽水井之间的距离不小于抽水井与观测井之间的距离时，通过求解第二方程得到承压－无压井流水位降深，第二方程为：
其中：h1(r，t)为无压态井流水头，m，即承压－无压井流水位降深。
一种承压含水层的参数反演的计算方法，包括如下步骤：
(1)、根据承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水时间、抽水井流函数、观测井水位和抽水井与观测井之间的距离计算出承压－无压转换界与抽水井之间的距离；
(2)、根据承压－无压转换界与抽水井之间的距离计算出无压区扩散率和无压区单位产水量。
进一步地，本发明较佳的实施例中，上述计算承压－无压转换界与抽水井之间的距离的具体方法为：判定观测井水位与含水层厚度的大小关系，计算出承压－无压转换界与抽水井之间的距离。
进一步地，本发明较佳的实施例中，上述当观测井水位不小于含水层厚度时，通过求解第三方程得到承压－无压转换界与抽水井之间的距离，第三方程为：
当观测井水位小于含水层厚度时，通过求解第四方程得到承压－无压转换界与抽水井之间的距离，第四方程为：
其中，h′(r1，t)为观测井水位，m；q为抽水速率，m3/s；t为抽水时间；r1为抽水井和观测井之间的距离，m；b为含水层厚度，m；kr为水平导水率m/s；h为初始水头；s为给水度；t为导水系数，m2/s；w()为抽水井流函数；r为承压－无压转换界距离抽水井距离，m。
进一步地，本发明较佳的实施例中，上述计算无压区单位产水量的具体方法是通过求解第二方程组得到，第二方程组为：
其中，
进一步地，本发明较佳的实施例中，上述计算无压区扩散系数的具体方法是通过求解第五方程得到，第五方程为：
其中，d为无压区扩散系数。
本发明的较佳实施例提供的不稳定承压－无压井流模型的计算方法的有益效果是：充分考虑了影响承压－无压转换界与抽水井之间的距离和无压区平均水头的各参数承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水速率、抽水时间、水平导水率、抽水井流函数和抽水井与观测井之间的距离的变化，使得到的结果更加准确，从而使计算出的承压－无压井流水位降深也更加准确，有利于煤矿区矿井的研究，有利于矿井的开采和挖掘。
本发明提供的承压含水层参数反演的计算方法的有益效果是：充分考虑了影响承压－无压转换界与抽水井之间的距离的承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水时间、抽水井流函数、观测井水位和抽水井与观测井之间的距离的变化，使得到的结果更加准确，从而使计算出的无压区扩散系数和无压区单位产水量也更加准确，有利于煤矿区矿井的研究，有利于矿井的开采和挖掘。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
实施例1
一种不稳定承压－无压井流模型的计算方法，包括如下步骤：
(1)、根据承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水速率、抽水时间、水平导水率、抽水井流函数和抽水井与观测井之间的距离计算出承压－无压转换界与抽水井之间的距离和无压区平均水头。
充分考虑了影响承压－无压转换界与抽水井之间的距离和无压区平均水头的各参数承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水速率、抽水时间、水平导水率、抽水井流函数和抽水井与观测井之间的距离的变化，使得到的结果更加准确。
优选地，计算承压－无压转换界与抽水井之间的距离和无压区平均水头的具体方法为：
通过求解第一方程组得到，第一方程组为：
其中，q为抽水速率，m3/s；t为抽水时间；sy为单位产水量；r为承压－无压转换界距离抽水井距离，m；r为抽水井和观测井之间的距离，m；b为含水层厚度，m；kr为水平导水率m/s；h为初始水头；s为给水度；t为导水系数，m2/s；w()为抽水井流函数；
第一方程组的每一个字母都为一个变量，每一个变量发生变化，承压－无压转换界距离抽水井距离r和无压区平均水头
(2)、根据承压－无压转换界与抽水井之间的距离和无压区平均水头计算出承压－无压井流水位降深。
由于上述承压－无压转换界与抽水井之间的距离和无压区平均水头的值计算更加准确，从而使计算出的承压－无压井流水位降深也更加准确。
为了进一步提高准确率，对观测井井流处于承压状态和观测井井流处于无压状态分别进行研究，即判定承压－无压转换界与抽水井之间的距离和抽水井与观测井之间的距离的大小关系，计算出承压－无压井流水位降深。
详细地，当承压－无压转换界与抽水井之间的距离小于抽水井与观测井之间的距离时，即r<r时，则观测井井流处于承压状态，通过求解第一方程得到承压－无压井流水位降深，第一方程为：
其中：h2(r，t)为承压态井流水头，m，即承压－无压井流水位降深。即计算出的h2(r，t)即为承压－无压井流水位降深。
当承压－无压转换界与抽水井之间的距离不小于抽水井与观测井之间的距离时，即r≥r时，则观测井井流处于无压状态，通过求解第二方程得到承压－无压井流水位降深，第二方程为：
其中：h1(r，t)为无压态井流水头，m，即承压－无压井流水位降深。即即计算出的h1(r，t)即为承压－无压井流水位降深。
在考虑在承压流转换到无压流的过程中导水系数、给水度以及扩散系数等水力参数的变化，分别采用二维渗流基本微分方程和潜水流动的布西涅斯克微分方程描述承压区和无压区的地下水运动，建立二维承压－无压井流流场的边界条件和初始条件。
实施例2
一种承压含水层参数反演的计算方法，包括如下步骤：
(1)、根据承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水时间、抽水井流函数、观测井水位和抽水井与观测井之间的距离计算出承压－无压转换界与抽水井之间的距离。
充分考虑了影响承压－无压转换界与抽水井之间的距离的承压含水层的初始水头、含水层厚度、导水系数和单位产水量以及抽水试验的抽水时间、抽水井流函数、观测井水位和抽水井与观测井之间的距离的变化，使得到的结果更加准确。
为了进一步提高准确率，对观测井井流处于承压状态和观测井井流处于无压状态分别进行研究，判定观测井水位与含水层厚度的大小关系，计算出承压－无压转换界与抽水井之间的距离。
详细地，当观测井水位不小于含水层厚度时，即当h′(r1，t)≥b时，观测井井流处于承压状态，通过求解第三方程得到承压－无压转换界与抽水井之间的距离，第三方程为：
当观测井水位小于含水层厚度时，即当h′(r1，t)＜b时，观测井井流处于无压状态，通过求解第四方程得到承压－无压转换界与抽水井之间的距离，第四方程为：
其中，h′(r1，t)为观测井水位，m；q为抽水速率，m3/s；t为抽水时间；r1为抽水井和观测井之间的距离，m；b为含水层厚度，m；kr为水平导水率m/s；h为初始水头；s为给水度；t为导水系数，m2/s；w()为抽水井流函数；r为承压－无压转换界距离抽水井距离，m。
在考虑在承压流转换到无压流的过程中导水系数、给水度以及扩散系数等水力参数的变化，分别采用二维渗流基本微分方程和潜水流动的布西涅斯克微分方程描述承压区和无压区的地下水运动，建立二维承压－无压井流流场的边界条件和初始条件。
(2)、根据承压－无压转换界与抽水井之间的距离计算出无压区扩散率和无压区单位产水量。由于上述承压－无压转换界与抽水井之间的距离的值计算更加准确，从而使计算出的承压－无压井流水位降深也更加准确。
计算无压区单位产水量的具体方法是通过求解第二方程组得到，第二方程组为：
其中，
第二方程组的每一个字母都为一个变量，每一个变量发生变化，承压－无压转换界距离抽水井距离r会发生变化，使其数据更加准确，充分考虑各参数的变化情况。同时，基于泰斯井流假设和裘布依假设，建立改进的二维不稳定承压－无压完整井流解析模型，模型较为简单，更加容易推广。
计算无压区扩散系数的具体方法是通过求解第五方程得到，第五方程为：
其中，d为无压区扩散系数。
本发明提供的一种承压－无压井流水位降深及承压含水层参数的计算方法，(1)在moench和prickett(1972)以及hu和chen(2008)等人的研究基础上，基于泰斯井流假设和裘布依假设，建立改进的二维不稳定承压－无压完整井流解析模型。首先，考虑在承压流转换到无压流的过程中导水系数、给水度以及扩散系数等水力参数的变化，分别采用二维渗流基本微分方程和潜水流动的布西涅斯克微分方程描述承压区和无压区的地下水运动，建立二维承压－无压井流流场的边界条件和初始条件。其次，利用bear(1972)提出的方法线性化布西涅斯克微分方程，采用玻耳兹曼变换(boltzmanntransform)取代吉林斯基势函数(girinskii’spotential)推导数学模型的解析解。
(2)基于降深数据和提出的解析模型，首次推导在单孔等速抽水实验中承压－无压流转换界面动态分布的预测方法以及无压流区扩散系数的计算方法。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。