非稳态气动力

2024-09-20

非稳态气动力（共7篇）

非稳态气动力篇1

摘要：随着日本新干线高速列车最高速度的增加, 隧道内由流体诱发的车体振动已成为研讨的课题之一。本文介绍了改变列车外形降低作用在列车上的非稳态气动力的研究结果。

关键词：高速列车,非稳态气动力,隧道,日本

1 概述

随着日本新干线列车最高速度的增加, 隧道内列车受非稳态气动力作用引发的横向振动, 尤其是横摆振动已成为乘坐舒适性方面一个重要的研讨课题。对各种类型的列车进行实测的结果表明, 气动力具有如下特征[1]:

(1) 非稳态气动力的幅值与相对于列车的气流速度的平方成正比;

(2) 面向隧道壁的列车侧壁出现的持续压力波动是产生非稳态气动力的主要原因;

(3) 压力波动在保持其形状的情况下, 沿顺流方向传播;

(4) 从头车到第6辆~第8辆车压力波动逐渐增大, 然后保持稳定, 最后在列车尾部急剧增大。

运用数值仿真对列车所受的气动力进行了研究。数值计算表明列车底部周围产生的涡流, 以及该涡流沿车辆侧面传播是沿整列车产生压力波动的原因[1]。该研究也揭示了由于列车尾部气流有效截面积的突扩而引起的气流分离是列车尾部产生剧烈压力波动的原因[2]。

已经运用了几种措施来提高隧道内列车的乘坐舒适度。在车辆之间采用了衰减力与车辆间角速度成正比的抗横摆减振器[3]。还开发了主动悬挂系统[4]。这些系统已经在几种系列的日本新干线列车上得以安装并且效果显著。然而, 这些措施只是权宜措施。列车的进一步提速需要减少气动力本身。

改变列车外形对减小作用在车辆上的气动力非常重要。本研究中, 通过线路试验和风洞试验来探究减小气动力的列车外形。首先, 对线路试验数据进行统计分析, 识别出引起非稳态气动力的列车外部部件。然后, 进行风洞试验对列车外形进行分析, 重点关注统计分析识别出的外部部件。最后, 进行线路试验验证根据风洞试验结果提出的列车外形的效果。

2 运用统计分析识别影响气动力的外部部件

2.1 方法

在过去的10年间, 分别对各种型式的列车进行了实车气动测试[1]。测试列车中车辆外形 (例如:车头、车顶和车底) 不同, 并且被测车辆在列车中的位置也不相同。由于有很多参数需要考虑, 因而根据这些数据很难识别出影响气动力的列车部件。

于是, 查找在同一列车上仅外部部件改变情况下的实车试验记录, 并找到了6辆编组列车为期2年的长期试验记录。试验期间, 受电弓和绝缘子罩 (以下称受电弓罩) 、转向架罩等发生了变化 (图1) 。测量了第5辆车两侧的压力。车辆每侧侧面积沿长度方向被分为四等份, 在每部分的中心处布置压力传感器。运用这些测试数据, 力图识别出影响气动力的列车部件。

当列车头车进入隧道时, 将产生压缩波, 当列车尾部进入隧道时, 将产生膨胀波。这些压力波在列车前方以声速传播, 在远端的隧道出口, 压力波被反射并且产生膨胀波朝着隧道入口方向传播。隧道内运行的列车遇到这些压力波时, 压力值将发生突变。在膨胀波经过列车瞬间和列车遇到反射的膨胀波瞬间这段时间内, 列车周围的气流变为准稳定状态。作者估算了该期间内作用在列车车体上的气动力。

运用布置在列车两侧的8个压力传感器的压力测试数据, 作用在车辆上的气动横摆力矩可近似表示如下:

式 (1) 中:Mp———气动横摆力矩;dp———车辆两侧的压力差;si———列车侧面的第i部分的特征面积;li———第i部分的力臂, 等于车辆中心与第i部分中心之间的距离。

以前的研究报道, 当不与迎面而来的列车相遇时, 非稳态气动力的均方根值 (RMS) 与列车速度的平方成正比[1]。运用该关系, 将列车速度250km/h以上的气动横摆力矩的均方根值换算为270km/h速度下的值。统计分析中采用了110组数据。

列车外部部件的改变大部分并不是为了气动研究, 而是有其他目的的。同时, 改变的部件是以各种方式组合在一起的 (表1) 。从而很难估算单个部件对气动力的影响。因此, 作者进行了被称为数量化理论Ⅰ类[5]的统计分析, 来识别引起作用在第5辆车上的非稳态气动力的列车外部部件。数量化理论Ⅰ类将定性数据转换为数字数据并进行多元回归分析。获得了一个作为解释变量的外部部件的组合, 以得出最佳回归线。

2.2 结果

通过采用下列解释变量 (图2) 获得了最佳回归线。

(1) 第5辆车上的受电弓罩;

(2) 中间车的转向架侧罩;

(3) 通道连接处的底部外部风挡;

(4) 转向架与车辆端部之间的侧罩。

自由度调整后的决定系数为0.95[5]。

3 通过风洞试验估测外部部件的气动效应

3.1 方法

如上所述, 根据在有限数量测点获得的压力数据, 近似得出了作用在线路试验真实列车上的气动力[1]。在以前的风洞试验中[6], 采用同样的程序来估测气动力。虽然该估测方法适用于局部压力波动相对缓和的中间车, 但是, 很难确定作用在列车尾部 (该处流场复杂) 的气动力, 也很难估测小型外部部件的影响。因此, 在风洞试验中, 作者力图运用荷重传感器直接测量气动力。

试验在日本铁道综合技术研究所的小型消声风洞进行。喷嘴宽720mm, 高600mm, 试验装置见图3。制作了一个由5辆车组成的新干线列车模型和一个双线隧道模型, 模型比例为1/40。测试在第4辆车和第5辆车上进行。转向架固定在隧道地面上。车内安装了2个荷重传感器 (Minebea U3C-5K-B) 。荷重传感器的一端与转向架连接, 以测量作用在车体上的横向力。气动横摆力矩ML由以下公式得出:

式 (2) 中:Lfront———前转向架处荷重传感器的测量值;Lrear———后转向架处荷重传感器的测量值;d———前转向架与后转向架间距离的一半。

第4辆车的两侧壁均设有4个测压孔。压力传感器 (Scanivalve ZOC23B) 安装在车内, 同时测量每个测压孔处的压力。

列车模型的外部部件, 如列车底部、尾部和罩均是可变的。标准列车模型为底部为方形 (见3.2节 (4) ) , 采用三维长型头部和尾部 (见3.2节 (5) ) , 无通道连接外部风挡、受电弓罩和转向架侧罩的列车模型。外形改变均在标准模型的基础上进行。

在隧道入口上方, 用皮托管测量气流速度。气流速度设定为30 m/s。根据流速和列车高度, 雷诺数 (Re) 为1.6×105。

以前的研究证实[6]:在真实情况正常观察到的振动加速度范围内, 列车振动对流场的影响很小。因此, 本试验中, 列车保持静止状态。

实际上, 列车和隧道之间存在相对运动, 然而风洞试验中两者无相对运动。尽管存在这种差异, 但在风洞试验中, 确定能定性地再现实际现象[6]。认为列车和隧道之间的气流不仅仅是库埃特状流, 而是库埃特流和泊肃叶流的组合[7]。人们认为泊肃叶状流对产生非稳态气动力起重要作用, 因此, 风洞试验能模拟该现象。

3.2 结果和讨论

如果模型的固有频率与气动力的峰值频率接近, 那么很难准确地测量其所受的气动力。因此, 首先要确定荷重传感器能否估测非稳态气动力。将估测的作用在第4辆车上的非稳态气动力与通过对压力数据积分得出的值进行比较 (图4) 。发现两者几乎一致, 事实证实荷重传感器测量系统能估测非稳态气动力。

接下来就运用该测试系统对外部部件, 尤其是统计分析识别出的外部部件的效果进行测试。

(1) 受电弓罩。

上节的统计分析证实受电弓罩对气动力的影响很大。于是估测了3种受电弓罩的影响:二维型、三维型和箱型受电弓罩 (图5) 。三维型受电弓罩的侧板是圆形的, 而二维型受电弓罩的侧板是直的。侧护板可安装在二维型和三维型受电弓罩的侧面。二维型和三维型受电弓罩安装在第2辆车和第4辆车的车顶, 箱型受电弓罩分为两个部分:前部分安装在第4辆车上, 后部分安装在第5辆车上 (图6) 。受电弓本身不安装。

首先研究受电弓罩对顺流场的影响, 测量了下列情况下作用在第4辆车上的气动力:

(a) 在第2辆车上安装二维型受电弓罩和侧护板;

(b) 在第4辆车上安装二维型受电弓罩和侧护板;

情况 (a) 的气动力几乎与无受电弓罩的情况相同。受电弓罩的影响局限在其邻近区域。情况 (c) 的气动力低于情况 (b) 。可能是因为第2辆车上的受电弓罩减小了逼近第4辆车受电弓罩的气流流速。

继而估测受电弓罩形状对作用在第4辆车上的气动力的影响。在第4辆车上分别安装二维型受电弓罩、二维型受电弓罩和侧护板、三维型受电弓罩、三维型受电弓罩和侧护板。箱型受电弓罩分开安装在第4辆车和第5辆车上。同时也对只安装侧护板但无受电弓罩的情况进行了研究。图8说明了受电弓罩形状的影响。安装箱型受电弓罩的情况作用在车辆上的气动力最大, 其次是仅安装侧护板的情况, 第三和第四分别为安装三维型受电弓罩+侧护板和安装二维型受电弓罩+侧护板。这说明侧护板引起了气动力增大。

(2) 转向架侧罩。

上节线路试验的分析说明转向架侧罩有助于减小气动力。为估测其影响, 测量了有无转向架侧罩时作用在第4辆车上的气动力。所有车辆均安装了转向架侧罩。图9证明了转向架侧罩能减小气动力。

(3) 通道连接外部风挡。

线路试验表明, 在通道连接周围存在较大的压力波动[8]。安装外部风挡似乎是一个降低压力波动的有效措施, 上节的统计分析也说明下部外部风挡对气动力有有利影响。因此, 整列车车辆间的间隙处安装由塑料膜制成的外部风挡。为安装方便, 只在模型列车的侧面安装外部风挡。图10显示了外部风挡对作用在第4辆车上的气动力的影响, 外部风挡降低了气动力。

(纵轴表示基于不安装情况气动力横摆力矩均方根值的无量纲值)

(纵轴表示基于无转向架侧罩情况横摆力矩均方根值的无量纲值)

(纵轴表示基于不安装情况气动横摆力矩均方根值的无量纲值)

(4) 车辆底部截面形状。

(纵轴表示基于无外部风挡情况气动横摆力矩均方根值的无量纲值)

数值仿真表明, 列车底部周围区域有涡流产生[1]。认为列车底部形状对于减小气动力起至关重要的作用。

对圆形的底部横截面形状进行了检验 (图11) 。在图11中, 小R的曲率半径为400mm, 大R的曲率半径为800mm, 这与实际列车情况是一样的, 列车所有车辆均采用同一形状。图11显示了作用在第4辆车上的气动力, 大R情况的气动力较小。

也在列车底部安装了导流板 (图12) , 导流板安装采用2种方式:列车的整个侧面安装导流板 (整体导流板) , 或除转向架侧面以外的侧面安装导流板 (部分导流板) 。导流板的高度对应于实际列车上的165 mm高度, 这是车辆限界允许的最大高度。在所有车辆上安装导流板, 测量作用在第4辆车上的气动力。图12说明整体导流板能有效减小气动力, 部分导流板对减小气动力也有一定的作用。

(纵轴表示基于底部方形角情况气动横摆力矩均方根值的无量纲值)

(纵轴表示基于无导流板情况气动横摆力矩均方根值的无量纲值)

(5) 列车的尾部形状。

线路试验说明列车尾端部分压力波动最大[1], 数值流体仿真证明这是由尾部形状影响的气流分离引起的[2]。因此, 测试了尾部形状对作用在尾车上的气动力的影响。试验中准备了5种不同类型的尾部形状:二维短型、二维长型、三维短型、三维长型和方角型 (图13) 。

(纵轴表示基于方角型情况气动横摆力矩均方根值的无量纲值)

短型尾部长度对应于实际列车的3.5m长, 长型尾部长度对应于实际列车的7.5 m长。二维型即所谓的楔形, 三维型的侧面是圆形的, 结果显示三维短形的气动力最大。以前数值研究[2]证明侧面的非稳态气流分离引起尾车的非稳态气动力。推荐运用二维形状抑制侧面的气流分离来降低气动力。现在的试验结果证实了此观点。长型尾部形状减小了由气流分离引起的压力波动, 因此, 二维长型和三维长型的影响几乎没有差异。

4 运用线路试验验证导流板的影响

在真实的列车上安装了上节风洞试验证明有效的导流板, 进行线路试验来验证其影响。

4.1 方法

在4辆车上安装了导流板, 即16辆编组列车的头车到第4辆车。虽然3.2节 (4) 中的风洞试验说明整列车侧面安装整体导流板对于减小气动力有显著作用, 然而, 考虑到转向架的维护, 在除转向架侧面外的侧面安装部分导流板。导流板的高度设为90mm, 这样在车体与转向架间的空气弹簧无气的情况下也不会侵犯车辆限界。通过对车辆每侧4个点测得的压力数据积分, 估算出作用在第4辆车上的气动力。在第4辆车转向架上方走廊地板上安装加速计测量车体的横向加速度。试验装置见图14。

(第5辆及以后车辆未加装导流板)

4.2 结果和讨论

图15显示了气动横摆力矩均方根值和典型功率谱密度分别随列车速度和频率的变化规律。可发现导流板有减小气动力的作用。

图16给出了车体横摆加速度的均方根值和典型功率谱密度。可以看到导流板对车体振动加速度均方根值没有明显效果。原因如下:气动力峰值频率约为3Hz~5Hz, 车体振动频率约为1Hz~2Hz。车体振动在3Hz~5 Hz区间能量很小, 而导流板只减小了约3Hz~5Hz区间气动力的功率谱密度, 因此, 对车体振动几乎没有影响。以前的线路试验验证了从第3辆车到尾车的气动力的峰值频率, 从第3辆车到第6辆~第8辆车峰值频率逐渐减小, 而后直到尾车保持在同一水平[1]。当列车以300km/h的速度运行, 第8辆车以后气动力峰值频率约为2Hz, 与车体振动频率相当。因此, 需要至少从头车到第8辆车安装导流板来探究导流板对车体振动的影响。

5 结束语

通过线路试验和风洞试验研究了降低隧道内作用在列车上的气动力的措施。首先, 对线路试验数据进行统计分析, 识别出引起非稳态气动力的列车外部部件。接着, 通过风洞试验对各种列车外形进行了检验, 重点关注由统计分析识别出的外部部件。风洞试验表明, 车体底部的导流板能有效减小气动力。最后, 运用线路试验来证实导流板的作用。

译自《Journal of Mechanical Systems for Transportation and Logistics》2009, №1, 1~12

参考文献

[1]Suzuki, M., Ido, A., Sakuma, Y.and Kajiyama, H..Full-scale measurement and numerical simulation of flow around high-speed train in tunnel[J].Journal of Mechanical Systems for Transportation and Logistics, 2008, 1 (3) :281-292.

[2]Suzuki, M., Arai, N.and Maeda, T..Numerical simulation of flow around train (2nd report, Unsteady aerodynamic force on tail car) [J].Transactions of the Japan Society of Mechanical Engineers B, 1996, 62 (595) :1069-1067 (in Japanese) .

[3]Fujimoto, H..Lateral Vibration and decreasing measure of it on a Shinkansen train (Decrease of train vibration with yaw damper between cars) [J].Quarterly Report of RTRI, 1995, 36 (3) :122-123.

[4]Sasaki, K., Kamoshita, S.and Simomura, T..Development and field results of semi-active suspension high speed train[J].RTRI Report, 1996, 10 (5) :25-30 (in Japanese) . (其他参考文献略)

利率传导的非稳态过程分析篇2

利率在金融市场中的传导规律是在利率市场化的前提下, 通过利差所形成的货币转移, 使得利率对市场主体产生影响。目前, 有关研究和建立金融市场中反映利率传导过程中非稳态传导基本方程, 从而利用这种传导方程来对利率变动的基本状态 (包括动态特性) 进行理论分析和研究的方法, 在国内有关文献中还几乎没有看到。但我们可以根据国内关于货币政策、利率市场化以及利率传导效应的相关文献和研究论述, 可以得到很多非常好的启发。例如, 夏斌等 (2001) 分析了中国货币政策的中介目标, 提出了由于货币供应量指标已经不适应货币政策中介目标的功能。中央银行要达到的货币政策目标是市场利率而不是货币供应量。从2004年开始, 央行在公开市场操作中的票据发行方式, 已经反映出货币政策的操作目标已经由货币数量转向利率[1];卢宝梅 (2002) 和吴曙明 (1997) 探讨了中国基准利率选择标准和基准利率的传导效应, 认为基准利率有效传导的关键环节是金融市场, 市场化利率体系中的利率传导途径应以基准利率为基础, 通过资本市场利率、信贷市场利率来影响总需求, 并且他们提出以同业拆借利率或国债利率作为基准利率的政策主张[2,3];国家计委经济研究所金融财政研究室课题组 (2003) 和李社环 (2000) 等的研究都提出了利率传导渠道的畅通是货币政策有效性的关键的论点[4,5,6]。总之, 国内已有的很多研究成果, 在理论和实证上为研究利率传导非稳态过程提供了基础和参考。

国外很多学者在对短期利率波动特性分析的基础上, 对于短期利率传导的非稳态特性有很多论述, 可以说, 短期利率的连续时间模型的定量估计, 特别是对于短期利率传导过程中的参数扩散函数的选择, 一直是近期实证金融研究中进展相对较快的领域之一。其中, Bali (1999) 分析、论证并比较了短期利率的连续时间的各类模型[7,8];Chan等 (1992) 提出了一系列不同的短期利率单因素期限结构模型的估计和比较[9]。相关的研究可参见其他学者的研究文献。从国外已有的研究成果中[10,11,12,13,14,15], 可以看出, 短期利率传导的非稳态过程的特点实际上反映了短期利率是随时间变化的长期利率的平均水平, 它体现了利率波动的期限结构, 表明了利率水平的变化取决于短期利率波动过程, 这种波动的当前水平和波动的随机特性都会影响未来利率水平的分布, 而未来利率水平分布则是确定衍生证券价格的一个重要因素。金融市场的实际运行需要确定市场内部的利率分布或确定市场利率达到某一数值所需要的时间。

本文提出和建立描述利率传导非稳态过程的基本方程, 分析利率传导非稳态过程的特性, 确定金融市场中许多实际问题都是一般利率传导过程基本方程所描述的问题的反映, 并且在利率非稳态传导方程的基础上, 通过方程的定解和特征解对利率传导非稳态过程中的利率分布以及利差随利率期限的变化规律进行解析论证。

2 利率传导非稳态过程的特性分析

2.1 利率非稳态传导基本方程

金融市场的实际运行需要确定市场内部的利率场或确定市场利率达到某一数值所需要的时间。利率场是指各时刻市场利率分布的总称, 对于利率传导非稳态过程而言, 利率场是一个非定常利率场, 利率分布随时间而变化。

假设, 在利率传导非稳态过程中, 以x、y、z分别代表不同市场媒介来表示影响利率传导过程的因素, 且相互正交。另外, 其流进市场主体的总货币流量加上市场主体内资金源的生成资金应等于市场主体内资金增量加上流出市场主体的总货币流量。因此, 可以建立一个关于描述利率传导非稳态过程的基本方程:

$\frac{\partial r}{\partial t} = \frac{λ}{ρ c} (\frac{\partial^{2} r}{\partial x^{2}} + \frac{\partial^{2} r}{\partial y^{2}} + \frac{\partial^{2} r}{\partial z^{2}}) + \frac{\frac{d Q}{d t}}{ρ c} (1)$

其中, $\frac{\partial r}{\partial t}$ 为利率变化率, t为时间, ρ为市场容量, c为资金量比, $ρ c \frac{d r}{d t} d x d y d z$ 为市场主体内资金增量, λ为传导率, $\frac{d Q}{d t}$ 为每单位市场主体资金源的生成资金, 因此, 市场主体内的生成资金应为 $\frac{d Q}{d t} d x d y d z . \frac{λ}{ρ c}$ 可以理解为资金扩散率, 它实际上起到了利率传导系数α的作用, α越大, 表明市场内利率均衡的能力就越大, 从而利率在市场中的传导能力就越强。因此, 方程 (1) 也可以表示为:

$\frac{\partial r}{\partial t} = α (\frac{\partial^{2} r}{\partial x^{2}} + \frac{\partial^{2} r}{\partial y^{2}} + \frac{\partial^{2} r}{\partial z^{2}}) + \frac{\frac{d Q}{d t}}{ρ c} (2)$

2.2 利率传导非稳态过程的特性

根据利率非稳态传导方程的基本特点, 可以得出利率传导非稳态过程存在以下3个特性:

①利率r随时间t变化, 当r<t1, 利率r的分布与初始利率有关;当r>t1, 利率r的变化取决于货币增量的条件;

②每一个市场主体的货币量是不同的, 在货币资本的传导过程中, 由于不同市场主体本身对利率变动的灵敏性不同, 因此, 利率变动本身也是积聚或消耗货币的过程。

③从微观层面上看, 与利率传导稳态过程相比, 利率传导非稳态过程表明了每一市场主体的货币流量是处处不相等的。

一般情况下, 传导率是作为利率的函数, 而不是常数。因此, 根据方程 (1) , 利率传导非稳态过程的基本方程也可以表示为:

$\frac{\partial}{\partial x} (λ \frac{\partial r}{\partial x}) + \frac{\partial}{\partial y} (λ \frac{\partial r}{\partial y}) + \frac{\partial}{\partial z} (λ \frac{\partial r}{\partial z}) + \frac{d Q}{d t} = ρ c \frac{\partial r}{\partial t} (3)$

所以, 利率传导过程的微分方程是描述利率传导过程中共性的数学表达式。对于任何传导过程, 无论是稳态的, 还是非稳态的, 也无论是一维的, 还是多维的, 利率传导方程都是普适的。正如F.Black等所述, 金融市场中许多实际问题都是一般利率传导过程基本方程所描述的问题的反映[16,17]。

3 利率传导非稳态过程中的利率分布解析分析

3.1 利率非稳态传导方程的分析解

利率非稳态传导方程的求解问题, 实质上是在规定的初始条件及边界条件下求解利率传导微分方程的过程[18]。初始条件表示了利率变动的初始时刻的利率分布。无论利率传导过程是稳态的还是非稳态的, 其定解条件都存在边界条件。边界条件可以归纳为2类:一类是规定了边界上的利率值, 这种边界为第一类边界条件;另一类是规定了边界上的货币流量密度, 这种边界为第二类边界条件。

解的唯一性定理:如果某一利率函数r (x, y, z, t) 满足利率传导微分方程及特定的初始条件和边界条件, 则此函数就是特定利率传导问题的唯一解。

解的唯一性表明一个解可以表示成等价的几种不同的函数形式[19]。然而, 对于利率传导过程中的具体经济问题, 不仅是求得方程的通解, 而且要求得定解。利率传导方程的定解是既满足表示利率传导非稳态过程的微分方程, 又满足根据问题给出的一系列约束条件的特解。

3.2 债券市场中的利率分布

(1) 利率分布方程及其分离变量求解定理

假设, 在债券市场上影响利率的因素为通货膨胀率i和债券价格P, r0为作为初始利率的债券票面利率, r为传导过程的市场利率, s=r (x, y, t) -r0为市场利差, 则在这个市场上的利率分布r (x, y, t) 可由下列利率传导微分方程和定解条件规定:

${\begin{cases} \frac{\partial s}{\partial t} = α (\frac{\partial^{2} s}{\partial x^{2}} + \frac{\partial^{2} s}{\partial y^{2}}) \\ s (x, y, o) = 1 \\ s (\frac{i}{2}, y, t) + \frac{λ}{α} \frac{\partial s (x, y, t)}{\partial x} |_{x = \frac{i}{2}} = 0 \\ s (x, \frac{Ρ}{2}, t) + \frac{λ}{α} \frac{\partial s (x, y, t)}{\partial y} |_{y = \frac{Ρ}{2}} = 0 \\ \frac{\partial s (x, y, t)}{\partial x} |_{x = 0} = 0 \\ \frac{\partial s (x, y, t)}{\partial y} |_{y = 0} = 0 \end{cases} (4)$

假设在利率影响因素分别为i和P的情况下, X (x, t) 和Y (y, t) 分别为分析解, 且具有同样的定解条件, 即

${\begin{cases} \frac{\partial X}{\partial t} = α \frac{\partial^{2} X}{\partial x^{2}} \\ X (x, o) = 1 \\ \frac{\partial X (x, t)}{\partial x} |_{x = 0} = 0 \\ X (\frac{i}{2}, t) + \frac{λ}{α} \frac{\partial X (x, t)}{\partial x} |_{x = 0} = 0 \end{cases} (5)$

和

${\begin{cases} \frac{\partial Y}{\partial t} = α \frac{\partial^{2} Y}{\partial y^{2}} \\ Y (y, o) = 1 \\ \frac{\partial Y (y, t)}{\partial y} |_{y = 0} = 0 \\ Y (\frac{Ρ}{2}, t) + \frac{λ}{α} \frac{\partial Y (y, t)}{\partial y} |_{y = 0} = 0 \end{cases} (6)$

利率分布方程可以通过分离变量求解定理进行求解。

分离变量求解定理:在第一类边界条件中的边界利率为定值且初始利率为常数的情况下, 利率传导二维微分方程的分析解就是二个利率传导一维微分方程分析解的乘积, 即:s (x, y, t) =X (x, t) ·Y (y, t) 。

(2) 分离变量求解定理的证明

假设s (x, y, t) =X (x, t) ·Y (y, t) 成立, 所以

$\frac{\partial s (x, y, t)}{\partial t} = \frac{\partial}{\partial t} [X (x, t) Y (y, t)] = X \frac{\partial Y}{\partial t} + Y \frac{\partial X}{\partial t}$

又

$\begin{array}{l} \frac{\partial s (x, y, t)}{\partial x} = Y \frac{\partial X}{\partial x} ‚ \frac{\partial^{2} s (x, y, t)}{\partial x^{2}} = Y \frac{\partial^{2} X}{\partial x^{2}} \\ \frac{\partial s (x, y, t)}{\partial y} = X \frac{\partial Y}{\partial y} ‚ \frac{\partial s^{2} (x, y, t)}{\partial y^{2}} = X \frac{\partial^{2} Y}{\partial y^{2}} \end{array}$

所以

$α (\frac{\partial^{2} s (x, y, t)}{\partial x^{2}} + \frac{\partial^{2} s (x, y, t)}{\partial y^{2}}) = α (Y \frac{\partial^{2} X}{\partial x^{2}} + X \frac{\partial^{2} Y}{\partial y^{2}})$

则有

$\begin{array}{l} \frac{\partial s (x, y, t)}{\partial t} - α (\frac{\partial^{2} s (x, y, t)}{\partial x^{2}} + \frac{\partial^{2} s (x, y, t)}{\partial y^{2}}) \\ = X \frac{\partial Y}{\partial t} + Y \frac{\partial X}{\partial t} - α (Y \frac{\partial^{2} X}{\partial x^{2}} + X \frac{\partial^{2} Y}{\partial y^{2}}) \\ = X (\frac{\partial Y}{\partial t} - α \frac{\partial^{2} Y}{\partial y^{2}}) + Y (\frac{\partial X}{\partial t} - α \frac{\partial^{2} X}{\partial x^{2}}) \end{array}$

因为X (x, t) 和Y (y, t) 分别是在 $\frac{i}{2}$ 和 $\frac{Ρ}{2}$ 情况下的一维非稳态利率传导微分方程的分析解, 根据前述, 则有

$\begin{array}{l} \frac{\partial X}{\partial t} = α \frac{\partial^{2} X}{\partial x^{2}} \\ \frac{\partial Y}{\partial t} = α \frac{\partial^{2} Y}{\partial y^{2}} \end{array}$

所以

$\frac{\partial s (x, y, t)}{\partial t} - α (\frac{\partial^{2} s (x, y, t)}{\partial x^{2}} + \frac{\partial^{2} s (x, y, t)}{\partial y^{2}}) = 0$

即

$\frac{\partial s (x, y, t)}{\partial t} = α (\frac{\partial^{2} s (x, y, t)}{\partial x^{2}} + \frac{\partial^{2} s (x, y, t)}{\partial y^{2}})$

又假设s (x, y, 0) =X (x, 0) ·Y (y, 0) 成立, 根据方程组 (5) 和方程组 (6) , 有X (x, o) =1, Y (y, o) =1, 所以s (x, y, 0) =1成立。又因为

$\begin{array}{l} X (\frac{i}{2}, t) Y (y, t) + \frac{λ}{α} Y (y, t) \frac{\partial X (x, t)}{\partial x} |_{x = \frac{i}{2}} \\ = Y (y, t) [X (\frac{i}{2}, t) + \frac{λ}{α} \frac{\partial X (x, t)}{\partial x} |_{x = \frac{i}{2}}] \end{array}$

由于X (x, t) 为一维非稳态解, 所以

$X (\frac{i}{2}, t) + \frac{λ}{α} \frac{\partial X (x, t)}{\partial x} |_{x = \frac{i}{2}} = 0$

所以

$X (\frac{i}{2}, t) \cdot Y (y, t) + \frac{λ}{α} Y (y, t) \frac{\partial X (x, t)}{\partial x} |_{x = \frac{i}{2}} = 0$

这样就有

$s (\frac{i}{2}, y, t) + \frac{λ}{α} Y (y, t) \frac{\partial X (x, t)}{\partial x} |_{x = \frac{i}{2}} = 0$

同理, $s (x, \frac{Ρ}{2}, t) + \frac{λ}{α} Y (y, t) \frac{\partial X (x, t)}{\partial x} |_{x = \frac{Ρ}{2}} = 0$ 。

因此, s (x, y, t) =X (x, t) ·Y (y, t) 可以表示为是方程组 (4) 的定解条件, 从而证明了分离变量求解定理的存在性。

分离变量求解定理表明了求解表示利率传导过程的偏微分方程, 可以通过分离变量法将含有n个变量的偏微分方程化解为n个常微分方程。

(3) 利率分布状况

假设, 市场内部各种利率r仅仅是通货膨胀率i的函数, 这样, 对于利率一维非稳态传导方程的分离变量分析, 实质上就是确定利率传导非稳态过程中市场利率的分布。根据方程组 (4) , 利率一维非稳态传导方程及定解条件如下:

${\begin{cases} \frac{\partial s}{\partial t} = α \frac{\partial^{2} s}{\partial x^{2}} (0 \leq x \leq \frac{i}{2}, t > 0) \\ s (x, o) = s_{0} (0 \leq x \leq \frac{i}{2}) \\ \frac{\partial s (x, t)}{\partial x} |_{x = 0} = 0 \\ α s [(\frac{i}{2}, t)] = λ \frac{\partial s (x, t)}{\partial x} |_{x = \frac{i}{2}} \end{cases} (7)$

根据分离变量求解定理, 方程 (7) 的分析解为s (x, t) =X (x) T (t) , 且有:

$X \frac{d Τ}{d t} = α Τ \frac{d^{2} X}{d x^{2}} (8)$

要使方程 (8) 在x和t的定义域内对任何一个x和t均成立, 则当且仅当

$\begin{array}{l} \frac{1}{α Τ} \frac{d Τ}{d t} = D \\ \frac{1}{X} \frac{d^{2} X}{d x^{2}} = D \end{array}$

成立, 其中, D为一个常数。

D的含义为当t→∞时, 利率传导过程达到稳态。这时, 债券市场所发行的债券的利率将不会随利率期限的变化而变化, 债券利率与基准利率相等。从经济意义上讲, D<0, 否则, 随着时间延长 (或期限的延长) , 债券市场中的利率将无限增高, 这在实际金融市场中是不现实的。

令D=-β2, 则有:

${\begin{cases} \frac{d Τ}{d t} = - α Τ β^{2} \\ \frac{d^{2} X}{d x^{2}} = - β^{2} X \end{cases} (9)$

这样, 方程 (9) 的通解如下:

$\begin{array}{l} Τ = C_{1} e^{- α β^{2} t} \\ X = C_{2} \cos (β x) + C_{3} \sin (β x) \end{array}$

于是有:

$s = e^{- α β^{2} t} [A c o s (β x) + B s i n (β x)]$

其中, A=C1C2, B=C1C3分别为积分常数, 不同的积分常数产生不同的特解。

根据方程组 (7) 所示的边界条件, 再考虑方程组 (9) 的通解, 则有:

$\begin{array}{l} \frac{\partial s (x, t)}{\partial x} |_{x = 0} = e^{- α β^{2} t} [- A β \sin (β x) + B β \cos (β x)] |_{x = 0} \\ = e^{- α β^{2} t} \cdot B β = 0 \end{array}$

所以, B=0。

这样, s (x, t) =Ae-αβ2tcos (βx) 。

再根据方程组 (7) 所示的边界条件, 则另有:

$\begin{array}{l} α s (\frac{i}{2}, t) = α A e^{- α β^{2} t} \cos (β \frac{i}{2}) = λ \frac{\partial s (x, t)}{\partial x} |_{x = \frac{i}{2}} \\ = λ A e^{- α β t} β \sin (β \frac{i}{2}) \end{array}$

所以

$\begin{array}{l} α \cos (β \frac{i}{2}) = λ β \sin (β \frac{i}{2}) \\ \frac{α}{λ β} = \tan (β \frac{i}{2}) (10) \end{array}$

β代表了一维利率传导方程的特征值, 它取决于特定的市场规模, 方程 (10) 本身就是特征方程。假设,

$y_{1} = t g (β \frac{i}{2})$ , 它代表了正切函数曲线

$y_{2} = \frac{α}{λ β}$ , 它代表了双曲函数曲线

由y1和y2的交点决定了β的取值。由于 $y_{1} = \tan (β \frac{i}{2})$ 是以π为周期的函数, 因此, β值将有无穷多个, 即:

${\begin{cases} s_{1} (x, t) = A_{1} e^{- α β_{1}^{2} t} \cos (β_{1} x) \\ s_{2} (x, t) = A_{2} e^{- α β_{2}^{2} t} \cos (β_{2} x) \\ \dots \dots \\ s_{n} (x, t) = A_{n} e^{- α β_{n}^{2} t} \cos (β_{n} x) \end{cases} (11)$

无论An (n=1, 2, …, n) 为何值, s1, s2, …, sn均满足利率传导微分方程及其边界条件。图1表示了利率一维非稳态传导方程特征解的轨迹, 它是一条随通货膨胀率递减的曲线。

由于e-αβ2t的存在, 所以市场中各点的利差是随时间t而降低, 也就是说, 随着利率期限的增大, 利差是呈递减态势。

另一方面, 根据方程组 (4) 和方程 (10) , $\frac{α}{λ}$ 越大, 市场中各点的利差也就越小, 这意味着, $\frac{α}{λ}$ 越大, 利率传导条件就越强, 其利率的传导速度也就越快。这种情况, 市场利率就比较容易趋于均衡和稳定。

4 结论

利率传导非稳态过程同样需要体现进出市场主体的总货币流量相等的概念, 其传导过程具有三个特性, 即, 利率在随时间变化的同时, 它的分布不仅与初始利率有关, 而且它的变化还取决于货币增量的条件;由于不同市场主体的利率灵敏性差异, 因此, 利率的变动也体现了货币积聚和消耗过程;利率传导非稳态过程在微观上表明了每一市场主体的货币流量是不相等的。

表达利率传导过程的微分方程可以说是普适的, 金融市场中的利率变动所引起的传导问题都是一般利率传导过程基本方程的变形。

由于描述利率传导非稳态过程是以微分的形式表示的, 因此其求解过程本身就是对利率传导微观机制的反映。非稳态利率传导问题的分离变量法实质上就是确定传导过程中市场利率的分布, 随着利率期限的增大, 利差是呈递减态势, 利率传导条件是和利率传导速度呈正相关关系, 它决定了利率是否能够较容易地趋于均衡和稳定。

非稳态气动力篇3

单裂隙最基本的概念模型是光滑平行板模型,由Navier-Stokes方程(NS方程)可以推导出裂隙内的稳态渗流符合立方定律[2].然而,实际上单裂隙并不是光滑的平行板,而是具有粗糙的表面,在局部是非平面的、非平行的,并有可能是接触在一起的[3]针对这一问题,许多学者采用非光滑平行板模型来研究单裂隙,这些模型的剖面可以为锯齿形[4]、正弦形[5]和阶梯形[6]裂隙的渗流特性用等效隙宽来表示.如果将粗糙裂隙中的渗流场近似为二维,控制方程可以使用Reynolds方程[7].基于此方程的数值模拟已经被用于求解复杂的单裂隙非稳态渗流场[8,9].但是,Reynolds方程只适合于低雷诺数的情况[10],对于高雷诺数、高流速的非达西和非稳态渗流(一般来说高速非达西渗流产生在高产油气井井筒周围附近区域内[11]),就必须采用三维NS方程来求解[12,13,14].

单裂隙稳态/非稳态渗流模型发展较完善,但真正有实用价值的裂隙网络渗流模型却发展较慢.现有的裂隙网络渗流模型也存在一些缺陷,如只能计算稳态渗流[15];无法解决不连通裂隙(孤立裂隙)所带来的计算不收敛问题[16];把渗透系数人为地变成流速的经验函数,忽略了非稳态渗流的真实物理过程[17].此外,对裂隙网络渗流的求解是复杂的,需要占用大量的计算时间和内存,在现有的计算能力下模拟工程尺度上的渗流是困难的,现有的分析基本局限于小尺度问题[3].

本文重点在于解决裂隙网络的非稳态渗流问题.将裂隙网络分解为许多个单裂隙,每个单裂隙所使用的渗流控制方程由三维NS方程简化而来,但不同于立方定律和Reynolds方程.使用有限差分法进行迭代求解,自由表面的处理采用流体体积法(即VOF法).单裂隙之间通过公共边进行渗流信息交换,经过数学推导,在数值模型中,这些公共边上的渗流都可以由专门的控制方程进行求解,保证数值模型的正确性.同传统方法相比,该非稳态渗流数值模型缩减了计算时间,可以模拟工程尺度上的渗流问题,而且避免了孤立裂隙所带来的影响.本文没有考虑单裂隙中粗糙度的影响,流体的雷诺数不高,属于达西渗流范畴内.

1 单裂隙非稳态渗流

1.1 控制方程

对于不可压缩黏性裂隙流而言,如果不考虑其它假设条件,裂隙内流体满足连续性方程和三维NS方程[18]

式中:u为速度矢量,t为时间,ρ为流体密度,p为压力,v为运动黏性系数.

对于实际工程问题而言,裂隙的长、宽一般都在米级(工程尺度:100～200m),而在厚度方向只有毫米级(工程尺度:<1 cm).数值模拟在解决此类问题时就会遇到多尺度问题,网格划分受到限制,运算速率降低,需要耗费很长的时间才能达到收敛.因此,需要对原模型作以下假设来解决耗时和收敛性问题:(1)假设渗流只在沿平行于裂隙面的方向发生;(2)流速沿厚度方向(z方向)呈抛物线分布,如图1所示.

则流速可以写成以下矢量形式

式中,umax(x,y,t)为最大流速,e为隙宽.

将式(2)代入式(1),并沿着厚度方向对速度积分,则单裂隙渗流控制方程可以写为

对于上述推导,有几点需要进行说明:

(1)公式(3)实质上是把公式(1)所代表的三维渗流问题简化为二维渗流问题,这样数值模拟就可以避免在厚度方向划网格,从而绕开网格划分的多尺度问题;

(2)流速沿厚度方向呈抛物线分布的假设适合于达西流状态.雷诺数Re是评价水流是否进入达西流状态的重要指标,其物理意义是指惯性项与黏滞项的比值.许多学者从试验和数值模拟的角度对此指标进行了研究[19],认为临界Re=1800～3000,当Re大于临界值时,水流即进入非达西流状态;

(3)公式(1)的物理意义表明:流体单位体积上的惯性力等于单位体积上的质量力加上单位体积上应力张量的散度.而公式(3)的物理意义是:把针对流体微元体的应力问题转化为微元体两端的合力平衡问题;

(4)公式(3)中的速度为抛物线速度分布(图1)中的最大速度,它是平均速度的1.5倍;

(5)单裂隙采用有限差分法进行迭代求解,自由表面的处理采用流体体积法,这样就可以较为准确地模拟非稳态渗流过程.

1.2 算例验证

首先要校核该数值模型的正确性和测试其计算效率.对于低雷诺数情况,单裂隙一维渗流达到稳态时,流速的理论解为

式中,utheory为稳态渗流速度,g为重力加速度,J为水力梯度.

王媛等[20]做了低流速单裂隙达西渗流试验,裂缝模型长、宽、厚分别为0.8m,0.1m,0.5mm;进水口压力P1=125 440Pa,出水口压力P0=0 Pa;稳态时水力梯度J=16,雷诺数Re=961.选取稳态后的流速u作为比较参数,数值模拟、理论解和试验结果对比如表1所示.

从表1可以看出,该数值模型的数值解同理论解和试验结果几乎相同,误差在5%以内,表明该数值模型正确.就同一个问题而言,如果用式(1)作为控制方程,依旧采用有限差分法和流体体积法进行迭代求解,计算效率会大幅降低30倍左右.其原因是若采用式(1)为控制方程,该问题就变为了一个三维问题,数值模型中网格尺寸在厚度方向会很小,为避免网格奇异性,就要求网格尺寸在长和宽方向也很小,这必然导致网格数量的增多.而若采用式(3)为控制方程,则网格尺寸只需在长和宽方向进行协调,这样就可以使每个网格的尺寸很大,进而减少计算所需网格.理论上,式(1)可能比式(3)的精度高,但采用式(3)为控制方程的数值模型,牺牲了一点精度却换来了计算效率的大幅提高,而且计算出来的结果,其精度对于工程问题来说足够了.

其次是验证该数值模型的适用性,主要是验证其处理自由表面问题的能力.图2为单裂隙二维非稳态渗流示意图,裂隙长和宽分别为10m和7.5m,厚度为1mm;上端中部和底端的压力分别为P1=29400Pa和P0=0Pa;考虑重力加速度的影响gy=-9.8m/s2;左右两端为封闭边界.应用该数值模型可求得非稳态渗流各个时刻的压力值(图3).

2 裂隙网络非稳态渗流

2.1 控制方程

在裂隙网络中,各单裂隙分别在自身局部坐标下进行迭代求解,在公共边处进行特殊处理.假设在公共边处水流处于瞬时稳定状态,式(3)中第2个方程的速度对时间偏导数项、对流项和扩散项可以忽略,那么公共边处控制方程可以简化为

假设在公共边处,只考虑法线方向的流量平衡,而忽略切向方向的流量交换,则可以得到

式中,I为公共边所连接的单裂隙个数,uk,vk和e分别为第k个单裂隙x方向速度、y方向速度和隙宽.

若公共边法线方向为x方向,则边上的待求压力可写为

式中,pij为公共边待求压力,pi-1,j为x方向相邻节点压力值,dx为x方向网格尺寸,下标k表示括弧内的参数为第k个单裂隙的参数值.

同理可得,当公共边法线方向为y方向时

式中,Pi,j-1为y方向相邻节点压力值,dy为y方向网格尺寸.

这种处理方法保证了公共边处压力的相等,并且使得待求压力与各单裂隙相邻压力节点建立了联系,这样使得公共边处的压力节点也能参与到整场迭代计算中,无需判断公共边处出入流关系.

2.2 算例验证

水流依次从水平薄裂隙(隙宽e=0.5mm)和厚裂隙(隙宽e=1.0 mm)通过,进水口和出水口压力分别为P1=19600Pa和P0=0Pa,如图4所示;图5为裂隙内流速随时间变化图;表2为稳定状态下数值解和理论解的对比.从图4和表2可以看出,该数值模型可以模拟整个非稳态渗流过程,各单裂隙流速和水力梯度的数值解同理论解很好地吻合,说明用该数值模型来处理裂隙网络非稳态渗流问题是符合要求的.

3 工程应用

3.1 水力耦合下岩体的渐进破坏

固体计算模型采用现成的基于连续介质离散元模型[21],结合本文的非稳态渗流数值模型,就可以很容易实现渗流-应力-破坏耦合,模拟岩体在水力耦合作用下的渐进破坏过程.

例如在10 m×10m的岩体区域内随机分布2组裂隙,裂隙的倾向角分别为45°和135°,隙宽服从正态分布,其各个参数如表3所示,由程序自动生成的裂隙网络如图6所示.给定裂隙网络的边界条件为:上部z=10m处压力边界p1=19600Pa,底部z=0m处压力边界p0=0Pa,左右两端为封闭边界.裂隙网络的压力分布如图7所示,随着水流的非稳态渗流,该耦合模型可以完整呈现裂纹扩展和渐进破坏过程,由于水流无法进入孤立裂隙,因此孤立裂隙的压力值为0.同Tang等[22]的研究成果相比,由于实现了裂隙岩体的非稳态渗流模拟,该模型比较真实地再现了岩体渐进破坏过程,为水压致裂或者油气开发等领域的相关问题打下了基础.

3.2 滑坡灾害

大量工程实例发现,对于库区堆积层滑坡来说,滑体物质组成松软,易渗水;基岩物质组成坚硬,不透水;中间有一条非常明显的滑带,易积水.库水涨落和降雨所产生的水流通过坡体表面的裂缝迅速沿着滑带流动,可以把滑带看成是裂隙介质,研究滑带的性质显得特别有意义.利用本文的非稳态渗流数值模型可以实现两个目的:(1)估算滑带的平均渗透系数;(2)反分析滑带内部压力分布.

如果知道降雨强度和渗流时间就可以估算滑带的平均渗透系数.以茅坪滑坡为例,根据文献[23]所示的滑坡剖面图,就可以建立其计算模型(如图8)-将滑带简化为多块光滑平行板组成的裂隙网络,假设这些平行板的厚度相同,调节这一厚度,使得其满足以下已知条件:降雨强度为50mm/d的情况下,水流沿着滑带从坡顶到坡底的渗流时间为3h.结果表明,当厚度为1mm时满足已知条件,这时的平均渗透系数为0.8mm/s.

在获得滑带平均渗透系数的基础上,根据外部的可测物理量可以反分析降雨停止后滑带内部压力分布.位于高程360m处的竖井揭示,在井深40m处开始出现低渗透层,并有承压水.井深40m处承压水约1m,至41～42m处承压水头上升至6m左右,基岩出现在42m处[24].竖井发现的6m承压水就是一个外部可测物理量,满足这一已知条件的结果就是所求的结果.

一般情况下,当降雨量较大时,底部出水口来不及将所有水都排出,导致坡体内部的压力上升.根据这一现象,数值模型中就需要微调出口处的渗透系数(主要是降低出水口流量),具体说来就是调节出口处板的厚度.由于只是调节出口处板的渗透系数,所以整体上不会影响平均渗透系数的计算结果.结算表明,当出口处板的厚度为0.9mm时,计算结果可以满足可测物理量的观测结果,这时滑带内部的压力分布如图9所示。

4 结论

(1)本文发展了一种新的裂隙岩体非稳态渗流数值模型.该模型将裂隙网络分解为多个单裂隙,每个单裂隙所使用的渗流控制方程由三维NS方程简化而来,简化后的控制方程将三维渗流问题转化为二维平面渗流问题,避免了网格划分的多尺度问题,提高了计算效率.有限差分法和流体体积法的采用使该模型可以模拟水流在裂隙内的非稳态渗流过程;

(2)经过数学推导,裂隙网络中公共边上的渗流可以由专门的控制方程进行求解,使得公共边处的压力节点能参与到整场迭代计算中,既能保证数值模型的正确性,又无需判断公共边处出入流关系;

非稳态气动力篇4

本文介绍了非稳态接触滚动力学的发展历程。在钢轨短波波磨的成因中,非稳态滚动接触起着重要的作用。

在Grassie和Kalousek报告的概述中[1],波长介于2 cm～10 cm的钢轨短波波磨被认为是一种所谓“波长固定机理”依然未知的波磨。钢轨短波波磨的波长总是在2 cm～10 cm的原因到底是什么呢?本文介绍了非稳态接触力学是如何解释钢轨短波波磨的“波长固定机理”,以及其他一些悬而未决的问题。

像许多滚动接触问题一样,非稳态接触力学也是由Kalker最先提出的,尽管他没有使用“非稳态滚动”这个词,而使用了“瞬时滚动”这一概念[2,3]。Kalker研究了车轮在法向载荷作用下,从静止到稳态滚动的过渡过程中接触斑的应力分布变化。这被认为是从Cattaneo[4]稳态解向Carter[5]稳态解的过渡,在第3节将针对这个问题进行详细的讨论。

然而,当人们对随时间变化极快的滚动接触现象感兴趣时,也处处需要使用非稳态接触力学。这方面的例子有:通过曲线时车轮的啸叫声,车轮在具有短波长不平顺的轨道上滚动等。

2 非稳态滚动接触

关于轮轨非稳态接触的最早论文是Kalker在1970年—1971年间发表的[2,3]。对于其他类似的非稳态过程,可追溯到更早的时候。早在 1924 年就发现了飞机空气动力学中的非稳态问题[6],在1931年发现了轮胎和路面的非稳态滚动接触问题[7]。这两种情况都很好地展示了非稳态的含义。

非稳态(包括非线性的非稳态)空气动力学最早是由Wagner着手研究的[6]。Teichmann的论文介绍了关于非稳态空气动力学引起机翼颤振的研究成果[8]。在飞机的空气动力学中,往往不是使用“非稳态”这个术语,而是使用“动态”(相对于“静态”,即相对于“稳态”)。在机翼的空气动力学研究领域,人们比在接触力学领域更早地了解了非稳态的含义(图 1)。当空气粒子沿着机翼流动时,如果空气的流线相对机翼几乎保持不变的情况下,这时的空气动力学问题可当作稳态过程来处理。然而,如果当空气粒子通过机翼时,机翼位置和空气的流线同时产生显著的变化,这种状态就称为非稳态过程。它的特征参数与运动的波长L(如流动速度、频率)以及机翼的厚度a有关。

类似的想法对于轮胎和路面的接触也是相通的(图2)[9,10]。如果接触斑的位置和形状以及接触斑附近的速度场保持不变,那么这样的滚动过程就可以当作稳态滚动来处理。但是,如果当接触斑内的单元在滚动接触过程中,接触斑的位置和形状以及接触应力场或蠕滑发生显著变化,那么所遇到的就是非稳态滚动接触问题。它的特征参数与运动的波长L以及接触斑的长度a有关。

这个定义可以直接引用到轮轨接触的情况。如果当接触斑内的一个质点通过接触斑的过程中,接触斑的形状或其他接触参量发生了明显的变化,那么这个过程就称之为非稳态滚动接触。在图3中,正弦变化分别用L/a=4或L/a=16表示。同样地,它的特征参数与波长L和滚动方向上的接触半径a有关。对于L/a=16来说,接触参量仅仅产生了微小的变化,相反地,对于较小的波长L/a=4,则可以看到明显的变化。对后一种情况,采用非稳态理论进行分析是十分必要的。

下面对一些典型的轮轨运动进行考察,以确定是否需要用到非稳态接触力学。对于踏面接触情况,接触斑在滚动方向上的长度在1 cm～2 cm之间,因此,接触半径a的最大值是1 cm。

(1) 对于曲线通过,在圆曲线上运动的有效波长L趋向无限大。

(2) 对于蛇行运动,以及轮对由直线进入曲线过程中即在缓和曲线上的瞬时运动,运动的波长在几米左右,典型的取L=10 m。

(3) 对于不圆车轮在钢轨上滚动,运动波长取决于踏面不平顺的周期。对于周期为1的不圆车轮,波长L大约是3 m;对于周期为6的多边形车轮,运动波长L=50 cm。

(4) 对于钢轨短波波磨,典型的运动波长L=4 cm。

(5) 啸叫噪声的波长L为1 cm(频率为2 000 Hz,车辆速度为72 km/h)。

综合上述情况:对曲线通过,运动波长比L/a =∞;对于蛇行运动,波长比L/a=1 000;对于多边形车轮,波长比L/a=50;对于短波钢轨波磨,波长比L/a=4;对于曲线啸叫,L/a=1(甚至更小)。根据Kalker和Groβ-Thebing的研究,当运动波长比L/a<10时,必须考虑非稳态效应。因此,对于在短波波磨的钢轨上运行的车轮以及曲线啸叫情况,必须使用非稳态滚动接触力学分析。

如果L/a值较大,那么理所当然接触斑和蠕滑率的变化也要考虑在内,但这时非稳态分析不是必要的。根据Kalker理论,把这种情况当作一种连续的稳定状态进行分析就足够了。

3 非稳态轮轨接触力学和波磨成因的研究现状

3.1 非稳态滚动接触力学

之前已经提到,Kalker最先分析了瞬态(也称非稳态)接触力学[2,3]。Kalker将其限于二维问题,并且仅考虑车轮在法向和切向载荷作用下,从静止到稳态滚动的过渡过程。最近,Gonzáles和Abascal使用边界元法和稍微不同的模型重新检验Kalker的研究成果,其结果几乎一致。下面介绍Gonzáles的研究成果(图4)。

上述分析中,假定切向力Q达到μN最大力值的75%。η(τ)表示质点的位置,它从τ=0的位置进入接触斑。η(τ)/aH= 2意味着质点通过了整个接触斑,和η(τ)/aH= ∞几乎没有任何差别。因此,进行非稳态分析时,考虑到η(τ)/aH= 4就足够了,这时其与稳态解几乎一样。

Groβ-Thebing完成了谐波激励的线性化非稳态分析[12,13],这对应于车轮在有微小谐波不平顺的钢轨上运行时的工况。线性化分析意味着所考虑的轨道不平顺(即波磨的幅值),以及所产生的法向力和蠕滑力都是微小变化的。通过线性化的假设,可分析波磨的初始成因——因为波磨的初始幅值也很小。但是参考状态未必是微小的。作为参考值,较大的法向力或者较大的蠕滑率(如曲线通过时)都是允许的。许多关于钢轨波磨成因的论文都是基于这个线性化理论发展的[14,15,16,17,18,19]。

在一份未发表的论文中,Kalker已经开始考虑简化理论(注:Kalker滚动接触的简化理论)[20]是否能够被用来解决非稳态滚动问题。Shen和Li开发了一种算法以加快解算速度[21]。但是后来,Kalker在一次口头交流中认为,使用通常的简化理论来解决非稳态问题,得不到可被接受的解。最近,Alonso和Giménes[22]再一次研究了这个问题,他们通过修改Winkler在简化理论基础研究工作中的柔性系数定义形式,获得了可以接受的结果。这个方法不仅可以分析钢轨波磨的初始成因——正如Groβ-Thebing和其他人所做的,也可以用来分析波磨的成长。

对于不圆车轮的滚动运动,在L/a值足够大的情况下,使用稳态接触力学已经足够了。这里,笔者仅提供4篇文献,即1份文献综述[23]和3份德文论文[24,25,26]。在曲线啸叫方面没有相应的论文,只查到使用稳态接触力学分析该问题的论文。

3.2 钢轨短波波磨的形成

人们研究钢轨波磨已经历经100多年,并且发表了大量的文献(参见Birmann[27]和Krabbendam[28]的综述论文)。

关于钢轨波磨成因的论文将会在随后的章节中介绍。这里,笔者给出2个最近的测试结果。在Baumann的博士论文中[29],对钢轨波磨的断面进行了冶金学研究。文中有两张图非常有趣。第一张图是钢轨的波磨形态图(图5),第二张是相应的轮廓曲线(图6)。波峰之间的距离大概是3.3 cm,波磨的深度小于30 μm。图5中明亮的区域(所谓的白色侵蚀层)对应的是波峰,而灰色区域对应的是波谷——它是由明亮区域被磨耗而形成的。短波波磨的波距可能是变化的,但它总是在2 cm～8 cm的范围内。波磨的波形不具有严格的周期性或者谐振性(见图6),但是它的准周期形态还是很明显的。

4 钢轨波磨的成因和线性接触模型的基本概念

众所周知,钢轨波磨的成因可以由短期的动力学和长期损伤的反馈环路来解释。对短波波磨,这种损伤机理就是磨耗[1](图7)。短期的动力学将轮轨结构动力学和接触力学作为线性模块处理,给出了轮轨接触面的接触几何变量、蠕滑率和接触力。

在接触力学中,接触点处的轮轨变形作为刚体位移来考虑(图8)。另一方面,接触力学的变量在结构力学中定义为一个接触点[13]。这个假设仅当结构振动的波长 LS相对于接触斑在滚动方向上的尺寸非常大的情况下才是有效的(2a是接触椭圆的直径)。接触斑的长度通常为1 cm,对直径为1 m的车轮,当其具有三个节线的模态振型(图9,左图), LS/(2a)的值大约是50;对于轨道动力学中主导的“pinned-pinned”模态振型(见图9,右图),LS/(2a)的值大约是120。需要明确指出的是:关系式LS/(2a)与之前提到的L/(2a)是两个概念,两者无关。

在反馈环路中,非稳态接触力学有两个重要的任务。首先,必须计算出接触力的波动——它造成轮轨的高频振动。其次,计算出车轮在具有不平顺的钢轨上滚动时,轮轨之间的磨耗率。不仅仅是接触应力的幅值,还有相位偏移以及对应于谐波轮廓不平顺所产生的磨耗,都必须使用正确的物理方式加以确定。

瞬态接触力和磨耗率可以由Kalker的非稳态和非线性理论确定[30]。但是由于接触斑必须被恰如其分地离散成网格,所以计算起来非常耗时。如果把轮轨关系的大量参数考虑在内,这些参数又会在很大的范围内变化,那么Kalker的完全精确理论对预测钢轨波磨就不再合适了。稳态、非线性的参考状态和非稳态线性化(图10)相分离的概念,能够使计算得以简化,因而对于波磨预测的参数研究可以在较少的计算次数里得到解决[13、16]。

5 非稳态接触力学的处理

5.1 法向接触问题

轮轨间的激励是通过接触面的不平顺产生的。在滚动方向上,表面几何形状用周期性的基本函数来描述,如图11所示。在此,使用独立的自由度来表达高度、梯度、曲率的变化。从这些基本函数中,可以推导出主曲率半径及轮轨之间的穿透深度,这些都是Hertz接触理论的输入参数[19]。

如果一个车轮在具有极短波长(相对于接触尺寸)的周期性不平顺轨道上滚动时,仅仅会引起非常小的法向力变化,这是因为在接触斑的一部分会产生法向压力的正波动,而另一部分会产生负波动。因此,轮轨接触面上极短波长的表面粗糙度对接触力的影响非常小。随着粗糙度波长的减小,这个影响逐渐减弱。因此,不平顺的幅值必须乘上Remington所给出的比例因子[31]。

另外有人用非稳态接触力学对水平方向曲率变化的影响进行研究,并得出了类似的结论。如图12所示,对于极小的轮廓不平顺,两者的效果相当于滤波器。所描述的滤波器支配着非稳态工况的法向接触问题,并且能够当作稳态工况来处理。这一点与非稳态的切向接触问题不同。

5.2 切向接触问题

在Hertz理论可以使用的条件下,法向问题依然可以得到解析解。对于切向接触力的计算,不管是稳态问题还是非稳态问题,只能得到数值解。切向接触问题表明了接触力(纵向、横向以及自旋力矩)和蠕滑率(纵向、横向以及自旋蠕滑)之间的关系。作为输入参数,参考状态的接触斑半轴长、车轮载荷以及摩擦系数,都是在求解法向力的过程中要用到的。切向接触的参考状态和非稳态问题,可由Kalker的滚动接触理论来确定[30]。对于参考状态,将用到非线性稳态方程以及线性化的非稳态方程来表达。

稳态条件下蠕滑力和蠕滑率的非线性关系见图13。诚然,对于稳态工况,这种关系可以在参考状态点进行线性化(Tξ、νξ)处理。蠕滑力的波动分量ΔTξ和蠕滑率的波动分量Δνξ是同相位的。曲线的切线斜率就是Kalker蠕滑系数。稳态滚动接触的线性模型可以用在接触点位置的阻尼器表示(见图14的左图)。

如图13所示,接触过程的非稳态处理对于蠕滑力和蠕滑率之间的线性化关系有着重要的影响。在蠕滑力ΔTξ和蠕滑率Δνξ之间出现了相位偏移,非稳态工况的蠕滑力振幅明显小于稳态工况的蠕滑力振幅。因为振幅和相位取决于激励频率f=v/L,通过确定速度v和波长L,就可以得到频率响应特性。使用关系式 a/L = (f a)/v来代替频率f是更理想的。例如,图15给出随a/L变化的横向蠕滑力的频率响应特性。类似于法向接触问题的Remington效应,当一个质点穿过接触区域时,它将先后表现为正蠕滑率和负蠕滑率。蠕滑力振幅随着波长的减小而持续减小。另外可以观察到,在波长很短时,蠕滑力和蠕滑率之间的相位偏移达到-90°。

频率响应特性的一个可能的力学模型是弹簧和阻尼的串联(见图14的右图)。对于稳态状态(a/L=0),只有阻尼发挥作用,并且能够补偿轮轨间的整体相对速度。对于短波长运动,阻尼器变成了刚性,而只有弹簧起作用。蠕滑力和蠕滑率波动分量之间的这种关系,可以解释为Kalker蠕滑系数与频率有关。在非稳态接触力学中还有更多这样的系数。例如,法向力的变化导致了切向力的波动。然而,人们不再使用蠕滑系数这个词,而说得到了切向接触力学的复杂线性系数。

5.3 非稳态接触力学的磨耗速率

车轮在微小波磨钢轨上滚动时,钢轨的磨耗是怎样计算的呢?摩擦功假说认为[32,33,34],根据接触力和蠕滑率,接触区域的磨耗率可以假定与摩擦功成比例。对于参考状态,在钢轨长度方向上的磨耗率是恒定不变的。对于非稳态线性化部分,则要考虑到摩擦功在纵向方向上随着接触变量而变化。对接触面单元的局部磨耗率,摩擦功密度的计算如下:对一个钢轨上的质点,将其通过整个接触斑历程中的摩擦功进行积分,即得到车轮通过一次时的摩擦功密度。既然已经用法向接触力学为波磨钢轨表面定义了基本函数,这个方法同样也可以被用来计算磨耗率。

磨耗率预测的结果可以用轨道不平顺的自由度来描述。钢轨不平顺的振幅、相位、梯度和曲率会引起频率响应,该值取决于L/a。对于长为1 cm的接触斑,3个特征函数(结合自由度)的磨耗率系数如图16所示。图16中,正值表示轮廓幅值的增长,负值则表示减少。当然,非稳态效应并不可能将材料添加到钢轨上。应当考虑到参考状态下,或多或少的材料总会被磨耗掉。显然,仅有波长为2 cm～10 cm范围内的周期性不平顺的幅值才会增强。因此,非稳态接触力学会产生类似波长滤波器的效果。

到此为止,这只关注了接触力学过程,轮轨之间的结构动力学同样需要关注。对于钢轨轨头的磨耗过程,轨道的结构动力学至关重要。例如,考虑离散支撑的钢轨(UIC60轨)受到作用在轨枕跨中或轨枕正上方的谐波垂向力激励[35]。图17给出了钢轨前三阶谐振和第一阶反谐振振型以及对应的阻抗。其中“pinned-pinned”振型是阻尼比相对较小的振型,这类振型的节点正好位于轨枕处,轨道似乎被固定在了轨枕处。在轨枕跨中距施加激励,很容易激起这类振型,并且频率响应出现一个显著的峰值。而在轨枕上方施加激励,则很难激起这类振型,并且频率响应对应于一个谷值(反谐振振型)。

波长与反谐振振型频率相对应的轨道不平顺,会引起极端的法向力波动,结果导致极端的准周期性磨耗。相反地,如果轨道不平顺的波长与谐振振型频率相对应,则钢轨的垂向振动不会引起上述问题。

与之相反的是与较高的基准横向蠕滑率(例如由轮对蛇行或通过曲线所引起)相关的钢轨横向振动的谐振。它们会引起较大的蠕滑率波动,另一方面也会产生较高的、准周期性的钢轨磨耗。

与车辆运行速度相关的磨耗系数的波长滤波器(见图16),可以转换为钢轨横向和垂向阻抗(图18(b)和图18(c))的频率滤波器(图18(a))。如果横向传递函数的谐振峰值或者垂向传递函数的反谐振谷值与磨耗率系数的最大值相吻合,那么临界情况就会出现。在速度为180 km/h和100 km/h时会出现“pinned-pinned”反谐振振型,第一阶“pinned-pinned”横向谐振振型对应于70 km/h,第二阶谐振振型对应于180 km/h。在这2个工况中,运动波长都在2 cm～10 cm之间。对于这些波长,磨损率系数出现了最大值,由此形成了钢轨波磨。因此,波长滤波器效应决定了钢轨波磨波长固定的机理。

5.4 线性接触单元的限制条件

在应用非稳态接触力学的线性化模型时,应首先检查线性模型是否以及何时超出了限制条件。以下列出了各变量的极值,超过这些范围就需要采用非线性模型:

(1) 接触斑半径的波动总是小于参考状态的接触斑半径;

(2) 轮廓曲率的波动总是小于车轮的曲率。如果超过了该限制,车轮就不再与波磨的波谷相适应。但这不是受限于线性化模型,而是Hertz接触力学的有效性界限;

(3) 法向力的波动和法向压强的波动必须小于参考状态的对应值,以确保轮轨接触;

(4) 当前的蠕滑力,即参考值和非稳态值的总和,必须小于摩擦系数和当前法向力的乘积。

5.5 线性化模型的验证

在第4节的结尾,提到Kalker非线性接触理论可以被用来验证线性化模型的结果。图19给出了L/a=4、谐波蠕滑率波动时,两种方法对应的计算结果。稳态切向接触力达到了车轮载荷和摩擦系数乘积的40%。正如前文对车轮从静止到滚动过渡的解释,仅当运动波动等于接触斑半径的两倍时,才需要使用非线性模型。上述结论对于从稳态状态向谐波变化的过渡过程也是同样有效的。

从图19可以看出,两种方法的结果几乎是相同的。对于由较高的波动蠕滑率产生的蠕滑力,不能被转化为磨耗计算。对于线性非稳态计算,在整个仿真过程中,参考状态接触斑都存在粘着和滑行区域。如果存在较高的蠕滑率波动,那么线性单元的磨耗预测就需用非线性计算来校核。

6 线路运行结果的比较

可以将图16的计算结果与测试结果进行比较。1982年,Widmayer[36]在德国北部靠近梅灵—Haspelmoor的一条长度将近3 km的德国铁路(DB)波磨试验线上,测得了波磨波长的统计分布(图20)。尽管数值分析的结果和各种测试结果之间的细节比较有待考证,但仍然可以注意到两者有着惊人的吻合。最大波磨对应的波长为4 cm,并且找不到波长小于1 cm的波磨,但存在极少数波长大于10 cm的波磨。长期的比较需要考虑轨道的参数和运行条件。其他的测试结果,如Frederick在英国铁路[34]所作的测试,也要被考虑在内。

7 理论和实践的总结

使用非稳态接触力学的分析结果,可以解释为什么钢轨短波波磨的波长分布在2 cm～10 cm之间。这种波长固定机理是一种波长滤波器,它来自非稳态接触效应。其中哪一波段的波长会得到最大程度的发展取决于轨道结构动力学。当钢轨垂向反谐振振型或者横向谐振振型与车速耦合时,会导致波长3 cm～4 cm的波磨成长,而波长滤波器的最大值就处于这个区间。

有什么补救方法呢?

(1) 改变波长滤波器是不太现实的,因为它是车轮在钢轨上滚动时所产生非稳态滚动接触效应所导致的。车轮和钢轨的参数,如载荷、车轮直径等只能在很小的范围内变化。因此,只能接受图16所示波长滤波器,并且其正值出现在2 cm～10 cm范围内。

(2) 由于钢轨短波波磨具有特定的形式,当考虑是否有可能消除在2 cm～10 cm波长范围内的钢轨剧烈磨损时,必须考虑结构力学。因为钢轨“pinned-pinned”振动的反谐振振型恰恰具有上述特性。

(3) 每个补救方法都要着眼于削弱反谐振的影响,并且减少产生波磨的敏感性。一种可能的方法是采用钢轨连续支承,或者减小轨枕的间距。每种方法是否能够改变钢轨在反谐振附近的频率响应特性,例如使用钢轨衬垫,都需要通过数值分析来验证。

(4) 波磨的形成与钢轨的最初不平顺相关,也与行车速度相关,反过来也限制了列车运行的速度。因此,在钢轨投入使用前进行打磨是很有必要的。对已经产生波磨的钢轨轨头打磨时需要特别当心,仅仅磨出光滑的几何表面是不够的。因为产生波磨的钢轨往往同时产生了塑性变形,这导致钢轨表面产生了残余应力。对残余应力存在的表层也必须打磨掉,否则打磨之前的不平顺又会重新生成。这个流程在许多铁路公司的钢轨打磨操作章程中都有表述。

8 待解决的问题

与波磨形式相关的非线性问题,在很大程度上并没有解决。在之前已经提到,这主要是需要考虑塑性变形的问题。根据钢轨的冶金学研究,即使中等波磨的钢轨也会产生塑性变形。因此,在波磨初始成因中,除了磨耗,塑性变形是第二大损伤机制。根据当前的认知水平,进行昂贵的有限元分析是不可避免的,其中的主要问题是为承受周期循环载荷作用的钢轨找到一种合适的材料,文献[42]的第一个研究成果令人鼓舞。

参考文献

[1] Grassie,S.L.and Kalousek,J..Rail corrugations.Characteristic,causes and treatments[J].Journal of Rail and Rapid Transit,Proc.Instn.Mech.Engrs,Part F,1993,207,57-68.

[2]Kalker,J.J..Transient phenomena in two elastic cylinders rollingover each other with dry friction[J].Journal of Applied Mechan-ics,1970,37,677-688.

[3]Kalker,J.J..Transient rolling contact phenomena[C].Proceed-ings of the 25th meeting of the American Society of LubricationEngineers(ASLE),Park Ridge,IL,1971,177-184.

非稳态气动力篇5

目前水平井是开发低渗透油藏的主要技术。正确预测低渗透油藏中水平井的产能对水平井的优化设计具有重要的意义。

目前低渗透油藏水平井非稳态产能的研究主要采用解析方法或者数值方法,在解析方法方面应用比较广泛的有源函数方法[1]、傅里叶变换与拉普拉斯变换结合的方法[2,3]、源函数与拉普拉斯变换相结合的方法[4]等。国内对水平井非稳态渗流研究较早的宋付权,等人[5]在假设水平井为线汇的情况下,利用水平井的椭球流态推导出水平井三维渗流时的无因次井底压力公式。杨正明,等人[6]采用数值方法对考虑启动压力梯度的直井的非线性渗流数学模型进行了求解。

低渗透油藏由于存在启动压力梯度,导致其渗流特征与中高渗油藏不同。目前低渗透油藏水平井非稳态产能研究中部分学者没有考虑启动压力梯度,部分学者考虑了启动压力梯度,但采用数值方法进行求解。现建立考虑启动压力梯度的低渗透油藏水平井不稳定渗流数学模型,采用拉普拉斯变换及数值反演技术对水平井的非稳态产能进行研究。

1 低渗透油藏中水平井开发时的数学模型

模型假设:(1)顶底封闭的油藏中心有一口水平井;(2)考虑启动压力梯度;(3)流体做单相流动;(4)不考虑井筒中的压力降。

考虑启动压力梯度的低渗透油藏中的运动方程为

$v = \frac{k}{μ} (\frac{d p}{d r} - G) (1)$

平面径向渗流时的连续方程为

$\frac{ρ v}{r} + \frac{\partial (ρ v)}{\partial r} = \frac{\partial (ρ φ)}{\partial t} (2)$

岩石和流体的状态方程分别为

$ϕ = ϕ_{0} + C_{f} (p - p_{i}) (3)$

ρ=ρ0exp[CL(p-pi)]=ρ0[1+CL(p-pi)] (4)

将式(1)、式(3)和式(4)代入式(2)中,得到基本的渗流微分方程

$\frac{\partial^{2} p}{\partial r^{2}} + (\frac{1}{r} - G C_{L}) \frac{\partial p}{\partial r} - \frac{G}{r} = \frac{1}{η} \frac{\partial p}{\partial t} (5)$

式(5)中,η=k/(ϕμCt),是导压系数。

模型的初始条件为: p|t=0=pi (6)

内边界条件为: $\frac{\partial p}{\partial r} |_{r = r_{w}} = \frac{q μ B}{2 π L r_{w} k} + G (7)$

外边界条件为: p(∞,t)=pi (8)

定义无量纲的量: $p_{D} = \frac{2 π k L}{q μ B} [p_{i} - p (r, t)]$ , $r_{D} = \frac{r}{r_{w}} ‚ t_{D} = η t / r_{w}^{2} ‚ G_{D} = \frac{2 π k L r_{w}}{q μ B} G$ ,并且由于 $G C_{L} < < \frac{1}{r}$ ,可以忽略不计。故上述模型写为无量纲量的形式为

${\begin{cases} \frac{\partial^{2} p_{D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{D}}{\partial r_{D}} + \frac{G_{D}}{r_{D}} = \frac{\partial p_{D}}{\partial t_{D}} \\ (r_{D} > = 0, t_{D} > 0) \\ p_{D} (r_{D}, 0) = 0 \\ \frac{\partial p_{D}}{\partial r_{D}} |_{r_{D} = 1} = - (1 + G_{D}) \\ p_{D} (\infty, t) = 0 \end{cases} (9)$

以上各式中,L为水平段长度,cm;v为渗流速度,cm/s;k为渗透率,μm2;μ为流体粘度,mPa·s;p为压力,10-1MPa;G为平均启动压力梯度,10-1MPa/cm;ρ为流体密度,g/cm3;rw为井径,cm;t为时间,s;Cf、CL分别为岩石和流体的压缩系数,1/10-1MPa。

2 低渗透油藏水平井非稳态产能计算

对式(9)的数学模型进行拉普拉斯变换[7],有:

$\begin{array}{l} \int_{0}^{\infty} [\frac{\partial^{2} p_{D}}{\partial r_{D}^{2}} + \frac{1}{r_{D}} \frac{\partial p_{D}}{\partial r_{D}} + \frac{G_{D}}{r_{D}}] e^{- s t_{D}} d t_{D} = \\ \int_{0}^{\infty} \frac{\partial p_{D}}{\partial t_{D}} e^{- s t_{D}} d t_{D} (10) \end{array}$

所以,渗流方程变换为:

$\frac{d^{2} \bar{p_{D}}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d \bar{p_{D}}}{d r_{D}} + \frac{G_{D}}{s r_{D}} = s \bar{p_{D}} (11)$

于是,拉氏空间的数学模型为

${\begin{cases} \frac{d^{2} \bar{p_{D}}}{d r_{D}^{2}} + \frac{1}{r_{D}} \frac{d \bar{p_{D}}}{d r_{D}} + \frac{G_{D}}{s r_{D}} = s \bar{p_{D}} \\ (r_{D} \geq 0, s > 0) \\ \frac{d \bar{p_{D}}}{d r_{D}} |_{r_{D} = 1} = - (1 + G_{D}) / s \\ \bar{p_{D}} (\infty, t) = 0 \end{cases} (12)$

求解得到拉氏空间的压降解为

$\bar{p_{D}} (r_{D}, s) = \frac{Κ_{0} (\sqrt{s} r_{D})}{Κ_{1} (\sqrt{s})} [\frac{1 + G_{D}}{\sqrt{s^{3}}} + F (1)] +$

$\bar{p_{D}}_{p} (r_{D}, s) (13)$

采用Laplace数值反演中的Stehfest方法[8]对式(13)进行求解,得到非稳态产能。

3 计算分析

3.1 水平井定产量生产时井底压力的变化

水平井定产量生产时,井底无因次压力的变化见图1。可以看出,在非稳态渗流过程中,保持井产量不变时,渗流初期,井筒压力下降较快。经过一段时间,压力出现拐点(b点),此后,压力一直保持平稳下降。随着无因次启动压力梯度的增加,无因次压力逐渐减小;这说明在产量一定的情况下,无因次启动压力梯度越大时,需要消耗更多的压力来克服启动压力梯度。

3.2 水平井定压生产时产量的变化

无因次产量与无因次时间的关系曲线见图2和图3。可以看出,产量随时间的延长而降低。在生产初期,产量下降速度较快,这一下降速度一直维持到拐点,拐点以后产量随时间延长平稳下降。从拐点以后产量变化规律(图3)可以看出启动压力梯度对水平井产量有较大的影响,启动压力梯度越大,水平井产量越小,这一影响在开发生产的后期显的尤为严重。

4 结论

(1)建立了考虑启动压力梯度的水平井非稳态渗流数学模型,采用拉普拉斯变换进行求解,得到了拉氏空间解,并采用数值反演技术得到其非稳态产能。

(2)计算分析表明,水平井产量随着开发的进行先是急剧下降,经过 “拐点” 之后,水平井的产量平稳下降。启动压力梯度对水平井产量有较大的影响,启动压力梯度越大,水平井产量越小。

摘要：用水平井开发低渗透油藏可以建立有效的驱替压差,提高单井产能。正确预测低渗透油藏中水平井的产能对水平井的优化设计具有重要意义。针对低渗油藏的特征,考虑启动压力梯度,建立了低渗透油藏中水平井开发时的不稳定渗流数学模型,采用拉普拉斯变换及数值反演技术,得到低渗透油藏中水平井的非稳态产能,并分析了无因次井底压力和无因次产量的变化规律及影响因素。研究成果对低渗透油藏水平井的开发具有重要的指导意义。

关键词：低渗透油藏,水平井,非稳态产能,拉普拉斯变换

参考文献

[1] Gringarten A C,Ramey H J Jr.The Use of Source and Green's Func-tions to Solve Unsteady-Flow Problems in Reservoirs,SPE3818

[2] Goode P A,Thambynayagam R K M.Pressure drawdown and buildupanalysis of horizontal wells in anisotropic media.SPE Formation Eval-uation,1987;2(4):683—697

[3] Jalal F O,Djebbar T.Transient pressure behavior of Bingham non-Newtonian fluids for horizontal wells.Journal of Petroleum Scienceand Engineering,2008;61,21—32

[4] Ozkan E,Raghavan,R.Some new solutions to solve problems in welltest analysis:part I-analytical considerations,SPE 18615

[5]宋付权,刘慈群.低渗油藏中含启动压力梯度水平井生产动态.西安石油学院学报,1999;14(3):11—14

[6]杨正明,于荣泽,苏致新,等.特低渗透油藏非线性渗流数值模拟.石油勘探与开发,2010;37(1):94—98

[7]张建国,雷光,张艳玉.油气层渗流力学.东营:中国石油大学出版社,2006:49

非稳态气动力篇6

巴什托普油田位于塔里木盆地西南坳陷麦盖提斜坡西段群苦恰克构造带, 主要含油层位为石炭系生屑灰岩段, 储集空间以裂缝—孔隙型为主, 储层横向分布稳定, 夹层不发育。生屑灰岩油藏为一准层状油气藏, 上油下水, 二者之间无隔层, 油水界面呈东高西低倾斜分布, 东西两端的油水界面差异达109m。

一般认为产生倾斜油水界面的原因主要有3种:水动力、毛细管压力和构造运动影响。

1.1 水动力

若水动力作用造成巴什托普油田生屑灰岩油藏倾斜油水界面, 则流过油藏东西方向截面的地层水流量为:

沿倾斜面的地层水渗流速度为:

Vw为地层水渗流速度, m/d;K为储集层渗透率, 10-3μm2;μw为地层水黏度, m Pa·s;ρw, ρo为地层水和地层原油密度, g/cm3;g为重力加速度, 9.8m/s2;θ为油水界面倾角, 度;qw为地层水流量, m3/d;Lw, Hw为地层水流过的油层横向长度和纵向厚度, m。

通过计算地层水流量需达到89.04×104m3/d才能形成倾斜油水界面, 该流量远远大于目前的日产液量48m3/d, 说明水动力不是造成巴什托普油田生屑灰岩油藏油水界面倾斜的原因。

1.2 毛细管压力

由毛细管力作用形成的油水过渡带厚度为:

△h为油水过渡带厚度, m;ps, pd为接近束缚水条件的毛管压力和排驱压力, MPa。

通过公式计算油水过渡带厚度为1.42~3.21m, 该厚度远远小于油水界面的倾斜幅度109m。因此, 毛细管力也不是造成巴什托普油田生屑灰岩油藏油水界面倾斜的原因。

1.3 新构造运动

精细古构造研究表明, 在更新世之前, 巴什托普油田为一西高东低的背斜构造, 构造高点位于西部 (见图1) , 有利于捕集油气形成古油藏。构造西部的群6X井实钻生屑灰岩段为水层, 但岩心中构造缝可见沥青充填, 由此可以证实先期存在古油藏。麦盖提构造主要存在3期成藏, 即海西早期、海西晚期和喜马拉雅期。其中海西早期油藏多被破坏, 海西晚期油藏保存有好有差, 后期又有调整, 喜马拉雅期主要为油气藏形成与古油藏调整期。后期受喜马拉雅造山运动的影响, 巴什托普油田构造东部迅速变陡, 地层整体由北西上倾转变为北东翘倾, 溢出点抬高, 古油藏遭到破坏, 油气开始由西向东调整运移。

现今构造呈东高西低的形态, 构造高点位于琼002井附近, 按照传统石油地质学理论, 高点部位的油气密度应该小于低部位的油气密度。但目前巴什托普油田生屑灰岩油藏油气密度分布却相反, 原油和天然气密度呈现由西向东逐渐增大的趋势, 西部的群5X井原油密度为0.807g/cm3, 天然气相对密度为0.72g/cm3, 而东面的琼002井原油密度为0.815 g/cm3, 天然气相对密度为0.89g/cm3, 这也说明目前巴什托普油田生屑灰岩油藏目前正处于调整中, 还未达到稳定状态。

综合分析认为, 现今巴什托普油田生屑灰岩油藏是海西期形成的古油藏受喜马拉雅晚期新构造运动影响遭受破坏, 同时高点逐渐向东偏移, 油气由西向东运移调整形成的, 受新构造运动影响, 生屑灰岩油藏目前仍处于调整中, 未达到稳定状态。

2 巴什托普油田生屑灰岩油藏是否是一个连续的油藏

地层对比生屑灰岩储层连片分布, 前期钻井资料认为生屑灰岩油藏是一个连续油藏。在群古1X井与琼002X井之间部署的群501X井实钻发现生屑灰岩段顶部储层物性非常好, 但试油却为水层。而构造低部位的群5X、群5-1X、群古1X井及构造高部位的琼002X井生屑灰岩段均为油层, 这与构造特征明显不符, 说明巴什托普油田生屑灰岩油藏不是一个连续的油藏。

原油分析表明, 生屑灰岩油藏西部区域蜡质含量少, 粘度低, 东部区域蜡质含量高, 粘度高。地层水矿化度分析表明, 西面的群6、群5井矿化度高 (12.78~20.59×104mg/L) , 地层水为Ca Cl2水型, 反映封闭的水体环境;东面的群古1X、群501X、琼002井矿化度低 (4.14~7.78×104mg/L) , 且地层水为Na2SO4水型, 反映开放的水体环境。流体性质的差异也反映巴什托普油田生屑灰岩油藏不是一个连续的油藏。

群古1X井的生屑灰岩储层物性差, 平均孔隙度在3.0%左右, 而其他井的平均孔隙度为7.2%~12.3%, 且投产后产量及地层能量迅速下降, 半年时间累积产油只有270t, 也说明群古1X井储层物性差。在生屑灰岩段均方根振幅图 (见图2) 上, 群古1X井以西存在一个均方根振幅较强的区域, 反映该区域存在一片致密带。综合分析认为在群古1X井以西区域存在一个物性遮挡带, 将巴什托普油田生屑灰岩油藏分成两部分。

3 巴什托普油田生屑灰岩油藏演化过程

从以上分析来看, 巴什托普油田生屑灰岩油藏是一个典型的非稳态油藏。受新构造运动的影响, 古油藏被破坏, 油气开始由西向东运移重新聚集成藏。在运移过程中, 由于受到群古1井以西物性遮挡带的影响, 油气在群古1井以西区域聚集成藏, 成为油气勘探开发的有利区域。西面的群5X井、群5-1X井累积产油分别为2.02×104t、3.63×104t, 生产过程中地层能量下降缓慢, 说明群古1X井以西区域油气资源丰富。群古1X井以东区域地层水为Na2SO4水型, 反映为开放的水体环境, 说明该区无有效遮挡, 油气散逸流失, 难以聚集成藏。构造低部位和储层物性好的区域, 油气散失快, 高部位和储层物性差的区域油气散失慢, 剩余部分残余油 (见图3) 。

参考文献

非稳态气动力篇7

常用的太阳能热水系统按介质流动形式主要分为封闭式和直流式。其中直流式太阳热水系统是由我国科学家率先提出的, 在国际上有着重要的影响。该系统是传热工质一次流过集热器加热后进入储水箱或用水点的非循环系统。直流式系统又分为热虹吸型和定温放水型。为了使热水温度达到用户的使用要求, 通常选用定温放水式。在出口处装有测温元件以保证出口热水温度。选用定温放水型直流式太阳热水系统中吸热板芯为翼管式的平板集热器为研究对象, 如图1、2所示。集中讨论了管之间的中心距和管径对于集热器的效率和出水量的影响。

1 性能分析的基本理论

1.1 非稳态下集热器的效率

undefined

其中 (Mc) w为水的总热容, kJ/kg;Tfa为水初始温度, ℃;Tfb为水的最终温度, ℃;A为集热器的采光面积, m2;undefined为平均辐照度, W/m2, Δτ为时间间隔, s。

1.2 直流式太阳能热水系统的产水量

产水量=集热器加热冷水的次数×每次加热的水量 (2)

2 建立模型及条件

在家用条件下出水温定为50℃, 可用于盥洗、沐浴和洗涤用, 完全可以满足用户对水温的要求。初始冷水温度为10℃, 模拟的是定温放水型直流式系统, 当集热器内水温低于50℃时, 冷水停留在集热器中吸收太阳能被加热直至达到指定温度, 此时集热器处于闷晒状态;当水温达到50℃时, 水阀打开, 热水流入储水箱。集热器吸热板厚度1mm, 盖板到吸热板距离40mm, 玻璃盖板发射率为0.88, 集热器长度为1m, 宽度随管间距变化。环境温度为293K。

3 模拟结果及分析

首先给出5组不同的管中心距, 在非稳态条件下, 看不同的管中心距对集热器瞬时效率的影响。计算结果见表1。太阳辐射强度如图3所示。

由表1可知, 管径为20mm时, 随着管中心距减小, 集热器效率上升, 当管中心距为50mm时, 集热器效率最高, 每平米的出水量最多。经过对数据的比较分析, 在不同的太阳辐照度下, 不同管中心距对应的集热板温度的变化趋势是一样的, 因此, 仅以太阳辐照度为700的情况下为例, 如图4所示。

集热器的集热板吸收太阳辐射, 通过管翼把热量传递给管中流体使流体升温。由于管径相同, 管内流体所能接受管翼向其传导的热量相同, 根据传热学中的傅里叶定律, 管翼长度减小, 使得吸热板与管中流体的温差也降低。因为管翼长度小接受太阳能辐射少, 在相同的时间下管中流体相对于长管翼的集热器管中流体的温升较小, 所以相应吸热板温度比较低。这样就出现了图4所示的情况, 随着管中心距的减小, 吸热板温度降低, 单次出水用时较长, 总出水量少。而集热板温度又对平板集热器的热损有重要的影响, 集热板温度低, 热损小;相反, 吸热板温度高, 集热器热损增加, 这就导致出现了管翼长度大的集热器效率低的情况。

从表1还可以看出, 当管间距相同时, 管径大的集热器效率高, 每平米出水量大。同理, 也以太阳辐照度为700的情况下为例来说明这种现象。如图5所示, 管径小的集热器集热板温度高于管径大的集热器集热板温度, 并且单次出水用时明显较短。出现这种现象的原因是管径小, 使得管内流体质量小, 那么管翼能向里传递的热量减小, 当集热器的管翼长度相同, 吸收太阳辐射量相同。这样造成了管径小的集热器集热板温度较高, 热损增加, 使得集热器效率降低。并且在相同时间情况下, 质量小的流体温升大, 管内流体温度分布均匀。

4 结论

主要讨论了直流式定温放水系统中在非稳态条件下管中心距和管径对集热器效率的影响。经过以上对数据的分析可以看出, 仅从热力学角度来讲, 集热器在管径相同的情况下, 管中心距越小, 对提高集热器效率越有利, 每平米集热器出水量多, 更能够满足用户的需求。提供的5组数据中, 管中心距50mm被认为最佳。当采用最佳管中心距时, 管径大的集热器效率高于管径小的集热器效率。

摘要：对平板型太阳能集热器的直流工作方式的非稳态传热作了分析, 建立了数值传热模拟模型, 利用fluent6.0软件对其作了数值模拟分析, 探讨了直流式平板集热器在非稳态传热中的排管间距、管径的集热器结构优化参数。模拟结果表明, 集热器在管径相同的情况下, 管中心距越小, 集热器效率越高, 本研究提供的5组数据中, 管中心距50mm被认为最佳。当采用最佳管中心距时, 管径大的集热器效率高于管径小的集热器效率。

【非稳态气动力】推荐阅读：

稳态模型07-17

稳态试验07-08

稳态模拟08-23

稳态诱发电位09-02

稳态分析论文06-19

稳态Kalman滤波08-08

稳态作业参考答案09-15

《人体的稳态》说课论文07-01

内环境的稳态教学设计05-15

对“内环境与稳态”一节教学的几点建议05-10

>> 查看更多相关文档