管道系统仿真

管道输送系统仿真
一实验题目及内容：
下图所示的管道，D＝0.2m，L＝3000m，Le＝3050m，H0＝60m，H1＝20m，自由流出，流动处于紊流区，m=0.125，f＝2.7 (S/m3)1.875，a=1000m/S。

试进行终端阀门瞬时关闭的水击计算。

二项目研究背景及意义：
随着我国石油工业的发展，油气长输管线的建设成为国家“十五”、“十一五”的重点建设项目。

未来几年内，我国将建成14跳油气输送通道，形成“两纵、亮横、四枢纽、五气库”，总长超过1万公里的油气逛到输送格局。

近年来，老的油气输送管道逐步进行自动化控制改造，新建的长输油气管道均采用先进的战地黄花控制系统，这对长输油气管道的管理和操作人员素质提出了更高的要求。

同时，由于生产运行必须保证安全、稳定、长期、优化地运行，因而在生产实际过程中不放兵依靠实际生产装置来提高和培训运行人员的操作技能，采用仿真系统则可以很有效地解决在岗人员的操作技能培训问题。

但是，到目前为止，我国还没有长输原油管道密闭传输自动化控制仿真系统。

通过使用长输油气管道自动控制仿真系统，能让操作人员更深入地了解生产装置的工艺流程；了解并测试各类自动检测仪表对油田生产过程的适应性；研究适合于长输管道系统特点的优化控制方案；提高生产岗位操作人员的操作技能和技术素质，为生产安全、稳定、长期、优化地运行提供保证。

三编程软件的选择及简介：
本次试验所选用的的编程软件为Visual Basic，Visual Basic是一种由微软公司开发的包含协助开发环境的事件驱动编程语言。

从任何标准来说，VB都是世界上使用人数最多的语言。

它源自于BASIC编
程语言。

VB拥有图形用户界面（GUI）和快速应用程序开发（RAD）系统，可以轻易的使用DAO、RDO、ADO连接数据库，或者轻松的创建ActiveX控件。

程序员可以轻松的使用VB提供的组件快速建立一个应用
四解题思路：
首先画x-t网格。

将管道分成三段（N＝3），先进行常数计算，接着计算稳态工况。

结束稳态、开始瞬变的时刻定为t＝0。

终端阀门于t=1s关闭，在此时刻只有右边界点受到扰动。

扰动逐渐向管道另一端传递，进而反弹，按照周期变化。

五实现过程；计算结果：
解：首先画x-t网格。

将管道分成三段（N＝3），每段长度为
Δx=L/N=1000m，节点编号为0、1、2、3，时步Δt=Δx/a=1s。

先进行常数计算：
由于L≠Le，计算差分管段的摩阻时应取fLe/L=2.745s1.875/m5.625=5.8949×10-7 h1.875/m5.625。

各节点的压头计算过程从略，只将结果标在水力坡降线的上方。

结束
稳态、开始瞬变的时刻定为t＝0。

终端阀门于t=1s关闭。

在此时刻只有右边界点受到扰动，按式（4—37）和式〔4一38〕，
QP3＝0
HP3＝33.33＋0.9013×210.08＝33.33+189.35＝222.7m
六计算结果:
七程序详细代码: Private Sub Form_click() Dim d As Single
Dim l As Single
Dim le As Single
Dim f As Single
Dim a As Single
Dim cw As Single
Dim g As Single
Dim dx As Integer
Dim m As Single
Dim t As Integer
Dim i As Integer
d = 0.2
l = 300
le = 3050
f = 2.7
a = 1000
g = 9.8
dx = 1000
m = 0.125
cw = a / (g * 3.1415 * d * d / 4) Dim q(3, 7) As Single
Dim h(3, 7) As Single
h(0, 0) = 60
h(1, 0) = 46.7
h(2, 0) = 33.3
h(3, 0) = 20
q(0, 1) = 102.1 / 3600
q(1, 1) = 210.1 / 3600
q(2, 1) = 210.1 / 3600
q(3, 1) = 0
h(0, 1) = 60
h(1, 1) = 46.7
h(2, 1) = 33.3
h(3, 1) = 222.7
For i = 0 To 3
For t = 1 To 6
q(i, 0) = 210.1 / 3600
If i = 0 Then
h(0, t + 1) = 60
q(0, t + 1) = (60 - h(1, t) + cw * q(1, t)) / (cw + dx * f * (Abs(q(1, t))) ^ (1 - m))
ElseIf i = 3 Then
q(3, t) = 0
h(3, t + 1) = h(2, t) + cw * q(2, t)
Else
q(i, t + 1) = ((-f * dx * (q(i - 1, t) * (Abs(q(i - 1, t))) ^ (1 - m) + q(i + 1, t) * (Abs(q(i + 1, t))) ^ (1 - m))) / cw + q(i - 1, t) + q(i + 1, t)) / 2
h(i, t + 1) = -f * dx * q(i - 1, t) * (Abs(q(i - 1, t))) ^ (1 - m) - cw * (q(i, t) - q(i + 1, t)) + h(i - 1, t)
End If
Next t
Next i
For i = 0 To 3
Print Fix(10 * q(i, 0) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 0)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 1) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 1)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 2) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 2)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 3) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 3)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 4) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 4)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 5) * 3600) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 5)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 6) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 6)) / 10,
Next i
Print
For i = 0 To 3
Print Fix(10 * q(i, 7) * 3600) / 10, Next i
Print
For i = 0 To 3
Print Fix(10 * h(i, 7)) / 10,
Next i
Print
Print End Sub
11。