首页> 中国专利> 一种基于流固耦合的4方程模型的改进方法

一种基于流固耦合的4方程模型的改进方法

页面导航

摘要
著录项
法律信息
说明书
相似文献

摘要

本发明涉及一种基于流固耦合的4方程模型的改进方法，属于流固耦合理论及计算领域。本发明首先采用

著录项

公开/公告号CN104573336A

专利类型发明专利
公开/公告日2015-04-29

原文格式PDF
申请/专利权人昆明理工大学;
展开▼

申请/专利号CN201410814639.3
发明设计人陈婷;苏志敏;朱建兵;皇甫飞华;冯赛;
展开▼

申请日2014-12-24
分类号G06F19/00;
代理机构
代理人
地址 650093 云南省昆明市五华区学府路253号
入库时间 2023-12-18 08:25:28

法律信息

法律状态公告日

法律状态信息

法律状态
2017-08-25

授权

授权
2015-05-27

实质审查的生效 IPC(主分类):G06F19/00 申请日:20141224

实质审查的生效
2015-04-29

公开

公开

说明书

技术领域

本发明涉及一种基于流固耦合的4方程模型的改进方法，属于流固耦合理论及计算领域。

背景技术

传统水击理论所用的连续性方程适用面较广，可用于任何恒定流或非恒定流的水力计算，但在发生水击时，管道内并存在液体压力波速、管道应力波速、流体流速，而经典的连续性方程并未能在微分方程中反映这个情况。另外，传统水击计算理论主要重点是研究流体的动力学行为对结构的影响分析，忽略了由于流体对结构运动状态改变而产生的流体运动变化，并且进行了大量的简化处理，这样导致一部分重要的系统信息丢失，不能更好的反应管道系统的实际运动状态。

发明内容

本发明提供了一种基于流固耦合的4方程模型的改进方法，以用反映管道内同时存在的液体压力波速、管道应力波速、流体流速以及管道自身的振动，所得模型在理论上的更严谨更为精确。

本发明的技术方案是：一种基于流固耦合的4方程模型的改进方法，首先采用对传统的4方程模型的连续性微分方程进行推导，得到改进的流体连续方程；接着对改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程；再对用于计算耦合水击的基本连续性方程进行简化；然后忽略简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程；再根据得到的连续性方程建立改进的4方程模型；接着对改进的4方程模型变换成矩阵表达形式；最后利用特征线法对矩阵表达形式进行求解，得其特征方程。

所述方法的具体步骤如下：

Step1、对传统的4方程模型的连续性微分方程进行推导时，定义得到改进的流体连续方程：

$2 (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} + (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f})}{\partial t} = 0$

式中，λ为改进后4方程中利用特征线法所求得的解，v_f为上游端断面的流速，ρ_f为上游端断面的密度，A_f为上游端断面面积，z为z轴坐标轴，t为时间，ds为微元段长度；

Step2、对步骤Step1中的改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程：

Step2.1、对步骤Step1中的改进的流体连续方程进行变换，再对变换后的方程进一步分析：

变换后的方程：

$[(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t}] + [(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f})}{\partial t} + (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z}] = 0$

进一步分析后的方程：

$[(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t}] + (λ \pm v_{f}) [ρ_{f} \frac{D A_{f}}{Dt} + A_{f} \frac{D ρ_{f}}{Dt} + ρ_{f} A_{f} v_{f}^{'}] = 0$

Step2.2、基于非恒定流的连续方程微分形式则得到用于计算耦合水击的基本连续性方程如下：

$(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} = 0$

Step3、对用于计算耦合水击的基本连续性方程进行简化，具体为：将用于计算耦合水击的基本连续性方程展开，接着将液体弹性模量K_f和液体截面面积A_f随时间的变化率、管道单元体轴向总应变、管道轴向内力对时间的导数方程代入用于计算耦合水击的基本连续性方程；接着忽略代入后的方程中各项的高次微量；再用压力水头代替平均压强p,p＝ρ_fgH+ρ_fgz>0；然后忽略压力梯度；最后采用由young推导的 $c_{f}^{2} = \frac{1}{ρ_{f} (\frac{1}{K_{f}} + \frac{2 R}{Eδ})}$ 简化，得到：

$v_{f}^{'} + \frac{g}{c_{f}^{2}} \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} + \frac{1}{(λ \pm v_{f})} {\overset{\cdot}{v}}_{f} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{{\overset{\cdot}{ρ}}_{f}}{ρ_{f}} + \frac{{\overset{\cdot}{A}}_{f}}{A_{f}}) = 0$

式中，g为重力加速度，H为液体压力水头，α为管轴线与水平面的夹角，p₀为液体压力，c_f为液体压力波速，R为管道内半径，E为管材弹性模量，δ为管壁厚度，μ为管材泊松比，A_p为管材横截面积，u_z为z轴方向的位移；

Step4、忽略步骤Step3简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程为：

$v_{f}^{'} + \frac{g}{c_{f}^{2}} \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} = 0$

Step5、根据步骤Step4得到的连续性方程建立改进的4方程模型如下：

流体动量方程： $\frac{\partial v_{f}}{\partial t} + g \frac{\partial H}{\partial z} = - \frac{f}{4 R} | v_{r} | v_{r}$

连续方程为： $\frac{\partial v_{f}}{\partial z} + \frac{g}{c_{f}^{2}} \frac{\partial H}{\partial t} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} \frac{{\overset{\cdot}{u}}_{z}}{\partial z} = 0$

管道动力方程： $\frac{\partial {\overset{\cdot}{u}}_{z}}{\partial t} - \frac{1}{ρ_{p}} \frac{\partial σ_{z}}{\partial z} = \frac{ρ_{f} A_{f}}{ρ_{p} A_{p}} \frac{f ρ_{f}}{4 R} v_{r} | v_{r} | - \frac{c_{d}}{ρ_{p} A_{p}} {\overset{\cdot}{u}}_{z}$

管道物理方程： $\frac{\partial {\overset{\cdot}{u}}_{z}}{\partial z} - \frac{1}{E} \frac{\partial σ_{z}}{\partial t} + \frac{g ρ_{f} μR}{Eδ} \frac{\partial H}{\partial t} = 0$

式中，f为液体与管壁之间的磨擦阻尼系数，v_r为液体与管道轴向相对流速，σ_z为管道在z轴方向的应力，c_d为管道在z轴方向的结构阻尼系数，ρ_p为管材密度；

Step6、对步骤Step5中改进的4方程模型变换成矩阵表达形式如下：

$A \frac{\partial}{\partial z} (\begin{matrix} v_{f} \\ H \\ {\overset{\cdot}{u}}_{z} \\ σ_{z} \end{matrix}) + B \frac{\partial}{\partial t} (\begin{matrix} v_{f} \\ H \\ {\overset{\cdot}{u}}_{z} \\ σ_{z} \end{matrix}) = F$

式中，A、B是系数矩阵，F是力向量，表达式如下：

$A = (\begin{matrix} 0 & g & 0 & 0 \\ 1 & 0 & - \frac{μR}{δ} \frac{A_{p}}{A_{f}} & 0 \\ 0 & 0 & 0 & - \frac{1}{ρ_{p}} \\ 0 & 0 & 1 & 0 \end{matrix}),$

$B = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{g}{c_{f}^{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & \frac{g ρ_{f} μR}{Eδ} & 0 & - \frac{1}{E} \end{matrix}),$

$F = (\begin{matrix} - \frac{f}{4 R} | v_{r} | v_{r} \\ 0 \\ - \frac{c_{d}}{ρ_{p} A_{p}} {\overset{\cdot}{u}}_{z} + \frac{f}{4 R} \frac{A_{f} ρ_{f}}{A_{p} ρ_{p}} | v_{r} | v_{r} \\ 0 \end{matrix});$

Step7、利用特征线法对步骤Step6求解，可得其特征方程为：

$λ^{4} - λ^{2} c_{f}^{2} [1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}}] + c_{f}^{2} c_{p}^{2} = 0;$

令 $p^{2} = 1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}},$ 求解出相应的特征值如下：

${\tilde{c}}_{f} = λ_{1} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{2} = - λ_{1};$

${\tilde{c}}_{p} = λ_{3} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{4} = - λ_{3};$

式中为液体耦合压力波速，为管道耦合应力波速，λ₁、λ₂、λ₃、λ₄为改进后4方程中利用特征线法所求得的解。

本发明的工作原理是：

①、取长度为ds微元段的流体为研究对象，如图1所示的1-1断面为上游端，2-2断面为下游端，假设v_f为1-1断面的流速、p为压强、ρ_f为液体密度、A_f为断面面积，则2-2断面相应的流速为压强为密度为断面面积为当下游端阀门关闭产生水击波向上游传播，以波速c_f从2-2断面向1-1断面传播，此时需要考虑流体自身的流动，经过dt时段后，流体微元流动的距离设为dΔ＝v_fdt，那么原流段到达新的位置即1′-1′断面与2′-2′断面之间。

其中，s为管道长度。

②、由1断面流进与从2断面流出的质量差可求得时段末流段1-2内液体质量为

$(\begin{matrix} \frac{1}{2} [ρ_{f} + (ρ_{f} + \frac{\partial ρ_{f}}{\partial s} ds) + \frac{\partial}{\partial t} (ρ_{f} + \frac{\partial ρ_{f}}{\partial s} ds) dt \\ \cdot \frac{1}{2} [A_{f} + (A_{f} + \frac{\partial A_{f}}{\partial s} ds) + \frac{\partial}{\partial t} (A_{f} + \frac{\partial A_{f}}{\partial s} ds) dt] ds \\ = \frac{1}{4} [2 ρ_{f} + \frac{\partial ρ_{f}}{\partial s} ds + \frac{\partial ρ_{f}}{\partial t} dt] [2 A + \frac{\partial A}{\partial s} ds + \frac{\partial A}{\partial t} dt] ds \\ = [ρ_{f} A_{f} + \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial s} ds + \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial t} dt] ds \end{matrix}) - - - (1)$

求得时段dt内流段1-2内液体质量的增量为

$(\begin{matrix} [ρ_{f} A_{f} + \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial s} ds + \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial t} dt] ds - (ρ_{f} A_{f} ds + \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial s} {ds}^{2}) \\ = \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial t} dt \cdot ds \end{matrix}) - - - (2)$

③、根据质量守恒原理，(1)式等于(2)式，得

$- [\frac{\partial (ρ_{f} A_{f} v_{f})}{\partial s} ds + \frac{1}{2} \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} dt] dt = \frac{1}{2} \frac{\partial (ρ_{f} A_{f})}{\partial t} dt \cdot ds - - - (3)$

定义整理得

$\frac{\partial (ρ_{f} A_{f} v_{f})}{\partial s} (c_{f} \pm v_{f}) + \frac{1}{2} \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} + \frac{1}{2} \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} (c_{f} \pm v_{f}) = 0$

由于c_f＞＞v_f，并将研究对象移至既定的坐标系中，使得管轴线与坐标z轴重合。即所取的管道长度ds与前面所建立的水击模型的在同一坐标系中，那么连续性方程可简化改为：

④、由所述步骤③可知，对传统的4-方程模型的连续性微分方程推导时，定义其中为Korteweg波速，ζ为Korteweg波速修正系数，K为弹性模量，R为管道内半径，E为管材弹性模量，δ为管壁厚度，z为z轴坐标轴。其仅反映了液体压力波和流体流速，对于流体与管道耦合振动的情况还应该有管道应力波的存在，所以，将③步骤推导过程中的将连续性方程推导中的更正为能够反映耦合水击特的管道与流体在纵横两向均耦合的耦合波速，即其中λ为特征波速，即为改进后4方程利用特征线法所求得的解。那么流体连续方程变换为：

$2 (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} + (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f})}{\partial t} = 0 - - - (5)$

⑤、对所述步骤④的(5)式进一步分析如下：

由于非恒定流的连续方程微分形式则得到用于计算耦合水击的基本连续性方程如下：

$(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} = 0 - - - (6)$

上式由于加入了λ使得本发明所采用的连续性方程能反映横纵向耦合水击特性。

⑥、为了从量化上来证明模型的准确性，将所述步骤⑤中(6)展开，得到：

$(λ \pm v_{f}) (ρ_{f} A_{f} \frac{\partial v_{f}}{\partial z} + ρ_{f} v_{f} \frac{\partial A_{f}}{\partial z} + A_{f} v_{f} \frac{\partial ρ_{f}}{\partial z}) + (ρ_{f} A_{f} \frac{\partial v_{f}}{\partial t} + ρ_{f} v_{f} \frac{\partial A_{f}}{\partial t} + A_{f} v_{f} \frac{\partial ρ_{f}}{\partial t}) = 0$

整理得到：

$(λ \pm v_{f}) (ρ_{f} A_{f} \frac{\partial v_{f}}{\partial z} + ρ_{f} v_{f} \frac{\partial A_{f}}{\partial z}) + (ρ_{f} A_{f} \frac{\partial v_{f}}{\partial t} + ρ_{f} v_{f} \frac{\partial A_{f}}{\partial t}) + [(λ \pm v_{f}) A_{f} v_{f} \frac{\partial ρ_{f}}{\partial z} + A_{f} v_{f} \frac{\partial ρ_{f}}{\partial t}] = 0$ 上式两边同除于λρ_fA_f得到：

$\frac{\partial v_{f}}{\partial z} + \frac{v_{f}}{(λ \pm v_{f}) A_{f}} \frac{\partial A_{f}}{\partial t} + \frac{v_{f}}{A_{f}} \frac{\partial A_{f}}{\partial z} + \frac{1}{(λ \pm v_{f})} \frac{\partial v_{f}}{\partial t} + \frac{1}{ρ_{f}} [\frac{v_{f}}{(λ \pm v_{f})} \frac{\partial ρ_{f}}{\partial t} + v_{f} \frac{\partial ρ_{f}}{\partial z}] = 0 - - - (7)$

已知液体弹性模量K_f和液体截面面积A_f随时间的变化率如下

$K_{f} = ρ_{f} (\frac{dp}{dt}) / (\frac{d ρ_{f}}{dt}),$ 也可以表达为 $\frac{\partial ρ_{f}}{\partial t} + v_{f} \frac{\partial ρ_{f}}{\partial s} = \frac{ρ_{f}}{K_{f}} \frac{dp}{dt}$

$\frac{1}{A_{f}} \frac{d A_{f}}{dt} = \frac{2 R}{Eδ} (\frac{dp}{dt} - \frac{μ}{2 A_{f}} \frac{dF}{dt}),$ 也可以表达为 $\frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t} + \frac{v_{f}}{A_{f}} \frac{\partial A_{f}}{\partial z} = \frac{2 R}{Eδ} (\frac{dp}{dt} - \frac{μ}{2 A_{f}} \frac{dF}{dt})$

式中，μ为管材泊松比；为管道轴向内力对时间的导数；

将以上两式代入(7)式，可得

$(\begin{matrix} \frac{\partial v_{f}}{\partial z} + \frac{v_{f}}{(λ \pm v_{f}) A_{f}} \frac{\partial A_{f}}{\partial t} + [\frac{2 R}{Eδ} (\frac{dp}{dt} - \frac{μ}{2 A_{f}} \frac{dF}{dt}) - \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t}] + \frac{1}{(λ \pm v_{f})} + \\ \frac{1}{ρ_{f}} [\frac{v_{f}}{(λ \pm v_{f})} \frac{\partial ρ_{f}}{\partial t} + (\frac{ρ_{f}}{K_{f}} \frac{dp}{dt} - \frac{\partial ρ_{f}}{\partial t})] = 0 \end{matrix})$

整理可得：

$(\begin{matrix} \frac{\partial v_{f}}{\partial z} + (\frac{v_{f}}{λ \pm v_{f}} - 1) \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t} + \frac{2 R}{Eδ} (\frac{dp}{dt} - \frac{μ}{2 A_{f}} \frac{dF}{dt}) + \frac{1}{λ \pm v_{f}} \frac{\partial v_{f}}{\partial t} + \\ \frac{1}{ρ_{f}} (\frac{v_{f}}{λ \pm v_{f}} - 1) \frac{\partial ρ_{f}}{\partial t} + \frac{1}{K_{f}} \frac{dp}{dt} = 0 \end{matrix})$

则连续方程变换为如下的形式：

$(\begin{matrix} \frac{\partial v_{f}}{\partial z} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) \frac{dp}{dt} - \frac{R}{Eδ} \frac{μ}{A_{f}} \frac{dF}{dt} + \\ \frac{1}{λ \pm v_{f}} \frac{\partial v_{f}}{\partial t} + (\frac{v_{f}}{λ \pm v_{f}} - 1) (\frac{1}{ρ_{f}} \frac{\partial ρ_{f}}{\partial t} + \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t}) = 0 \end{matrix}) - - - (8)$

管道单元体轴向总应变为：

$ϵ_{z} \approx u_{z}^{'} + ({u_{x}^{'}}^{2} + {u_{y}^{'}}^{2}) / 2 - 3 u_{z}^{'} ({u_{x}^{'}}^{2} + {u^{' 2}}_{y}) / 2$

其对时间的导数方程为：

$\frac{d ϵ_{z}}{dt} = \frac{\partial ϵ_{z}}{\partial t} + v_{f} \frac{\partial ϵ_{z}}{\partial z} \approx {\overset{\cdot}{u}}_{z}^{'} + {\overset{\cdot}{u}}^{'} u_{x}^{'} + {\overset{\cdot}{u}}^{'} u_{y}^{'} + v_{y} (u_{z}^{''} + u_{x}^{''} u_{x}^{'} + u_{y}^{''} u_{y}^{'})$

式中，u_x、u_y、u_z为管道在x、y、z轴方向的位移

管道轴向内力对时间的导数方程为：

$\frac{dF}{dt} = \frac{d (F_{0} + E A_{p} ϵ_{z})}{dt} = E A_{p} \frac{d ϵ_{z}}{dt}$

则 $\frac{dF}{dt} = E A_{p} [{\overset{\cdot}{u}}_{z}^{'} + {\overset{\cdot}{u}}^{'} u_{x}^{'} + {\overset{\cdot}{u}}^{'} u_{y}^{'} + v_{f} {(u_{z}^{''} + u_{x}^{''} u_{x}^{'} + u_{y}^{''} u^{'})}_{y}]$

式中，F₀为作用在管道轴向的初始外力，A_p为管道截面面积；

将上两式代入(8)式可得：

$\frac{\partial v_{f}}{\partial z} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) \frac{dp}{dt} - \frac{R}{Eδ} \frac{μ}{A_{f}} \frac{dF}{dt} + \frac{1}{(λ \pm v_{f})} \frac{\partial v_{f}}{\partial t} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{1}{ρ_{f}} \frac{\partial ρ_{f}}{\partial t} + \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t}) = 0$

将所有参数都展开可得：

$(\begin{matrix} \frac{\partial v_{f}}{\partial z} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) \frac{dp}{dt} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} [{\overset{\cdot}{u}}_{z}^{'} + {\overset{\cdot}{u}}^{'} + {\overset{\cdot}{u}}^{'} u_{y}^{'} + v_{f} {(u_{z}^{''} + u_{x}^{''} u_{x}^{'} + u_{y}^{''} u^{'})}_{y}] + \\ \frac{1}{(λ \pm v_{f})} \frac{\partial v_{f}}{\partial t} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{1}{ρ_{f}} \frac{\partial ρ_{f}}{\partial t} + \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t}) = 0 \end{matrix})$

将上式忽略高次微量简化可得：

$\frac{\partial v_{f}}{\partial z} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) \frac{dp}{dt} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} + \frac{1}{(λ \pm v_{f})} \frac{\partial v_{f}}{\partial t} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{1}{ρ_{f}} \frac{\partial ρ_{f}}{\partial t} + \frac{1}{A_{f}} \frac{\partial A_{f}}{\partial t}) = 0$

将上式整理并用指标表示：

$v_{f}^{'} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) (\overset{\cdot}{p} + v_{f} p^{'}) - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} + \frac{1}{(λ \pm v_{f})} {\overset{\cdot}{v}}_{f} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{{\overset{\cdot}{ρ}}_{f}}{ρ_{f}} + \frac{{\overset{\cdot}{A}}_{f}}{A_{f}}) = 0$

忽略压力梯度，则有

$v_{f}^{'} + (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) \overset{\cdot}{p} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} + \frac{1}{(λ \pm v_{f})} {\overset{\cdot}{v}}_{f} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{{\overset{\cdot}{ρ}}_{f}}{ρ_{f}} + \frac{{\overset{\cdot}{A}}_{f}}{A_{f}}) = 0$

接着，用压力水头代替平均压强将p,p＝ρ_fgH+ρ_fgz>0，并忽略压力梯度则上式变为：

$v_{f}^{'} = (\frac{2 R}{Eδ} + \frac{1}{K_{f}}) ρ_{f} g \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} + \frac{1}{(λ \pm v_{f})} {\overset{\cdot}{v}}_{f} + [\frac{v_{f}}{(λ \pm v_{f})} - 1] (\frac{{\overset{\cdot}{ρ}}_{f}}{ρ_{f}} + \frac{{\overset{\cdot}{A}}_{f}}{A_{f}}) = 0$

采用由young推导的进行变换，可得：

式中，g为重力加速度，H为液体压力水头，α为管轴线与水平面的夹角，p₀为液体压力，c_f为液体压力波速，u_z为z轴方向的位移；

⑦、对于4方程而言，由于未知量过多导致方程无法求解，故将所述步骤⑥中(9)式忽略流体密度及液体截面面积随时间的变化率与液体速度变化小量，则连续方程为：

$v_{f}^{'} + \frac{g}{c_{f}^{2}} \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} = 0$

⑧、与经典4方程中简化后的流体动量方程、管道运动方程及物理方程构成了改进的轴向4-方程模型如下：

流体动量方程： $\frac{\partial v_{f}}{\partial t} + g \frac{\partial H}{\partial z} = - \frac{f}{4 R} | v_{r} | v_{r} - - - (10)$

连续方程为： $\frac{\partial v_{f}}{\partial z} + \frac{g}{c_{f}^{2}} \frac{\partial H}{\partial t} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} \frac{{\overset{\cdot}{u}}_{z}}{\partial z} = 0 - - - (11)$

管道物理方程： $\frac{\partial {\overset{\cdot}{u}}_{z}}{\partial z} - \frac{1}{E} \frac{\partial σ_{z}}{\partial t} + \frac{g ρ_{f} μR}{Eδ} \frac{\partial H}{\partial t} = 0 - - - (13)$

式中，f为液体与管壁之间的磨擦阻尼系数，v_r为液体与管道轴向相对流速，σ_z为管道在z轴方向的应力，c_d为管道在z轴方向的结构阻尼系数，ρ_p为液体密度；

⑨、所述步骤⑧中改进轴向4方程(10)、(11)、(12)、(13)方程的另一种矩阵表达形式如下：

式中，A、B是系数矩阵，F是力向量，表达式如下：

$A = (\begin{matrix} 0 & g & 0 & 0 \\ 1 & 0 & - \frac{μR}{δ} \frac{A_{p}}{A_{f}} & 0 \\ 0 & 0 & 0 & - \frac{1}{ρ_{p}} \\ 0 & 0 & 1 & 0 \end{matrix}), B = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{g}{c_{f}^{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & \frac{g ρ_{f} μR}{Eδ} & 0 & - \frac{1}{E} \end{matrix}),$

⑩、利用特征线法对所述步骤⑨求解可得其特征方程为

$λ^{4} - λ^{2} c_{f}^{2} [1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}}] + c_{f}^{2} c_{p}^{2} = 0 - - - (15)$

令 $p^{2} = 1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}},$ 求解出相应的特征值如下：

${\tilde{c}}_{f} = λ_{1} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{2} = - λ_{1}; - - - (16)$

${\tilde{c}}_{p} = λ_{3} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{4} = - λ_{3} - - - (17)$

其中，为液体耦合压力波速，为管道耦合应力波速，λ₁、λ₂、λ₃、λ₄为改进后4方程中利用特征线法所求得的解。

其中，传统的4方程模型参照：杨超.非恒定流充液管系统耦合振动特性及振动抑制[D].武汉:华中科技大学，2007。

本发明的有益效果是：更能反映管道内同时存在的液体压力波速、管道应力波速、流体流速以及管道自身的振动，所得模型在理论上的更严谨更为精确。

附图说明

图1为本发明中管道简图；

图2为本发明中不同管壁厚的液体耦合压力波速一；

图3为本发明中不同管壁厚的液体耦合压力波速二；

图4为本发明中不同管壁厚的管道耦合应力波速；

图5为本发明中不同管道半径的液体耦合压力波速一；

图6为本发明中不同管道半径的液体耦合压力波速二；

图7为本发明中不同管道半径的管道耦合应力波速。

具体实施方式

实施例1：如图1-7所示，一种基于流固耦合的4方程模型的改进方法，首先采用对传统的4方程模型的连续性微分方程进行推导，得到改进的流体连续方程；接着对改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程；再对用于计算耦合水击的基本连续性方程进行简化；然后忽略简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程；再根据得到的连续性方程建立改进的4方程模型；接着对改进的4方程模型变换成矩阵表达形式；最后利用特征线法对矩阵表达形式进行求解，得其特征方程。

所述方法的具体步骤如下：

Step1、对传统的4方程模型的连续性微分方程进行推导时，定义得到改进的流体连续方程：

$2 (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} + (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f})}{\partial t} = 0$

Step2、对步骤Step1中的改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程：

Step2.1、对步骤Step1中的改进的流体连续方程进行变换，再对变换后的方程进一步分析：

变换后的方程：

进一步分析后的方程：

Step2.2、基于非恒定流的连续方程微分形式则得到用于计算耦合水击的基本连续性方程如下：

$(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} = 0$

Step4、忽略步骤Step3简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程为：

$v_{f}^{'} + \frac{g}{c_{f}^{2}} \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} = 0$

Step5、根据步骤Step4得到的连续性方程建立改进的4方程模型如下：

流体动量方程： $\frac{\partial v_{f}}{\partial t} + g \frac{\partial H}{\partial z} = - \frac{f}{4 R} | v_{r} | v_{r}$

连续方程为： $\frac{\partial v_{f}}{\partial z} + \frac{g}{c_{f}^{2}} \frac{\partial H}{\partial t} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} \frac{{\overset{\cdot}{u}}_{z}}{\partial z} = 0$

管道物理方程： $\frac{\partial {\overset{\cdot}{u}}_{z}}{\partial z} - \frac{1}{E} \frac{\partial σ_{z}}{\partial t} + \frac{g ρ_{f} μR}{Eδ} \frac{\partial H}{\partial t} = 0$

Step6、对步骤Step5中改进的4方程模型变换成矩阵表达形式如下：

式中，A、B是系数矩阵，F是力向量，表达式如下：

$A = (\begin{matrix} 0 & g & 0 & 0 \\ 1 & 0 & - \frac{μR}{δ} \frac{A_{p}}{A_{f}} & 0 \\ 0 & 0 & 0 & - \frac{1}{ρ_{p}} \\ 0 & 0 & 1 & 0 \end{matrix}),$

$B = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{g}{c_{f}^{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & \frac{g ρ_{f} μR}{Eδ} & 0 & - \frac{1}{E} \end{matrix}),$

Step7、利用特征线法对步骤Step6求解，可得其特征方程为：

$λ^{4} - λ^{2} c_{f}^{2} [1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}}] + c_{f}^{2} c_{p}^{2} = 0;$

令 $p^{2} = 1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}},$ 求解出相应的特征值如下：

${\tilde{c}}_{f} = λ_{1} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{2} = - λ_{1};$

${\tilde{c}}_{p} = λ_{3} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{4} = - λ_{3};$

式中，为液体耦合压力波速，为管道耦合应力波速，λ₁、λ₂、λ₃、λ₄为改进后4方程中利用特征线法所求得的解。

实施例2：如图1-7所示，一种基于流固耦合的4方程模型的改进方法，首先采用对传统的4方程模型的连续性微分方程进行推导，得到改进的流体连续方程；接着对改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程；再对用于计算耦合水击的基本连续性方程进行简化；然后忽略简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程；再根据得到的连续性方程建立改进的4方程模型；接着对改进的4方程模型变换成矩阵表达形式；最后利用特征线法对矩阵表达形式进行求解，得其特征方程。

实施例3：如图1-7所示，一种基于流固耦合的4方程模型的改进方法，首先采用对传统的4方程模型的连续性微分方程进行推导，得到改进的流体连续方程；接着对改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程；再对用于计算耦合水击的基本连续性方程进行简化；然后忽略简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程；再根据得到的连续性方程建立改进的4方程模型；接着对改进的4方程模型变换成矩阵表达形式；最后利用特征线法对矩阵表达形式进行求解，得其特征方程。

所述方法的具体步骤如下：

Step1、对传统的4方程模型的连续性微分方程进行推导时，定义得到改进的流体连续方程：

$2 (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} + (λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f})}{\partial t} = 0$

Step2、对步骤Step1中的改进的流体连续方程进行分析，得到用于计算耦合水击的基本连续性方程：

Step2.1、对步骤Step1中的改进的流体连续方程进行变换，再对变换后的方程进一步分析：

变换后的方程：

进一步分析后的方程：

Step2.2、基于非恒定流的连续方程微分形式则得到用于计算耦合水击的基本连续性方程如下：

$(λ \pm v_{f}) \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial z} + \frac{\partial (ρ_{f} A_{f} v_{f})}{\partial t} = 0$

Step4、忽略步骤Step3简化后的方程中的流体密度及流体截面面积随时间的变化率、液体速度变化小量，得到连续性方程为：

$v_{f}^{'} + \frac{g}{c_{f}^{2}} \overset{\cdot}{H} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} {\overset{\cdot}{u}}_{z}^{'} = 0$

Step5、根据步骤Step4得到的连续性方程建立改进的4方程模型如下：

流体动量方程： $\frac{\partial v_{f}}{\partial t} + g \frac{\partial H}{\partial z} = - \frac{f}{4 R} | v_{r} | v_{r}$

连续方程为： $\frac{\partial v_{f}}{\partial z} + \frac{g}{c_{f}^{2}} \frac{\partial H}{\partial t} - \frac{μR}{δ} \frac{A_{p}}{A_{f}} \frac{{\overset{\cdot}{u}}_{z}}{\partial z} = 0$

管道物理方程： $\frac{\partial {\overset{\cdot}{u}}_{z}}{\partial z} - \frac{1}{E} \frac{\partial σ_{z}}{\partial t} + \frac{g ρ_{f} μR}{Eδ} \frac{\partial H}{\partial t} = 0$

Step6、对步骤Step5中改进的4方程模型变换成矩阵表达形式如下：

式中，A、B是系数矩阵，F是力向量，表达式如下：

$A = (\begin{matrix} 0 & g & 0 & 0 \\ 1 & 0 & - \frac{μR}{δ} \frac{A_{p}}{A_{f}} & 0 \\ 0 & 0 & 0 & - \frac{1}{ρ_{p}} \\ 0 & 0 & 1 & 0 \end{matrix}),$

$B = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & \frac{g}{c_{f}^{2}} & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & \frac{g ρ_{f} μR}{Eδ} & 0 & - \frac{1}{E} \end{matrix}),$

Step7、利用特征线法对步骤Step6求解，可得其特征方程为：

$λ^{4} - λ^{2} c_{f}^{2} [1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}}] + c_{f}^{2} c_{p}^{2} = 0;$

令 $p^{2} = 1 + \frac{c_{p}^{2}}{c_{f}^{2}} + \frac{ρ_{f} μ^{2} R^{2}}{ρ_{p} δ^{2}} \frac{A_{p}}{A_{f}},$ 求解出相应的特征值如下：

${\tilde{c}}_{f} = λ_{1} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{2} = - λ_{1};$

${\tilde{c}}_{p} = λ_{3} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; λ_{4} = - λ_{3};$

式中，为液体耦合压力波速，为管道耦合应力波速，λ₁、λ₂、λ₃、λ₄为改进后4方程中利用特征线法所求得的解。

具体参数如下：

已知附图1钢管密度取7600kg/m³，管材弹性模量为201GPa，泊松比为0.27，管壁厚取8mm≤δ≤20mm，管道半径200mm≤R≤800mm，管道内流体为水。

由 ${\tilde{c}}_{f} = λ_{1} = c_{f} \sqrt{\frac{1}{2} [p^{2} - \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}; {\tilde{c}}_{p} = λ_{3} = c_{f} \sqrt{\frac{1}{2} [p^{2} + \sqrt{p^{4} - 4 \frac{c_{p}^{2}}{c_{f}^{2}}}]}$ 利用MATLAB编程可以分别计算得到不同连续方程下不同管壁厚、管道半径的液体的液体耦合压力波速及管道耦合应力波速，之后与杨玲霞、侯咏梅、Tijessing编写的文献《Fluid–structure>

表1不同管壁厚的液体耦合压力波速(m/s)

管壁厚(m)0.0080.010.0120.0140.0160.0180.02改进值836.028897.982947.867989.0691023.7661053.4351079.127文献值836.132898.112948.019989.2431023.9591053.6461079.354相对误差-0.104-0.13-0.152-0.174-0.193-0.211-0.227经典值857.763919.581969.0831009.7741043.8971072.9681098.058

由上表1可绘出附图2、3，可知，对于不同的管道壁厚，改进的轴向4方程与采用文献计算的得到的液体耦合压力波速非常接近，另外，因流体与管道间的耦合作用，使得改进值与文献值都比经典值低约20m/s。可以看到，附图中改进值与文献值几乎重叠，但从表1中可以看出改进值略比文献值降低，说明对4方程模型的连续方程改进后对水击波的影响比文献的略大，且随着管壁的厚度不断的增加，改进的4方程模型计算得到的液体耦合压力波速比文献模型所的压力波速的影响也逐渐变大。

表2不同管壁厚的管道耦合应力波速(m/s)

由上表2可绘出从附图4，可知，改进的轴向4方程与采用文献计算的得到的管道耦合应力波速同样很接近，而且都高于经典值，知经典管道耦合应力波速不受管道几何参数的影响。另外，可以看到，随着管壁厚度的增大，其管道耦合应力波速也在不断的减小。

表3不同管道半径的液体耦合压力波速(m/s)

管道半径0.20.30.40.50.60.70.8改进值1182.1491094.4331023.766965.2482915.751873.169836.029文献值1182.4501094.6701023.959965.410915.888873.288836.133相对误差-0.301-0.237-0.193-0.1618-0.137-0.119-0.104经典值1197.8591112.9691043.897986.270937.241894.866857.763

由上表3可绘出从附图5、6可知，对于不同的管道半径，改进的轴向4-方程与采用文献计算的得到的液体耦合压力波速其差值很小，同样，计算得到的液体耦合压力波速都小于经典值。另外，可以看到，随着管道半径的增大，其液体耦合压力波速不断的减小。

表4不同管道半径的管道耦合应力波速(m/s)

由上表4可绘出从附图7，可知，改进的轴向4方程与采用文献计算的得到的管道耦合应力波速其差值很小，计算所得管道耦合应力波速都高于经典值。随着管道半径的增大，其管道应力速波速增大。

综上可知，改进的轴向4方程与文献的数学模型所得结果是一致的，而且改进以后，流体与管道的耦合特性的反映比文献略有加强。另外，对于不同的管道壁厚，随着管道壁厚的增加，流体压力波不断增大，而对应的管道应力波则不断减小，对于不同管道半径而言，管道半径不断增大，液体耦合压力波速不断的减小，管道耦合应力速波速增大。从流固耦合的特性来讲，改进的连续性方程所得的结果更精确一些，所以改进的4方程模型用于耦合水击波的计算分析是可靠的、合理的。当然，对管道振动特性的影响还有待进一步的研究分析。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。

去获取专利，查看全文>

相似文献

专利
中文文献
外文文献

1. 一种基于流固耦合的4方程模型的改进方法 [P] . 中国专利： CN104573336B . 2017.08.25
2. 一种基于流固耦合的4方程模型的改进方法 [P] . 中国专利： CN104573336A . 2015-04-29
3. Airway ventilation analysis system using fluid structure coupled analysis of maxillofacial area [P] . 日本专利： JP6579472B2 . 2019-09-25

机译：基于颌面面积流固耦合分析的气道通气分析系统
4. METHOD, APPARATUS AND COMPUTER PROGRAM PRODUCT PROVIDING SIMULATOR FOR ENHANCED OIL RECOVERY BASED ON MICRON AND SUBMICRON SCALE FLUID-SOLID INTERACTIONS [P] . 美国专利： US2017298714A1 . 2017-10-19

机译：基于微米和亚微米级流固耦合的用于提高采收率的方法，装置和计算机程序产品模拟器
5. Airway ventilation analysis system using fluid structure coupled analysis of maxillofacial area [P] . 日本专利： JPWO2016006633A1 . 2017-05-25

机译：基于颌面面积流固耦合分析的气道通气分析系统