0 引言
水驱是中国油田最基础、应用最广泛的主体开发方式。随着开发的不断深入，主力油田普遍进入开发中后期，处于产量快速递减阶段，开发效益变差，面临剩余油有效开采需求与技术水平不适应等一系列世界级难题[1⇓-3]，亟需构建“油藏-工程”高效生产体系，将精细化理念贯穿各生产环节，实现分层注水实时监测和层段调整，利用实时监测的生产数据结合油藏静态数据不断强化对地质的再认识[4]。刘合等[5-6]突破了井下高压永置式流量检测、长效动密封等核心技术，研制出中国第4代分层注水工艺技术，包括缆控式、波码通信两套主体技术，实现了中国分层注水井“边注边测边调”的数字化升级跨越，现已规模应用。波码通信分层注水工艺技术无需在油管外预置电缆，可满足带压作业要求，以注入流体为信号载体，通过人工干预地面控制阀或井下层段调节阀实现井筒无线双向通信，调节层段注入量并获取生产数据[7]，该工艺在长庆油田水驱开发的低渗透—特低渗透油藏取得了显著效果。
相关学者开展了井筒流体脉冲响应模型的研究[8]。其中国外学者针对管道流体系统进行了深入分析，明确了长距离悬垂管道流体输送机理[9-10]，相关成果在深海油气立管输送过程中得以应用[11]，但其工况与分层注水存在较大差异。中国相关研究主要关注作业过程，分析了井筒流体脉冲传输参数与井身结构、层段吸水特性之间的内在规律[12-13]，定量描述数据上传、下传的脉冲响应特性，然而研究成果侧重于理想流态下稳态压力、流量工作点的解释，未考虑流体黏度、井筒结构参数与储集层渗流因素对流态的影响[14]，且缺乏井筒流体动态响应过程的描述与分析，井筒空间各点的压力、流量状态响应机理尚不明确[15]。波码通信式分层注水工艺技术编解码核心参数仍主要依赖于工程经验，需要开展注水井井筒流体动态响应机理研究，为进一步提高通信效率提供理论依据。
针对分层注水井井筒实际工况，本文对井筒流体微元进行受力及运动状态的分析，建立流体微元的压力-流速关系模型，将微元模型拓展至全井空间流场，结合井筒流体波传输的时域特性，求解井筒流体流速、压力等参数的动态时空分布。在此基础上开展注水井井筒流体动态响应分析，与无阻力损失条件下的传统伯努利方程模型进行对比，并且明确沿程阻力损失、传输延迟、信号衰减等因素对注水井筒无线波码通信过程的影响，进一步验证该模型的准确性。
1 注水井井筒动态流体力学模型
1.1 基于流体力学的压力-流速微元分析
首先根据流体微元的空间位置与受力状态，推导流体微元的合应力矢量矩阵，然后根据流体微元间的相对运动状态，结合牛顿内摩擦定律，引入流体黏度，得出应力张量与变形速度张量的本构关系，最后根据牛顿力学定律求解流体微元流速场与压力场的关系。
1.1.1 流体微元受力分析
将宏观管路流体进行网格微元化处理。以井筒流体流动方向为横轴建立三维直角坐标系，定义微元六面体的速度矩阵与加速度矩阵分别为：
$\begin{matrix} ~u={{\left[ \begin{matrix} {{u}_{x}} & {{u}_{y}} & {{u}_{z}} \\ \end{matrix} \right]}^{\text{T}}} \\ \end{matrix}$
$a={{\left[ \begin{matrix} {{a}_{x}} & {{a}_{y}} & {{a}_{z}} \\ \end{matrix} \right]}^{\text{T}}}$
根据牛顿第二定律，可知流体微元的合力矩阵为：
$\begin{matrix} ~F=ma=\rho dxdydz\left[ \begin{matrix} {{a}_{x}} \\ {{a}_{y}} \\ {{a}_{z}} \\\end{matrix} \right]=\rho dV\left[ \begin{matrix} {d{{u}_{x}}}/{dt}\; \\ {d{{u}_{y}}}/{dt}\; \\ {d{{u}_{z}}}/{dt}\; \\\end{matrix} \right] \\\end{matrix}$
设定向井流体微元与水平方向夹角为α，且流动为均匀流动，流体微元仅受重力和表面力作用，其中微元自身重力施加到微元各个面的压力矩阵为：
$\begin{matrix} ~G=\rho gdV{{\left[ \begin{matrix} sin\alpha & 0 & -\cos \alpha \\ \end{matrix} \right]}^{\text{T}}} \\ \end{matrix}$
流体微元在各方向上都受到相邻微元的作用力，以流动方向为主方向，流体微元受力分析如图1 所示。其中σxb、σxf分别表示轴向后方和前方微元产生的垂直于受力面的法向压应力，τxy、τxz分别代表x方向上的形变对y、z方向产生的切向应力，其他受力面及各方向的受力状态以此类推。进一步开展流体微元表面力的构成分析，以x轴方向为例，双侧法向应力梯度的方向相反，可得出σxf和σxb的表达式如下：
同理可推导出y和z方向流体微元所受的法向应力σyf，σyb，σzf，σzb，进而得出流体微元的法向应力矩阵σp为（7）式：
${{\sigma }_{p}}=\left[ \begin{matrix} {{\sigma }_{xb}}-{{\sigma }_{xf}} & 0 & 0 \\ 0 & {{\sigma }_{yb}}-{{\sigma }_{yf}} & 0 \\ 0 & 0 & {{\sigma }_{zb}}-{{\sigma }_{zf}} \\ \end{matrix} \right]= \mathbf{diag}\left( dx\frac{\partial {{\sigma }_{x}}}{\partial x},dy\frac{\partial {{\sigma }_{y}}}{\partial y},dz\frac{\partial {{\sigma }_{z}}}{\partial z} \right)$
流体微元沿x轴方向受到相邻微元的切向应力，双侧应力梯度方向相反，顶面与底面的切向应力τuyz和τlyz如（8）式和（9）式所示：
$~\begin{matrix} {{\tau }_{uyz}}={{\tau }_{yz}}-\frac{\partial {{\tau }_{yz}}dy}{2\partial y} \\ \end{matrix}$
${{\tau }_{lyz}}={{\tau }_{yz}}+\frac{\partial {{\tau }_{yz}}dy}{2\partial y}$
则流体微元的切向应力矩阵τq为：
${{\tau }_{q}}=\left[ \begin{matrix} 0 & {{\tau }_{lyx}}-{{\tau }_{uyx}} & {{\tau }_{lzx}}-{{\tau }_{uzx}} \\ {{\tau }_{lxy}}-{{\tau }_{uxy}} & 0 & {{\tau }_{lzy}}-{{\tau }_{uzy}} \\ {{\tau }_{lxz}}-{{\tau }_{uxz}} & {{\tau }_{lyz}}-{{\tau }_{uyz}} & 0 \\\end{matrix} \right]=$ $\left[ \begin{matrix} \text{d}x\text{ d}y\text{ d}z \\ \text{d}x\text{ d}y\text{ d}z \\ \text{d}x\text{ d}y\text{ d}z \\\end{matrix} \right]\left[ \begin{matrix} 0 & \frac{\partial {{\tau }_{yx}}}{\partial y} & \frac{\partial {{\tau }_{zx}}}{\partial z} \\ \frac{\partial {{\tau }_{xy}}}{\partial x} & 0 & \frac{\partial {{\tau }_{zy}}}{\partial z} \\ \frac{\partial {{\tau }_{xz}}}{\partial x} & \frac{\partial {{\tau }_{yz}}}{\partial y} & 0 \\\end{matrix} \right]$
流体微元所受合应力由法向应力和切向应力构成，其矢量矩阵eστ为：
${{e}_{\sigma \tau }}={{\sigma }_{p}}+{{\tau }_{q}}=\left[ \begin{matrix} \text{d}x\text{ d}y\text{ d}z \\ \text{d}x\text{ d}y\text{ d}z \\ \text{d}x\text{ d}y\text{ d}z \\\end{matrix} \right]\left[ \begin{matrix} \frac{\partial {{\sigma }_{x}}}{\partial x} & \frac{\partial {{\tau }_{yx}}}{\partial y} & \frac{\partial {{\tau }_{zx}}}{\partial z} \\ \frac{\partial {{\tau }_{xy}}}{\partial x} & \frac{\partial {{\sigma }_{y}}}{\partial y} & \frac{\partial {{\tau }_{zy}}}{\partial z} \\ \frac{\partial {{\tau }_{xz}}}{\partial x} & \frac{\partial {{\tau }_{yz}}}{\partial y} & \frac{\partial {{\sigma }_{z}}}{\partial z} \\\end{matrix} \right]$
依据牛顿力学第二定律将（11）式与（4）式联立可得流体微元所受合力F：
$\boldsymbol{F}=\left(\left[\begin{array}{l} \mathrm{d} y \mathrm{~d} z \\ \mathrm{~d} x \mathrm{~d} z \\ \mathrm{~d} x \mathrm{~d} y \end{array}\right]^{\mathrm{T}} \boldsymbol{e}_{\sigma \tau}\right)^{\mathrm{T}}+\boldsymbol{G}=\left[\begin{array}{c} \frac{\partial \sigma_{x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}+\rho g \sin \alpha \\ \frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \sigma_{y}}{\partial y}+\frac{\partial \tau_{z y}}{\partial z} \\ \frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \sigma_{z}}{\partial z}-\rho g \cos \alpha \end{array}\right] \mathrm{d} V$
1.1.2 流体微元速度-压力关系求解
井筒流体输送过程中，各流体微元存在相对运动，从而产生内摩擦力。为构建压力分布与速度分布之间的关系，需要明确应力张量与变形速度张量之间的本构关系。根据牛顿内摩擦定律，可得流体微元的本构切向应力矩阵eτ：
$\begin{matrix} {{e}_{\tau }}=\mu \left[ \begin{matrix} 0 & \frac{\partial {{u}_{x}}}{\partial y}+\frac{\partial {{u}_{y}}}{\partial x} & \frac{\partial {{u}_{x}}}{\partial z}+\frac{\partial {{u}_{z}}}{\partial x} \\ \frac{\partial {{u}_{x}}}{\partial y}+\frac{\partial {{u}_{y}}}{\partial x} & 0 & \frac{\partial {{u}_{z}}}{\partial y}+\frac{\partial {{u}_{y}}}{\partial z} \\ \frac{\partial {{u}_{x}}}{\partial z}+\frac{\partial {{u}_{z}}}{\partial x} & \frac{\partial {{u}_{z}}}{\partial y}+\frac{\partial {{u}_{y}}}{\partial z} & 0 \\\end{matrix} \right] \\ \end{matrix}$
定义eσ、$\overline{{{e}_{\sigma }}}$分别为流体微元运动、静止状态下的法向应力矩阵，根据牛顿流体特征可知法向应力与法向变形速度张量成线性关系，微元静止状态下各受力面所受法向应力相同（σp），则：
$\left\{ \begin{array}{*{35}{l}} {{e}_{\sigma }}={{a}_{0}}\left( {{e}_{\sigma x}}+{{e}_{\sigma y}}+{{e}_{\sigma z}} \right)I+{{b}_{0}}\mu \mathsf{diag}\left( \frac{\partial {{u}_{x}}}{\partial x},\frac{\partial {{u}_{y}}}{\partial y},\frac{\partial {{u}_{z}}}{\partial z} \right) \\ \overline{{{e}_{\sigma }}}=-{{\sigma }_{p}}\mathbf{I} \\\end{array} \right.$
静止状态下流体微元不受切向应力作用且流速的散度为零，根据（14）式推导可得（15）式：
$\mathsf{diag}\left[ \left( 1-{{a}_{0}} \right){{e}_{\sigma x}}-a\left( {{e}_{\sigma y}}+{{e}_{\sigma z}} \right),\left( 1-{{a}_{0}} \right){{e}_{\sigma y}}- \right.$ $\left. {{a}_{0}}\left( {{e}_{\sigma x}}+{{e}_{\sigma z}} \right),\left( 1-{{a}_{0}} \right){{e}_{\sigma z}}-{{a}_{0}}\left( {{e}_{\sigma x}}+{{e}_{\sigma y}} \right) \right]=\mathsf{0}$
求解（15）式，将结果代入（14）式。由于流体微元的合应力由法向应力和切向应力构成，结合（13）式得：
${{e}_{\sigma \tau }}={{e}_{\tau }}+{{e}_{\sigma }}=\mu \mathsf{diag}\left[ \frac{4\partial {{u}_{x}}}{3\partial x}-\frac{2}{3}\left( \frac{\partial {{u}_{x}}}{\partial y}+\frac{\partial {{u}_{x}}}{\partial z} \right)-\frac{{{\sigma }_{p}}}{\mu }, \right.$$\frac{4\partial {{u}_{y}}}{3\partial y}-\frac{2}{3}\left( \frac{\partial {{u}_{y}}}{\partial x}+\frac{\partial {{u}_{y}}}{\partial z} \right)-\frac{{{\sigma }_{p}}}{\mu },\text{ }\frac{4\partial {{u}_{z}}}{3\partial z}-$ $\left. \frac{2}{3}\left( \frac{\partial {{u}_{z}}}{\partial x}+\frac{\partial {{u}_{z}}}{\partial y} \right)-\frac{{{\sigma }_{p}}}{\mu } \right]+{{e}_{\tau }}$
将（16）式与（3）式、（12）式联立，求解得到流体微元流速场与压力场的关系（见（17）式），为后续井筒流场的精确分析提供理论基础。
$F=\rho dV\frac{\partial u}{\partial t}=$ $dV\left[ \begin{matrix} \mu gra{{d}^{2}}u+\frac{\mu divu}{3\partial x}-\frac{\partial {{\sigma }_{p}}}{\partial x}+\rho gsin\alpha \\ \mu gra{{d}^{2}}u+\frac{\mu divu}{3\partial y}-\frac{\partial {{\sigma }_{p}}}{\partial y} \\ \mu gra{{d}^{2}}u+\frac{\mu divu}{3\partial z}-\frac{\partial {{\sigma }_{p}}}{\partial z}-\rho gcos\alpha \\ \end{matrix} \right]$
1.2 井筒流场空间的压力-流速分析
考虑注水作业实际工况，首先将三维流体微元流速-压力关系表示为柱坐标形式，并根据层流的流动连续性方程，结合井筒结构与地层参数，求解井筒流场的压力分布与流速分布。
1.2.1 基于流体力学的空间流场分析
首先进行流态分析，黏性流体流动的雷诺数方程如（18）式所示。实际注水流量一般较小，流动速度较慢，流体流动雷诺数远小于湍流临界值，可视为层流。层流过程流体密度变化极小，可忽略不计，因此本文不考虑管柱流体的可压缩性。由此得出注水工况下流态满足以下条件：①由于压力脉冲幅度远小于背景压力，连续管道中为定常流场；②管柱流体为层流流动；③忽略生产时管柱流体的可压缩性。
$\begin{matrix} Re=\frac{\rho uD}{\mu }=\frac{4\rho {{q}_{\text{v}}}}{\mu \pi D} \\ \end{matrix}$
由于微元内流体流入、流出质量差与单位时间段内微元质量变化相同，进一步推导流体流动的连续性方程：
$\frac{\partial \left( \rho {{u}_{x}} \right)}{\partial x}dV+\frac{\partial \left( \rho {{u}_{y}} \right)}{\partial y}dV+\frac{\partial \left( \rho {{u}_{z}} \right)}{\partial z}dV=$ $\left( \frac{\partial {{u}_{x}}}{\partial x}+\frac{\partial {{u}_{y}}}{\partial y}+\frac{\partial {{u}_{z}}}{\partial z} \right)\rho dV$
注水井井筒流道为圆形回旋体，为准确分析实际工况下流体在流道中的流动状态，首先需要将三维直角坐标系转化为柱坐标系（见图2 ），根据直角坐标系与柱坐标系下的速度转换关系，得出柱坐标系下切向（θ）、径向（r）、轴向（x）流速梯度的表达式：
将连续性方程代入流速-压力关系式，得到忽略流体可压缩性时流体微元流速梯度的表达式：
$gradu={{\left[ \begin{matrix} {{u}_{x}}\frac{\partial }{\partial r}\left( \frac{\partial r}{\partial x} \right) \\ {{u}_{r}}\left[ \frac{\partial }{\partial r}\left( \frac{\partial r}{\partial y} \right)+\frac{\partial }{\partial \theta }\left( \frac{\partial \theta }{\partial y} \right) \right]cos\theta \\ {{u}_{r}}\left[ \frac{\partial }{\partial r}\left( \frac{\partial r}{\partial z} \right)+\frac{\partial }{\partial \theta }\left( \frac{\partial \theta }{\partial z} \right) \right]sin\theta \\ \end{matrix} \right]}^{T}}=$ ${{\left[ \begin{matrix} {{u}_{x}}\frac{\partial }{\partial r}\left( \frac{\partial r}{\partial x} \right) \\ {{u}_{r}}\left[ \frac{\partial cos\theta }{\partial r}-\frac{\partial }{\partial \theta }\left( \frac{sin\theta }{\partial y} \right) \right]cos\theta \\ {{u}_{r}}\left[ \frac{\partial sin\theta }{\partial r}+\frac{\partial }{\partial \theta }\left( \frac{cos\theta }{\partial z} \right) \right]sin\theta \\ \end{matrix} \right]}^{T}}$
将（21）式与（20）式联立，代入（7）式，可得柱坐标系下流体微元的压力-流速关系式：
$\left[ \begin{align} & \frac{\partial {{u}_{x}}}{\partial t} \\ & \frac{\partial {{u}_{r}}}{\partial t} \\ & \frac{\partial {{u}_{\theta }}}{\partial t} \\ \end{align} \right]=\frac{\mu }{\rho }\left[ \begin{matrix} \frac{{{\partial }^{2}}{{u}_{x}}}{\partial {{r}^{2}}}+\frac{{{\partial }^{2}}{{u}_{x}}}{\partial {{x}^{2}}}+\frac{\partial {{u}_{x}}}{r\partial r}+\frac{{{\partial }^{2}}{{u}_{x}}}{{{r}^{2}}\partial {{\theta }^{2}}}-\frac{{{\sigma }_{p}}}{\mu \partial x}+\frac{\rho gsin\alpha }{\mu } \\ \frac{{{\partial }^{2}}{{u}_{r}}}{\partial {{r}^{2}}}+\frac{{{\partial }^{2}}{{u}_{r}}}{\partial {{x}^{2}}}+\frac{\partial {{u}_{r}}}{r\partial r}+\frac{{{\partial }^{2}}{{u}_{\text{r}}}}{{{r}^{2}}\partial {{\theta }^{2}}}-\frac{2\partial {{u}_{\theta }}}{{{r}^{2}}\partial \theta }-\frac{{{u}_{\text{r}}}}{{{r}^{2}}}-\frac{{{\sigma }_{p}}}{\mu \partial r}-\frac{\rho gcos\alpha }{\mu } \\ \frac{{{\partial }^{2}}{{u}_{\theta }}}{\partial {{r}^{2}}}+\frac{{{\partial }^{2}}{{u}_{\theta }}}{\partial {{x}^{2}}}+\frac{\partial {{u}_{\theta }}}{r\partial r}+\frac{{{\partial }^{2}}{{u}_{\theta }}}{{{r}^{2}}\partial {{\theta }^{2}}}+\frac{2\partial {{u}_{r}}}{{{r}^{2}}\partial \theta }-\frac{{{\sigma }_{p}}}{\mu r\partial \theta } \\ \end{matrix} \right]$
结合流体的流动状态及前文的流态假设，流体微元不存在径向与切向的速度，进一步推导出柱坐标系下流体流动的连续性方程：
$\begin{matrix} ~\frac{\partial {{u}_{r}}}{\partial r}+\frac{{{u}_{r}}}{r}+\frac{\partial {{u}_{x}}}{\partial x}+\frac{\partial {{u}_{\theta }}}{r\partial \theta }=0 \\ \end{matrix}$
将（21）式、（23）式代入（22）式，根据圆管层流条件可知流体匀速稳定流动条件如下：
$\begin{matrix} ~\left\{ \begin{array}{*{35}{l}} {\partial p}/{\partial r}\;=\rho gcos\alpha \\ {\partial p}/{\partial \theta }\;=0 \\ \frac{{{\partial }^{2}}{{u}_{x}}}{\partial {{r}^{2}}}+\frac{\partial {{u}_{x}}}{r\partial r}=\frac{\partial p}{\mu \partial x}-\frac{\rho gsin\alpha }{\mu } \\ \end{array} \right. \\ \end{matrix}$
由（24）式可知，管路内流体连续流动工况下，流速与压力分别在径向、轴向上呈一维分布，求解管路内流体的流速分布得：
$\begin{matrix} ~u=\frac{1}{4\mu }\left( {{r}^{2}}-\frac{{{D}^{2}}}{4} \right)\left[ \frac{\partial \left( p-\rho grcos\alpha \right)}{\partial x}+\rho gsin\alpha \right] \\ \end{matrix}$
层流状态下，管路流体压力在流向上呈均匀递减分布，设单位长度井筒管路的流体压降为Δpe，进而求取管路的等效流量为：
$\begin{matrix} ~{{q}_{\text{v}}}=\mathop{\int }_{0}^{2\pi }d\theta \mathop{\int }_{0}^{\frac{D}{2}}urdr=\frac{\pi {{D}^{4}}}{128\mu } \\ \end{matrix}\left( \Delta {{p}_{e}}-\rho gsin\alpha \right)$
令连续管道的流导系数Z0=πD4/128μ，实际工程中流导系数与管壁材质粗糙程度有关：
$\begin{matrix} ~Z={{k}_{\varepsilon }}{{Z}_{0}} \\ \end{matrix}$
管路中流体的平均流速为：
$~\begin{matrix} \bar{u}=\frac{4{{q}_{\text{v}}}}{\pi {{D}^{2}}}=\frac{{{D}^{2}}}{32\mu }\left( \Delta {{p}_{e}}-\rho gsin\alpha \right)=\frac{{{u}_{max}}}{2} \\ \end{matrix}$
1.2.2 注水井筒空间流场分析
结合实际油田注水井井筒结构，建立井筒空间流场模型。设井深lx有一长度为dl、垂深为dh的微井段，对应垂深hx，水平角为α，全井井深为l，垂深为h，微井段入口压力为p1，出口压力为p2。根据（26）式可得微井段流量qv为：
$\begin{matrix} ~{{q}_{\text{v}}}=\frac{\pi {{D}^{4}}}{128\mu }\left( \frac{{{p}_{1}}-{{p}_{2}}}{\mathrm{d}l}-\rho g\mathrm{sin}\alpha \right) \\ \end{matrix}$
设全井井眼轨迹呈不规则弧状，则井段与水平线夹角为关于井深及垂深的函数$\varphi \left( {{l}_{x}},{{h}_{x}} \right)$，井口压力为p0，井底压力为ph，将压力损失进行全井段积分：
$~\begin{matrix} {{p}_{0}}-{{p}_{h}}=\mathop{\int }_{0}^{l}\left\{ \frac{128\mu {{q}_{\text{v}}}}{\pi {{D}^{4}}}-\rho gsin\left[ \varphi \left( {{l}_{x}},{{h}_{x}} \right) \right] \right\}dl \\ \end{matrix}$
当井深趋于无限小时，设入口处切线方向与水平线夹角为α1，井段出口处切线方向与水平线夹角为α2，微井段的空间结构如图3 所示。
由于dl极小，对应的井斜角变化量可视为无穷小，故可将微井段视为直段，即全井由无限多个长度为dl的井段组成，根据几何关系可得：
$\begin{matrix} ~h=\underset{{{l}_{x}}=0}{\overset{l}{\mathop \sum }}\,\varphi \left( {{l}_{x}},{{h}_{x}} \right)dl=\mathop{\int }_{0}^{l}sin\left[ \varphi \left( {{l}_{x}},{{h}_{x}} \right) \right]dl \\ \end{matrix}$
结合储集层模型计算全井筒的压力和流量特性。井底末端与储集层相通，可将井底压力ph视为地层压力，设定地层的吸水指数为k，开启压力为pb，根据层段吸水特性模型得：
$\begin{matrix} ~{{p}_{h}}={{k}^{-1}}{{q}_{\text{v}}}+{{p}_{b}}~ \\ \end{matrix}$
将（29）式代入（32）式，建立井筒流量与井筒结构、地层特性参数、井口压力的关系如（33）式所示：
$\begin{matrix} ~{{q}_{\text{v}}}=\frac{\pi {{D}^{4}}k\left( {{p}_{0}}+\rho gh-{{p}_{b}} \right)}{\pi {{D}^{4}}+128\mu kl} \\ \end{matrix}$
上式与（25）式和（28）式联立，可以获得井筒内流体流速分布：
$\begin{matrix} ~u=\frac{8k\left( {{p}_{0}}+\rho gh-{{p}_{b}} \right)}{\pi {{D}^{4}}+128\mu kl}\left( {{r}^{2}}-\frac{{{D}^{2}}}{4} \right)~ \\ \end{matrix}$
根据层流状态下压力梯度分布特性，联立（29）式和（33）式，得到距井口长度为lx处的压力px，即井筒动态压力空间分布：
$\begin{matrix} ~{{p}_{x}}=\left( 1-\frac{128\mu {{l}_{x}}k}{\pi {{D}^{4}}+128\mu kl} \right){{p}_{0}}+ \\ \end{matrix}\frac{128\mu {{l}_{x}}k\left( {{p}_{b}}-\rho gh \right)}{\pi {{D}^{4}}+128\mu kl}+\rho g{{h}_{x}}$
1.3 注水井井筒流体波传输时域特性
井筒流体波生成的时域延迟特性是设计无线通信编解码的重要依据。流体波由节流元件开度周期性变化或泵压规律变化时激发的水击效应产生（见图4 ）：井筒内某处流速改变，由于流段存在惯性，会对设备产生冲击波，并沿井筒轴向传播。设定流段体积为ΔV，A为管道截面积，管道内径为D，流体波传播速度为c，管壁厚度为e。此段流体由于阀门开度减小而被阻挡和压缩，压缩后总压力增加、密度发生变化，流速由u0变为u1，根据动量原理可得：
根据材料力学理论，脉冲生成的物理过程与油管材料、流体物性有关。根据流体流动的连续性与质量守恒关系可计算流体弹性模量E（见（37）式）；管材的弹性模量E0为流体压缩后径向应力与内径变化之比，设管材的应力变化为dσ，推导得管材弹性模量与应变关系如（38）式所示。
$\begin{matrix} ~E=V\frac{\text{d}p}{dV}=\rho \frac{dp}{d\rho }~ \\ \end{matrix}$
$\begin{matrix} ~{{E}_{0}}=\frac{d\sigma }{{dD}/{D}\;}=\frac{{{D}^{2}}dp}{2edD} \\ \end{matrix}$
设定在时间Δt内，水击波传播距离长度为Δls（Δls=cΔt）。井筒远端未受到水击压力作用的流体由于惯性作用，仍以初始速度u0流入该流段，以速度u1流出，导致流段内质量变化。受管材弹性作用，该流段内管道截面积由A增至A+dA，流体密度由ρ增至ρ+dρ，则该流段的质量增值为：
$\left( \rho +d\rho \right)\left( A+dA \right)c\Delta t-\rho Ac\Delta t=\rho A\left( {{u}_{1}}-{{u}_{0}} \right)\Delta t$
化简得到速度、面积与密度变化的关系式：
$\begin{matrix} ~\frac{{{u}_{1}}-{{u}_{0}}}{c}=\frac{dA}{A}+\frac{d\rho }{\rho } \\ \end{matrix}$
由管道内径变化与截面积关系可得：
$\begin{matrix} ~\frac{dA}{A}=\frac{d{{D}^{2}}}{{{D}^{2}}}=\frac{2dD}{D}=\frac{Ddp}{e{{E}_{0}}} \\ \end{matrix}$
上式与（36）式联立，化简后得：
$\begin{matrix} ~c=\sqrt{\frac{eE{{E}_{0}}}{\rho \left( e{{E}_{0}}+DE \right)}\left( 1-\frac{1}{2g}\frac{\partial {{u}^{2}}}{\partial {{l}_{x}}} \right)} \\ \end{matrix}$
由注水层流条件可知，流体流速与水击波速相比变化幅度极小，且井筒压力变化幅度较大，因此水击波传输速度可等价为：
$\begin{matrix} ~c=\sqrt{\frac{eE{{E}_{0}}}{\rho \left( e{{E}_{0}}+DE \right)}}\quad \\ \end{matrix}$ $\left( \frac{\partial {{u}^{2}}}{2g\partial {{l}_{x}}}\approx 0 \right)$
部分水击信号能量被管壁摩擦阻力消耗，从而产生阻尼衰减。设定节流元件处产生的压力波幅值为Δpl，lx处的压力波幅值为Δpx，定义二者之比为脉冲传输信号衰减系数δx，则根据指数衰减规律可得：
$\begin{matrix} ~{{\delta }_{x}}=\frac{\Delta {{p}_{x}}}{\Delta {{p}_{l}}}=exp\left( -\frac{{{l}_{x}}}{{{l}_{\delta ={{\text{e}}^{-1}}}}} \right) \\ \end{matrix}$
其中，${{l}_{\delta ={{\text{e}}^{-1}}}}$为信号衰减到原信号强度Δpl的e−1倍时的传输距离，与信号生成频率f有关，推导可得：
$\begin{matrix} ~{{l}_{\delta ={{\text{e}}^{-1}}}}=\frac{cD}{2}\sqrt{\frac{\rho }{\pi f\mu }}=D\sqrt{\frac{eE{{E}_{0}}T}{2\pi \mu \left( e{{E}_{0}}+DE \right)}}~ \\ \end{matrix}$
脉冲信号频率f为节流元件开关的速度，一次脉冲为节流元件一组开度往复过程，脉冲时间（T）须大于等于传输延迟时间2l/c，且有T=1/2f，可得：
$\begin{matrix} ~f\le \frac{c}{4l} \\ \end{matrix}$
经信号衰减后lx处的压力波幅值Δpx为：
$\begin{matrix} \Delta {{p}_{x}}=\Delta {{p}_{l}}{{\delta }_{x}}=\Delta {{p}_{l}}exp\left[ -\frac{{{l}_{x}}}{D}\sqrt{\frac{2\pi \mu \left( e{{E}_{0}}+DE \right)}{\rho eE{{E}_{0}}T}} \right] \\ \end{matrix}$
2 仿真实验与结果分析
选取中国某油田一口典型注水井进行仿真分析，该井井深和垂深分别为1 600，1 400 m，为斜井结构，储集层开启压力1.5 MPa，吸水特性指数25.2 m3/（d·MPa）。此工况适用波码通信式分层注水工艺技术。
2.1 恒定注入流量工况注水井井筒压力动态分布
利用本文提出的流体微元压力-流速关系模型，结合实际生产工况，进行井筒压力动态分布模拟。设定全井注水量为80 m3/d，井口压力分别为6，8，10 MPa，模拟计算全井筒压力沿井深方向和垂深方向的动态分布，结果如图5 所示。
2.2 恒定注入压力工况注水井井筒压力动态分布
设定井筒注入压力为8 MPa，注入量分别为60，80，100 m3/d，进行井筒压力动态分布数值模拟。模拟结果如图6 所示，虚线为停注状态下忽略沿程阻力损失的压力动态分布，与传统伯努利方程模型求解结果一致，验证了模型的准确性。井口注入量为60，80，100 m3/d对应的井底压力分别为21.72，21.29，20.87 MPa（见图6a ），与静态压力值分别相差0.18，0.61，1.03 MPa；井筒压力沿垂深方向呈近似线性分布（见图6b ）。对比分析可知，注入压力恒定的工况下，井口注入流量越大，井筒沿程阻力损失越大。利用该模型可定量分析井筒通信过程中地面压力与井下压力绝对值的对应关系。
2.3 注水井井筒波码传输特性
井筒波码通信编解码方案设计中两个核心参数为信号延时和信号幅值。鉴于此，利用本文水击波传输速度模型研究不同油管尺寸下流体波信号的延时规律，选取3种常规油管对流体波信号的动态响应特性进行仿真计算。内径为50.3，62.0，75.9 mm的常规油管对应井深1 000 m处的传输延迟分别为1.473，1.470，1.462 s，对应井深1 600 m处的传输延迟分别达2.361，2.356，2.343 s（见图7 ）。对比分析可知，内径越小延时越长，但其影响为毫秒级，可以忽略。因此，延时幅值主要取决于井筒的长度，即针对较长的井筒，设计井筒波码通信编解码方案需要考虑延时的影响：指令下达与数据回传的时间间隔应大于延时时间。需进一步在考虑通信效率和误码率的条件下，开展该环节工程量化关系研究。
研究不同内径油管内流体波信号传输衰减规律，可有效指导井筒波码无线通信压力识别幅值的设计。识别幅值过大会导致无法检测，过小会造成硬件资源过度损耗。内径分别为50.3，62.0，75.9 mm的油管对应井深1 600 m处的信号衰减幅度分别为5.50%，5.05%，4.68%（见图8 ）。对比分析可知，管径越小，流体波信号衰减幅度越大。
现场应用中，压力检测通常具有一定误差，需要分析压力幅值动态变化以实现井筒无线信号传输。针对现有工程及监测条件下的压差识别有效值，大幅值压差更有利于井筒无线通信；对于井口升压和低渗透油藏层段压差变化小的情况，工程实现难度大，需要针对不同工况量化压力变化，根据不同的地层特性计算压差识别有效值（即传感器的压差阈值），进而提高通信效率。本文模拟实验中给定了压差识别有效值，利用模型开展了不同井口压力波幅值下井筒压力动态分布模拟，定量分析适用条件（不考虑层间干扰）。
仿真结果显示，井深越大、井口压力波幅值越大，压力波信号衰减越大，井口压力波幅值为3 MPa时，井底压力波幅值为2.849 MPa，衰减0.151 MPa，因此，需要通过增加井口压力变化幅值和压力识别精度提高通信可靠性。图9 模拟计算了不同井深处压差识别有效值分别为0.5，1.0，1.5 MPa时，井口所需的最小压力波幅值，特征曲线的右侧区域为在井下可识别压差与井深约束下满足流体波无线通信的井口压力波幅值范围，该特征曲线表征了井下压差识别有效值、井深与井口所需最小压力波幅值的关系，为后续不同井深井况下井筒波码高效无线通信编解码方案的制定提供理论数据支撑。
3 结论
基于流体力学微元分析法，建立流体微元压力-流速关系模型并拓展至全井空间流场，结合井筒流体波传输的时域特性，得到井筒空间的流场流速、压力的动态时空分布模型。基于某典型注水井生产数据进行数值模拟分析，发现井筒沿程压力损失与绝对压力值无关，设计波码通信编解码方案时无须考虑注入的绝对压力值；注入压力恒定时，井口注入流量越大，井筒沿程损失越大。流体波信号延时幅值主要取决于井筒的长度；管径越小，流体波信号衰减幅度越大。同井深下生成目标波码幅值（压差识别有效值）越高，为克服井筒沿程损失所需的井口压力波幅值越大。该模型可为分层注水工艺波码通信效率的提高优化提供机理模型和理论数据支撑。
符号注释：
a——流体微元运动加速度，m/s2；a——微元六面体的加速度矩阵，m/s2；a0——流体微元法向应力项系数，无因次；b0——流体微元本构关系项系数，无因次；A——管道截面积，m2；c——流体波传播速度，m/s；D——管道内径，m；e——井筒管道壁厚，m；eσ——流体微元运动状态下的法向应力矩阵，Pa；$\overline{{{e}_{\sigma }}}$——流体微元静止状态下的法向应力矩阵，Pa；eτ——流体微元的本构切向应力矩阵，Pa；eστ——流体微元所受合应力，Pa；E——流体弹性模量，Pa；E0——管材弹性模量，Pa；f——流体波生成频率，Hz；F——流体微元所受合力，N；g——重力加速度，m/s2；G——流体微元所受重力，N；G——流体微元的重力矩阵，N；h——全井垂深，m；hx——井筒空间点对应垂深，m；I——单位矩阵；k——地层吸水指数，m3/（s·Pa）；kε——材质系数，无因次；l——全井井深，m；lx——井筒空间点对应井深，m；Δls——水击效应的流段长度，m；m——微元质量，kg；p——流体压力，Pa；pb——地层开启压力，Pa；ph——井底压力，Pa；px——距井口长度为lx处的压力，Pa；p0——井口压力，Pa；p1——微井段入口压力，Pa；p2——微井段出口压力，Pa；qv——体积流量，m3/s；Δp——压力波动，Pa；Δpe——单位长度井筒管路流体压降，Pa/m；Δpx——lx处的压力波幅值，Pa；Δpl——节流元件处的压力波幅值，Pa；r——柱坐标系中的径向距离，m；Re——雷诺数，无因次；t——时间，s；Δt——时间段，s；T——信号脉冲周期，s；u——流体流速，m/s；u——微元六面体的速度矩阵，m/s；$\bar{u}$——平均流速，m/s；u0，u1——流体压缩前后的流速，m/s；V——流体体积，m3；ΔV——流段体积，m3；x，y，z——直角坐标系，m；Z——实际工程中流导系数，m4/（s·Pa）；Z0——连续管道流导系数，m4/（s·Pa）；α——流体微元与水平线夹角，（°）；α1——井段入口处切线与水平线夹角，（°）；α2——井段出口处切线与水平方向夹角，（°）；δx——信号衰减系数，无因次；θ——柱坐标系中的极角，（°）；μ——流体黏度，Pa·s；ρ——流体密度，kg/m3；σ——流体微元所受法向应力，Pa；σp——流体微元的法向应力矩阵，Pa；σp——静止状态下流体微元各受力面所受法向应力，Pa；σxb，σyb，σzb——x、y、z轴向后方微元产生的垂直于受力面的法向压应力，Pa；σxf，σyf，σzf——x、y、z轴向前方微元产生的垂直于受力面的法向压应力，Pa；τ——流体微元所受切向应力，Pa；τq——流体微元的切向应力矩阵，Pa；τxy，τxz——x方向上形变对y、z方向产生的切向应力，Pa；τyx，τyz——y方向上形变对x、z方向产生的切向应力，Pa；τzx，τzy——z方向上形变对x、y方向产生的切向应力，Pa；φ——井段与水平线夹角，为井深及垂深的函数，（°）。下标：u，l——与轴向垂直方向上的微元顶面和底面；max——最大值。