输流管道系统广泛应用于工业生产生活的各个领域，其流量调节、阀门快关或意外事故如断电停泵等均有可能触发水击即水锤现象，导致流体压强大幅波动，危害管道系统的安全性。根据是否考虑流固耦合(fluid-structure interaction, FSI)效应，求解模型可以分为经典水击理论^[1]和流固耦合水击理论。

忽略热量传递，管道系统的流固耦合效应是指在非定常流动过程中，管道与所含流体之间动量和力的双向传递。流固耦合方式主要有Poisson耦合、摩擦耦合和结合部耦合3种。Poisson耦合由管壁径向胀缩产生^[2]，耦合程度与管道的呼吸模态和Poisson比有关^[3]；摩擦耦合存在于黏性流体与管道内壁之间，由二者的相对运动产生^[4-7]；结合部耦合主要发生在截面或边界条件突变的管道连接处如阀门、弯管、分支管和支撑等，是一种强烈的耦合方式。张挺等^[8]改变阀门处边界条件，发现结合部耦合会影响系统特征参数的振动幅值和频率。Riedelmeier等^[9-10]建立了弯头系统实验台架并进行了水击特性研究，通过频域方法进行了评价，结果表明结合部耦合在长管道和中管道弯头系统中较为显著。Ahmadi等^[11]采用特征线法和有限元法，对分支管道系统和泵系统中的结合部耦合效应进行了深入研究，指出结合部耦合高度依赖于管道系统的刚度，且考虑结合部耦合后的压强会略高于经典水击理论的压强预测值。Ferras等^[12]研究了支撑处的结合部耦合对瞬变压强波传播的影响。Keramat等^[13]开展了水击实验，分析了压强信号在时域和频域的表现，研究表明支撑的位置和刚度对系统瞬态特性影响较大。Hosseini^[14]考虑了Poisson耦合和结合部耦合，采用特征线法研究了黏弹性支撑对水击现象的影响。对于结合部耦合，以上研究大多设置某一监测点，通过实验测量或理论模型求解得到监测点处的表征参数如流体的压强、管道的轴向速度和位移等，进一步对表征参数展开时域和频域的分析。但该类研究方法只能分析系统某一点处物理量的变化，未能解释系统在瞬变过程中的整体响应，难以定量化描述流固耦合的性质，因此本文基于能量分析方法，引入了多种流体能量和管道能量，从系统整体层面研究分析结合部耦合对水击瞬态响应特性的影响。

针对经典水击理论，Karney^[15]引入了能量的概念，建立了最早的能量分析方法。之后Kung^[16]考虑了对流项，推广了能量方程，对系统能量转换的机制进行了更充分的讨论。基于能量分析方法，Ghidaoui等^[17]量化了波速调整法、时间插值和空间插值方案的离散误差，同时指出了减小离散误差的方法。Meniconi等^[18]将能量分析方法应用在反瞬态分析中，指出同一衰减曲线可能对应不同的瞬态扰动条件。朱炎等^[19-20]针对气液两相瞬变流过程，研究了不同含气率对瞬变流能量衰减的影响。Pan等^[21]分析了弹性管道系统和黏弹性管道系统的能量转换和耗散机制，研究了拟稳态摩阻、非恒定摩阻和管壁黏弹性的影响。但上述研究仅适用于经典水击理论，对于流固耦合水击理论不适用。因此本文针对流固耦合模型，提出了流固耦合能量分析方法。

本文建立了基于流固耦合水击理论的能量分析方法。首先采用特征线法对流固耦合4方程模型进行了数值求解，通过基准算例验证了理论模型。之后对模型进行数学推导，获得了管道和流体能量的表达式及对应的能量方程。最后分析了能量在流体和管道之间的传递和转换过程，通过无量纲因子描述了结合部耦合对水击瞬变响应特性的影响。

1 流固耦合4方程模型

对于线弹性的薄壁管道，忽略管壁的径向惯性和对流项，同时假定瞬变流动过程中不发生空化现象，推导得到流固耦合4方程模型^[22]。该模型重点研究流体的轴向运动和管道的轴向运动，包含流体的连续性方程、流体的动量方程、管道的动量方程和本构关系^[23]：

$ \frac{\partial V}{\partial t}+g \frac{\partial H}{\partial z}=-f \frac{V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}, $

(1)

$ \frac{\partial V}{\partial z}+\frac{g}{c_{\rm{f}}^2} \frac{\partial H}{\partial t}=2 \nu \frac{\partial U}{\partial z}, $

(2)

$ \frac{\partial U}{\partial t}-\frac{1}{\rho_{\mathrm{t}}} \frac{\partial \sigma}{\partial z}=f \frac{\rho_{\mathrm{f}} A_{\mathrm{f}}}{\rho_{\mathrm{t}} A_{\mathrm{t}}} \frac{V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}+g \sin \alpha, $

(3)

$ \frac{\partial U}{\partial z}-\frac{1}{\rho_{\mathrm{t}} c_{\mathrm{t}}^2} \frac{\partial \sigma}{\partial t}=-\rho_{\mathrm{f}} g \frac{\nu R}{E e} \frac{\partial H}{\partial t} . $

(4)

其中：下标f代表流体参数，下标t代表管道参数，下标r代表流体物理量与管道物理量之差；g为重力加速度，本文取9.8 m/s²；V和U分别为流体和管道的轴向运动速度，V_r为两者速度差，m/s；H为流体压强水头，m；σ为管道的轴向应力，Pa；ρ_f和ρ_t分别为流体和管道材料的密度，kg/m³；A_f和A_t分别为流体和管道截面的面积，m²；K为流体体积模量，Pa；E为管道弹性模量，Pa；z为轴向长度，m；t为时间，s；R为管道内半径，e为管壁厚度，m；f为摩擦阻力因数；ν为Poisson比；α为管道轴线与水平面的夹角，以管道提升高度降低为正方向。方程中压强波速c_f和应力波速c_t计算如下：

$c_{\mathrm{f}}^2=\left[\rho_{\mathrm{f}}\left(\frac{1}{K}+\left(1-\nu^2\right) \frac{2 R}{E e}\right)\right]^{-1}, $

(5)

$ c_{\mathrm{t}}^2=\frac{E}{\rho_{\mathrm{t}}} . $

(6)

1.1 数值求解

将流固耦合4方程模型写成矩阵形式：

$ \boldsymbol{A} \frac{\partial}{\partial t} \boldsymbol{\varphi}(z, t)+\boldsymbol{B} \frac{\partial}{\partial z} \boldsymbol{\varphi}(z, t)=\boldsymbol{r}. $

(7)

其中：A和B是系数矩阵，φ(z, t)=(V, H, U, σ)^T，r为1×4维的向量。

分析矩阵系数可得，特征方程|B-λA|有4个不同的实根，等价于矩阵A^-1 B有4个不同的实特征值：

$ \lambda_1=\widetilde{c}_{\mathrm{f}}=\left\{\frac{1}{2}\left[\gamma^2-\left(\gamma^4-4 c_{\mathrm{f}}^2 c_{\mathrm{t}}^2\right)^{\frac{1}{2}}\right]\right\}^{\frac{1}{2}}, $

(8)

$ \lambda_2=-\widetilde{c}_{\mathrm{f}}=-\left\{\frac{1}{2}\left[\gamma^2-\left(\gamma^4-4 c_{\mathrm{f}}^2 c_{\mathrm{t}}^2\right)^{\frac{1}{2}}\right]\right\}^{\frac{1}{2}}, $

(9)

$ \lambda_3=\widetilde{c}_{\mathrm{t}}=\left\{\frac{1}{2}\left[\gamma^2+\left(\gamma^4-4 c_{\mathrm{f}}^2 c_{\mathrm{t}}^2\right)^{\frac{1}{2}}\right]\right\}^{\frac{1}{2}}, $

(10)

$ \lambda_4=-\widetilde{c}_{\mathrm{t}}=-\left\{\frac{1}{2}\left[\gamma^2+\left(\gamma^4-4 c_{\mathrm{f}}^2 c_{\mathrm{t}}^2\right)^{\frac{1}{2}}\right]\right\}^{\frac{1}{2}} . $

(11)

其中：

$ \gamma^2=\left(1+2 \nu^2 \frac{\rho_{\mathrm{f}}}{\rho_{\mathrm{t}}} \frac{R}{e}\right) c_{\mathrm{f}}^2+c_{\mathrm{t}}^2. $

(12)

使用特征线法，结合4个特征值，将式(7)转换为4个常微分方程进行求解。在图 1中的特征线A₁、A₂、A₃和A₄上有相容方程式(13)—(16)成立。

$ {\xi _{{\rm{f}}\mathit{H}}}\frac{{{\rm{d}}H}}{{{\rm{d}}t}} + {\xi _{{\rm{f}}\mathit{V}}}\frac{{{\rm{d}}V}}{{{\rm{d}}t}} - {\xi _{{\rm{f}}\sigma }}\frac{{{\rm{d}}\sigma }}{{{\rm{d}}t}} + {\xi _{{\rm{f}}\mathit{U}}}\frac{{{\rm{d}}U}}{{{\rm{d}}t}} = {q_{\rm{f}}}, $

(13)

$ - {\xi _{{\rm{f}}\mathit{H}}}\frac{{{\rm{d}}H}}{{{\rm{d}}t}} + {\xi _{{\rm{F}}\mathit{V}}}\frac{{{\rm{d}}V}}{{{\rm{d}}t}} + {\xi _{{\rm{f}}\sigma }}\frac{{{\rm{d}}\sigma }}{{{\rm{d}}t}} + {\xi _{{\rm{f}}U}}\frac{{{\rm{d}}U}}{{{\rm{d}}t}} = {q_{\rm{f}}}, $

(14)

$ - {\xi _{{\rm{t}}H}}\frac{{{\rm{d}}H}}{{{\rm{d}}t}} - {\xi _{{\rm{t}}\mathit{V}}}\frac{{{\rm{d}}V}}{{{\rm{d}}t}} - {\xi _{t\sigma }}\frac{{{\rm{d}}\sigma }}{{{\rm{d}}t}} + {\xi _{{\rm{t}}\mathit{U}}}\frac{{{\rm{d}}U}}{{{\rm{d}}t}} = {q_{\rm{t}}}, $

(15)

$ {\xi _{{\rm{t}}\mathit{H}}}\frac{{{\rm{d}}H}}{{{\rm{d}}t}} - {\xi _{{\rm{t}}\mathit{V}}}\frac{{{\rm{d}}V}}{{{\rm{d}}t}} + {\xi _{{\rm{t}}\sigma }}\frac{{{\rm{d}}\sigma }}{{{\rm{d}}t}} + {\xi _{{\rm{t}}U}}\frac{{{\rm{d}}U}}{{{\rm{d}}t}} = {q_{\rm{t}}}. $

(16)

其中下标H、V、σ和U代表对应的物理量。

$ \begin{gathered} \xi_{\mathrm{f}H}=\frac{g}{c_{\mathrm{f}}}\left[\left(\frac{\widetilde{c}_{\mathrm{f}}}{c_{\mathrm{f}}}\right)+\frac{2 \nu^2 R \rho_{\mathrm{f}}}{e \rho_{\mathrm{t}}} \frac{\left(c_{\mathrm{f}} \widetilde{c}_{\mathrm{f}}\right) / c_{\mathrm{t}}^2}{1-\left(\widetilde{c}_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}\right], \\ \xi_{\mathrm{f}V}=1, \\ \xi_{\mathrm{f}\sigma }=\frac{2 \nu}{\rho_{\mathrm{t}} \widetilde{c}_{\mathrm{f}}} \frac{\left(\widetilde{c}_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}{1-\left(\widetilde{c}_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}, \\ \xi_{\mathrm{f}U}=2 \nu \frac{\left(\widetilde{c}_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}{1-\left({\left.\widetilde{c}_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}\right.}, \end{gathered} $

$ \begin{gathered} q_{\mathrm{f}}=\xi_{\mathrm{f}U} g \sin \alpha+\left(\frac{\rho_{\mathrm{f}} A_{\mathrm{f}}}{\rho_{\mathrm{t}} A_{\mathrm{t}}} \xi_{\mathrm{f}U}-\xi_{\mathrm{f}V}\right) \frac{f V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}, \\ \xi_{\mathrm{t} H}=\frac{\nu g R \rho_{\mathrm{f}}}{e c_{\mathrm{t}} \rho_{\mathrm{t}}}\left(\frac{\widetilde{c}_{\mathrm{t}}}{c_{\mathrm{t}}}\right) \frac{\left(c_{\mathrm{f}} / \widetilde{c}_{\mathrm{t}}\right)^2}{1-\left(c_{\mathrm{f}} / \widetilde{c}_{\mathrm{t}}\right)^2}, \\ \xi_{\mathrm{t} V}=\frac{\nu R \rho_{\mathrm{f}}}{e \rho_{\mathrm{t}}} \frac{\left(c_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}{1-\left(c_{\mathrm{f}} / \widetilde{c}_{\mathrm{t}}\right)^2}, \\ \xi_{\mathrm{t} \sigma}=\frac{1}{\rho_{\mathrm{t}} c_{\mathrm{t}}}\left(\frac{\widetilde{c}_{\mathrm{t}}}{c_{\mathrm{t}}}\right), \end{gathered} $

$ \begin{gathered} \xi_{\mathrm{t}U}=1+\frac{2 \nu^2 R \rho_{\mathrm{f}}}{e \rho_{\mathrm{t}}} \frac{\left(c_{\mathrm{f}} / c_{\mathrm{t}}\right)^2}{1-\left(c_{\mathrm{f}} / \widetilde{c}_{\mathrm{t}}\right)^2}, \\ q_{\mathrm{t}}=\xi_{\mathrm{t}U} g \sin \alpha+\left(\frac{\rho_{\mathrm{f}} A_{\mathrm{f}}}{\rho_{\mathrm{t}} A_{\mathrm{t}}} \xi_{\mathrm{t}U}+\xi_{\mathrm{t}V}\right) \frac{f V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R} . \end{gathered} $

对式(13)—(16)沿着对应的特征线积分并联立求解，即可获得下一时刻内部节点的物理量。

1.2 边界条件

本文研究水箱、管道和阀门组成的管道系统，上游水箱和下游阀门处有4个未知物理量，但只有2条特征线，因此需要补充2个辅助方程。

1.2.1 水箱

水箱处的压强水头H_res保持恒定，同时连接处的管道刚性固定在地面上，不发生运动，补充式(17)和(18)为辅助方程。

$ H=H_{\mathrm{res}}, $

(17)

$ U=0. $

(18)

1.2.2 阀门

阀门处的边界条件有固支和简支2种。

在阀门固支工况中，瞬时关闭的阀门完全固定在地面上，补充式(19)和(20)为辅助方程：

$ V=0, $

(19)

$ U=0. $

(20)

在阀门简支工况中，瞬间关闭的阀门在轴向上能自由移动，补充式(21)和(22)为辅助方程：

$ V=U, $

(21)

$ A_{\mathrm{t}} \sigma=\rho_{\mathrm{f}} g A_{\mathrm{f}} \Delta H. $

(22)

其中ΔH为阀门两侧的压强水头差。

1.3 模型验证

代尔夫特基准测试A(Deft Benchmark problem A)^[22]是校核流固耦合4方程模型的基准算例。系统结构布置如图 2所示，管道长度为20 m，内径为0.398 5 m，壁厚为0.008 m，弹性模量为210 GPa，Poisson比为0.3，管壁材料的密度为7 900 kg/m³，流体密度为1 000 kg/m³，流体体积模量为2.1 GPa，摩擦阻力因数为0.02，初始流体速度为1 m/s。阀门瞬时关闭，边界条件包含固支和简支两种，阀门后侧压强为0 Pa。

在阀门固支和简支的工况下，对比模型仿真值与基准算例值，二者符合很好(见图 3)，验证了流固耦合4方程模型数值求解和编程实现的正确性。

2 流固耦合能量分析方法

对流固耦合模型进行推导和变换，获得了管道和流体的能量方程，叠加得到了总能量方程，最终建立了流固耦合能量分析方法，获得了系统在瞬态过程中的整体响应，有助于理解水击瞬变过程，评价水击瞬变响应特性。

2.1 管道能量方程

将式(3)与U相乘，得到

$ \frac{1}{2} \frac{\partial U^2}{\partial t}-\frac{U}{\rho_{\mathrm{t}}} \frac{\partial \sigma}{\partial z}=U g \sin \alpha+f \frac{\rho_{\mathrm{f}} A_{\mathrm{f}}}{\rho_{\mathrm{t}} A_{\mathrm{t}}} \frac{U V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}. $

(23)

将式(4)与-σ/ρ_t相乘，得到

$ -\frac{\sigma}{\rho_{\mathrm{t}}} \frac{\partial U}{\partial z}+\frac{\sigma}{\rho_{\mathrm{t}}^2 c_{\mathrm{t}}^2} \frac{\partial \sigma}{\partial t}=\frac{\sigma}{\rho_{\mathrm{t}}} \frac{\rho_{\mathrm{f}} g \nu R}{E e} \frac{\partial H}{\partial t} . $

(24)

将式(23)和(24)相加，得到

$ \begin{gathered} \frac{1}{2} \frac{\partial U^2}{\partial t}+\frac{\sigma}{\rho_{\mathrm{t}}^2 c_{\mathrm{t}}^2} \frac{\partial \sigma}{\partial t}=U g \sin \alpha+\frac{1}{\rho_{\mathrm{t}}} \frac{\partial(\sigma U)}{\partial z}+ \\ f \frac{\rho_{\mathrm{f}} A_{\mathrm{f}}}{\rho_{\mathrm{t}} A_{\mathrm{t}}} \frac{U V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}+\frac{\sigma}{\rho_{\mathrm{t}}} \frac{\rho_{\mathrm{f}} g \nu R}{E e} \frac{\partial H}{\partial t} . \end{gathered} $

(25)

对于截面均匀的单管，将式(25)与ρ_tA_t相乘，并且在管道轴向z上积分，得到

$ \frac{\partial K_{\mathrm{t}}}{\partial t}+\frac{\partial U_{\mathrm{t}}}{\partial t}=G_{\mathrm{t}}+W_{\mathrm{t}}+F_{\mathrm{t}}+P_{\mathrm{t}}. $

(26)

其中变量表达式和含义如表 1所示。

式(26)与热力学第一定律相对应，因此称为能量方程。等号左侧代表管道能量，包括动能和应变能；右侧代表做功项，包括重力做功、外部系统输入功、流体对管道施加的摩擦力做功和Poisson耦合效应做功。

2.2 流体能量方程

将式(1)与V相乘，得到

$ \frac{1}{2} \frac{\partial V^2}{\partial t}+V g \frac{\partial H}{\partial z}+f \frac{V V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}=0. $

(27)

将式(2)与gH相乘，得到

$ g H \frac{\partial V}{\partial z}+\frac{g^2}{2 c_{\mathrm{f}}^2} \frac{\partial H^2}{\partial t}=2 \nu g H \frac{\partial U}{\partial z} . $

(28)

将式(27)和式(28)相加，得到

$ \begin{gathered} \frac{g^2}{2 c_{\mathrm{f}}^2} \frac{\partial H^2}{\partial t}+\frac{1}{2} \frac{\partial V^2}{\partial t}=-g \frac{\partial(H V)}{\partial z}- \\ f \frac{V V_{\mathrm{r}}\left|V_{\mathrm{r}}\right|}{4 R}+2 \nu g H \frac{\partial U}{\partial z} . \end{gathered} $

(29)

对于截面均匀的单管，将式(29)与ρ_fA_f相乘，并且在管道轴向z上积分，得到

$ \frac{\partial I_{\mathrm{f}}}{\partial t}+\frac{\partial K_{\mathrm{f}}}{\partial t}=W_{\mathrm{f}}+F_{\mathrm{f}}+P_{\mathrm{f}}. $

(30)

其中变量表达式和含义如表 2所示。

等号左侧是流体能量，包括动能和内能；右侧是做功项，包括外部系统输入功、管道对流体施加的摩擦力做功和Poisson耦合效应做功。

2.3 总能量方程

将管道能量方程与流体能量方程相加，得到总能量方程：

$ \begin{gathered} \left(\frac{\partial K_{\mathrm{t}}}{\partial t}+\frac{\partial U_{\mathrm{t}}}{\partial t}\right)+\left(\frac{\partial I_{\mathrm{f}}}{\partial t}+\frac{\partial K_{\mathrm{f}}}{\partial t}\right)= \\ G_{\mathrm{t}}+\left(W_{\mathrm{t}}+W_{\mathrm{f}}\right)+\left(F_{\mathrm{t}}+F_{\mathrm{f}}\right)+\left(P_{\mathrm{t}}+P_{\mathrm{f}}\right). \end{gathered} $

(31)

定义管道总能量E_t为动能K_t和应变能U_t之和，流体总能量E_f为动能K_f和内能I_f之和。同时考虑到流体速度远小于波速，忽略对流项的影响，偏导数近似为全导数，总能量方程为

$ \begin{gathered} \frac{\mathrm{d} E_{\mathrm{t}}}{\mathrm{d} t}+\frac{\mathrm{d} E_{\mathrm{f}}}{\mathrm{d} t}= \\ G_{\mathrm{t}}+\left(W_{\mathrm{t}}+W_{\mathrm{f}}\right)+\left(F_{\mathrm{t}}+F_{\mathrm{f}}\right)+\left(P_{\mathrm{t}}+P_{\mathrm{f}}\right) . \end{gathered} $

(32)

摩擦耦合和Poisson耦合效应直接作用在整根管道上，所以在流固耦合模型有对应的表达式。结合部耦合只作用在管道边界条件突变处，通过边界条件对水击瞬变流动过程施加影响，不直接体现在模型中。图 4展示了流体能量、管道能量与流固耦合效应的关系。

3 能量分析结果

将能量分析作为流固耦合水击仿真的后处理部分，通过对节点的物理量进行数值积分后，分析能量在流体和管道之间的传递和转换过程。在阀门固支和简支工况中，流体能量和管道能量变化趋势相似。以下以阀门固支工况为例，采用代尔夫特基准测试A的管路布置和参数，为了贴合实际工程管道系统的运行工况，仅修改阀门后的压强为2 MPa，使整个水击过程中流体压强最小值始终大于0。在其他工况中，流固耦合水击计算和能量分析过程与此相似。

3.1 流体内能和动能

能量分析方法清晰地展现了关阀水击过程中流体动能和内能的转换关系，如图 5所示。根据能量的变化趋势，一个水击周期可以划分为4个阶段，关键的时间节点为0、t₁、t₂、t₃和t₄。在第1阶段(0 < t≤t₁)中，阀门瞬间关闭后，附近的流体首先受到压缩而变为静止状态，形成一个从下游阀门向上游水箱传播的升压压强波，这个过程中，流体动能逐渐减小，内能逐渐增加。在t₁时刻，管道内的流体全部静止，流体动能为0，内能达到最大值。随后在第2阶段(t₁ < t≤t₂)中，由于管道内的压强高于上游压强，整个系统处于一种不平衡状态，流体要倒流回上游水箱，这个过程中流体内能逐渐减小，动能逐渐增加，于t₂时刻达到最大值。因为下游阀门处于关闭状态，在阀门处不可能有流体来一直维持这个流动，因此在第3阶段(t₂ < t≤t₃)中，阀门附近出现了一个低压区域，且低压区域向上游不断延伸，使流体出流的速度减小。在这一阶段中，流体内能逐渐减小，动能也逐渐减小。在t₃时刻，流体再次处于静止状态，动能和内能均处于最小值。这样又形成了上游边界处不平衡的流动条件，流体从上游水箱流入管道内，动能和内能都逐渐增大，直至t₄时刻进入下一个水击周期。基于上述分析，能量分析方法提供了自然和直接的角度来理解系统的动态响应过程，易于描述水力系统中的能量传递过程，可视化性强。

3.2 流体能量和管道能量

图 6展示了流体能量和管道能量随响应时间的变化情况，各能量按照平均值从大到小排序，依次为I_f、K_f、U_t和K_t。K_t仅在0~16 J变化，远小于其他能量，可以忽略。因此流体能量中占主导地位的是I_f，管道能量中占主导地位的是U_t。

图 6还表明I_f与U_t的周期相等，均为K_f周期的2倍。以管道中点处物理量为例进行说明，流体速度、压强水头和管道应力的变化周期相同，变化规律如图 7所示。其中流体压强水头和管道应力始终大于0，而流体速度在0位置附近上下波动，这一规律对所有节点均成立。因此流体内能和管道应变能在1个周期中只达到1次最大值，动能达到2次最大值，即流体内能与管道应变能的周期相等，均为流体动能周期的2倍。

在流体能量和管道能量的基础上，相加得到流体总能量与管道总能量，见图 8，变化趋势和周期几乎完全重合。当发生正水击时，二者都增加，发生负水击时，二者都减小。机理分析如下：对比能量方程中各项做功功率的量级，发现主导流体总能量变化的是外界输入功率W_f，主导管道总能量变化的是流体对管道的Poisson耦合效应功率P_t。阀门关闭后，Q(L, t)=0，对表 2中W_f的表达式积分，化简得到W_f=ρ_fgH(0, t)Q(0, t)。当系统中发生正水击时，上游水箱中的流体流入管内，Q(0, t)>0，因此W_f为正，流体的总能量增加；同时根据P_t的表达式，正水击时压强升高导致$\partial H / \partial t>0$，因此P_t为正，管道总能量增加。

4 能量因子分析

结合初始条件和边界条件，流固耦合4方程模型计算得到不同时刻所有节点的流体压强、流速、管道应力和速度。在此基础上结合流体能量和管道能量的表达式，采用复合的Simpson积分方法，对各节点物理量进行数值积分，得到整根管道的能量数值。进一步分析各能量的最大波动幅值和能量因子，研究各能量间的传递和转化机制，定量化描述流固耦合水击响应特性。

当阀门固支时，系统中存在Poisson耦合和摩擦耦合；当阀门简支时，阀门附近管道可以轴向自由位移，系统还存在结合部耦合。2种工况中流体能量和管道能量的变化趋势相似，但是波动幅值发生了较大的改变。因此定义流体内能、流体动能、管道应变能、管道动能、流体总能量和管道总能量的最大波动幅值分别为ΔI_f、ΔK_f、ΔU_t、ΔK_t、ΔE_f和ΔE_t，具体数值如表 3所示。

基于流固耦合4方程模型和能量方程，考虑结合部耦合后，管道下游的约束条件减弱，流固耦合的强度显著增大，能量传递和转换的剧烈程度陡增。在这个过程中，流体能量小幅增大，其中ΔI_f增大了约3%，ΔK_f增大了11%；管道能量大幅增大，其中ΔU_t增大了35%，ΔK_t增大到了阀门固定工况的68.62倍。分析结果表明，阀门处的固支和简支边界条件对流体能量的影响相对较小，而对管道能量的影响较大。同时结合表 3可得，无论是否考虑结合部耦合，能量波动幅值从大到小排序依次为ΔI_f、ΔK_f、ΔU_t和ΔK_t，再次验证了在流体能量中，内能占主导地位；在管道能量中，应变能占主导地位，因此存在ΔE_f≈ΔI_f，ΔE_t≈ΔU_t。

进一步定义无量纲的能量因子来描述流体和管道整体的水击瞬变响应特性，包括流固耦合因子ϕ_E、管道振动因子ϕ_t和液压脉动因子ϕ_f：

$ \phi_{\mathrm{E}}=\frac{\Delta E_{\mathrm{t}}}{\Delta E_{\mathrm{f}}}, $

(33)

$ \phi_{\mathrm{t}}=\frac{\Delta K_t}{\Delta U_{\mathrm{t}}}, $

(34)

$ \phi_{\mathrm{f}}=\frac{\Delta I_{\mathrm{f}}}{\Delta K_{\mathrm{f}}}. $

(35)

考虑结合部耦合后，ϕ_E从4.35%增大到5.78%，表明系统整体的流固耦合剧烈程度有较大提升；管道能量大幅增加，ϕ_t从0.91%显著增大到46.38%，表明管道运动的剧烈程度陡然增大；ϕ_f从8.04减小到7.46，表明流体压强脉动的剧烈程度相对有所减小，但变化幅度不大。

上述分析表明，结合部耦合是一种强烈的耦合方式，对管道的应力、速度、应变能和动能变化具有显著的影响，在弹性管路振动和压强波瞬态响应中起到了重要的作用，直接关系到管路系统动力学性能预测的精度和管道工作过程中的安全性。在此基础上将能量因子作为安全评价指标，则可以采取相关防护措施来提高管道系统的安全性。当ϕ_E较大时，需要减弱流固耦合的强度，在不改变材料性质的条件下，可以从结合部耦合入手，减少管道截面或边界条件突变的管道附件；当ϕ_t较大时，可以加强管道的固定约束，如增加固定支座的数量和约束的刚度；当ϕ_f较大时，可以在系统中增加消能器，减小水击能量造成的直接危害。

5 结论

本文建立了基于流固耦合水击理论的能量分析方法，引入了多样的能量种类，推导得到流体能量和管道能量的表达式，以及能量方程。在此基础上改变阀门处的边界条件，定义无量纲的能量因子，描述了结合部耦合对水击过程瞬变流动响应特性的影响。

能量分析方法从系统整体层面解释了水击过程中能量传递和转换的过程，提供了自然而直接的角度来理解系统的动态响应过程，这是传统的波传输和反射理论难以表现的。描述了流体能量和管道能量的变化关系，指出其中的主导能量分别为流体内能和管道应变能。同时揭示出流体总能量与管道总能量的同向变化关系，厘清了二者能量的来源。考虑结合部耦合后，加剧了能量的传递和转换，流体能量小幅增大，管道能量大幅增大，导致流固耦合因子有所增大，管道振动因子显著增大，液压脉动因子略微减小。进一步将三个能量因子作为安全评价指标，提出了管道系统的防护措施。能量方程是由流固耦合水击模型推导获得的，因此与流固耦合模型是线性相关的，但是能量分析方法提供了一种理解和对比水击过程动态响应的新方法，为定量化描述流固耦合水击响应特性提供了方向和依据。