Лабораторные работы / Faz - 2006 / lab2
.docРХТУ им. Д.И. Менделеева
Кафедра информатики и компьютерного моделирования
Практическая работа №2
по курсу
«Компьютерное моделирование химических процессов»
Компьютерное моделирование простой гидравлической системы в динамическом режиме
Группа: О-41
Студент: Faz Birmutt :)
Вариант: 3
Преподаватель: Тамбовцев И.И.
Москва 2006
Схема гидравлической системы:
Известные величины:
Давления на входе и выходе системы: P1 - P6;
Коэффициенты сужающих устройств: k1 - k7;
Геометрические высоты емкостей: H1, H2.
Неизвестные величины:
Высоты столбов жидкости в емкостях: h1, h2;
Давления жидкости на дно сосуда: P7, P8;
Давления газа над жидкостью: P9, P10;
Расходы жидкости v1 - v7.
Основные допущения, накладываемые на модель:
1) Система работает в стационарном режиме;
2) Рассматривается только гидравлика процесса;
3) Сопротивление трения потоку пренебрежимо мало;
4) Жидкость несжимаема и имеет постоянную плотность;
5) Газ - идеальный pV=nRT;
6) Режим - изотермический T=const;
7) Трубопроводы находятся на одном уровне;
8) Цилиндрическая форма закрытой емкости с площадью поперечного сечения
S и геометрической высотой H;
9) Одинаковое давление газа Pн в емкостях, не заполненных жидкостью;
Система уравнений математического описания:
A) Уравнения расходов жидкости через клапаны:
1. v1 = k1 (p1 - p7)1/2
2. v2 = k2 (p2 - p7) 1/2
3. v3 = k3 (p3 - p7) 1/2
4. v4 = k4 (p8 - p4) 1/2
5. v5 = k5 (p8 - p5) 1/2
6. v6 = k6 (p8 - p6) 1/2
7. v7 = k7 (p7 - p8) 1/2
B) Определение давлений жидкости и газа в емкостях:
8. p9 = Pн H1 / (H1 - h1)
9. p7 = p9 + ρ g h1
10. p10 = Pн H2 / (H2 - h2)
11. p8 = p10 + ρ g h2
C) Дифференциальные уравнения заполнения емкостей:
12. dh1/dt = (v1 + v2 + v3 – v7) / s1
13. dh2/dt = (v7 – v4 – v5 – v6) / s2
D) Начальные условия:
14. h1(t=0) = h10
15. h2(t=0) = h20
Информационная матрица системы уравнений:
|
N п/п |
v1 |
v2 |
v3 |
v4 |
v5 |
v6 |
v7 |
p7 |
p8 |
p9 |
p10 |
h1 |
h01 |
h2 |
h02 |
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
4 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
5 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
6 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
7 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
8 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
9 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
10 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
11 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
12 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
13 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
14 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
15 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Блок-схема алгоритма расчета:
Текст программы:
Option Explicit
Option Base 1
Public culc As Integer
Const m As Byte = 2 ' Количество решаемых ОДУ
Const np% = 10, nk% = 7, nv% = 7, g! = 9.8
Public del! ' øàã
Private k!
Private y0!(m), x0! ' Начальное значение
Private n% ' Число шагов
Private kv% ' Кратность вывода
Private x1!, y1!(m), i%
Private ki!, ki1!, ki2!, vm!(nk), v!(nk), ak!(nk), p!(np)
Private hg!(m), h!(m), s!(m), pr!(m), ro!, pn!, x!, ipr%, ipr1%
Sub dydx(x As Single, y() As Single, pr!())
' расчет производных для текущих(pr) x,y
Dim j%
For j = 1 To m: h(j) = y(j): Next j
p(9) = pn * hg(1) / (hg(1) - h(1)): p(10) = pn * hg(2) / (hg(2) - h(2)) ''8 ''10
p(7) = p(9) + ro * g * h(1) * 0.000001: p(8) = p(10) + ro * g * h(2) * 0.000001 ''9 ''11
v(1) = ak(1) * Sgn(p(1) - p(7)) * Sqr(Abs(p(1) - p(7))) ''1
v(2) = ak(2) * Sgn(p(2) - p(7)) * Sqr(Abs(p(2) - p(7))) ''2
v(3) = ak(3) * Sgn(p(3) - p(7)) * Sqr(Abs(p(3) - p(7))) ''3
v(4) = ak(4) * Sgn(p(8) - p(4)) * Sqr(Abs(p(8) - p(4))) ''4
v(5) = ak(5) * Sgn(p(8) - p(5)) * Sqr(Abs(p(8) - p(5))) ''5
v(6) = ak(6) * Sgn(p(8) - p(6)) * Sqr(Abs(p(8) - p(6))) ''6
v(7) = ak(7) * Sgn(p(7) - p(8)) * Sqr(Abs(p(7) - p(8))) ''7
pr(1) = (v(1) + v(2) + v(3) - v(7)) / s(1): pr(2) = (v(7) - v(4) - v(5) - v(6)) / s(2) ''12 ''13
For j = 1 To nk: vm(j) = v(j) * ro: Next j
With Worksheets("TMP")
If ki > 0 And ki < 3 Then
If ki = 2 Then
.Cells(ipr, 1) = i: .Cells(ipr, 2) = "p(5-7)": .Cells(ipr, 6) = "vm"
.Cells(ipr, 3) = p(5): .Cells(ipr, 4) = p(6): .Cells(ipr, 5) = p(7)
.Cells(ipr, 7) = vm(1): .Cells(ipr, 8) = vm(2): .Cells(ipr, 9) = vm(3)
.Cells(ipr, 10) = vm(4): .Cells(ipr, 11) = vm(5): ipr = ipr + 1
End If
.Cells(ipr, 2) = "x": .Cells(ipr, 4) = "y(1-2)": .Cells(ipr, 7) = "pr(1-2)"
.Cells(ipr, 3) = x: .Cells(ipr, 5) = y(1): .Cells(ipr, 6) = y(2)
.Cells(ipr, 8) = pr(1): .Cells(ipr, 9) = pr(2): ipr = ipr + 1
End If
If ipr = 180 Then ipr = 1
End With
End Sub
Sub step(x As Single, y() As Single, del As Single, x1 As Single, y1() As Single)
' шаг интегрирования по X и Y расчет X1 и Y1
Dim y12!(m), j%
Call dydx(x, y, pr)
For j = 1 To m: y12(j) = y(j) + del * pr(j) / 2: Next j
Call dydx(x + del / 2, y12, pr): x1 = x + del
For j = 1 To m: y1(j) = y(j) + del * pr(j): Next j
End Sub
Public Sub gidra()
Dim j%, contr As String
Static x0s!, y0s!(2)
ipr = 1: ipr1 = 1
With Worksheets("TASK")
hg(1) = .Cells(4, 5): hg(2) = .Cells(5, 5): s(1) = .Cells(4, 9): s(2) = .Cells(5, 9)
ro = .Cells(6, 5): pn = .Cells(6, 9)
' Давление (1-6) МПа
For i = 1 To 6: p(i) = .Cells(8, i + 4): Next i
' Коэф. пpопускной способности (1-7)
For i = 1 To nk: ak(i) = .Cells(9, i + 4): Next i
' Начальные условия x0,y0(1),y0(2)
x0 = .Cells(10, 5): y0(1) = .Cells(10, 6): y0(2) = .Cells(10, 7)
If culc = 4 Then x0 = x0s: y0(1) = y0s(1): y0(2) = y0s(2)
' Число шагов, шаг, кратность вывода
n = .Cells(11, 5): del = .Cells(11, 6): kv = .Cells(11, 7)
' Относительная локальная погpешность ( % ) и Кpатность шагов на вывод
ki1 = .Cells(12, 5): ki2 = .Cells(13, 5)
End With
Worksheets("TMP").Activate
Cells.Select
Selection.Clear
Range("a1").Select
Worksheets("RESULTS").Activate
Cells.Select
Selection.Clear
Range("a1").Select
For i = 1 To n
ki = ki1: Call step(x0, y0, del, x1, y1): x0 = x1: x0s = x0
For j = 1 To m: y0(j) = y1(j): y0s(j) = y0(j): Next j
If (i \ kv) = (i / kv) Then
If ki2 = 1 Then
If ipr1 = 1 Then
Worksheets("RESULTS").Cells(ipr1, 1) = "x0"
Worksheets("RESULTS").Cells(ipr1, 2) = "y0(1)"
Worksheets("RESULTS").Cells(ipr1, 3) = "y0(2)": ipr1 = ipr1 + 1
End If
Worksheets("RESULTS").Cells(ipr1 + 1, 1) = x0
Worksheets("RESULTS").Cells(ipr1 + 1, 2) = y0(1)
Worksheets("RESULTS").Cells(ipr1 + 1, 3) = y0(2)
ipr1 = ipr1 + 1
End If
If ki2 = 2 Then ki = ki2: Call dydx(x0, y0, pr)
End If
Next i
End Sub
Результаты моделирования:
1) Номинальный режим
|
K1-6 |
0,01 |
|
|
|
|
|
|
|
|
|
|
|
x0 |
y0(1) |
y0(2) |
|
|
|
|
|
|
|
|
|
|
0 |
5 |
5 |
|
|
|
|
|
|
|
|
|
|
100 |
8,460128 |
7,170963 |
|
|
|
|
|
|
|
|
|
|
200 |
9,382276 |
8,657288 |
|
|
|
|
|
|
|
|
|
|
300 |
9,392796 |
8,831817 |
|
|
|
|
|
|
|
|
|
|
400 |
9,393432 |
8,839948 |
|
|
|
|
|
|
|
|
|
|
500 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
600 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
700 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
800 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
900 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
1000 |
9,393452 |
8,840228 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2) Режим заполнения
|
K4-6 |
0 |
|
|
|
|
|
|
|
|
|
|
|
x0 |
y0(1) |
y0(2) |
|
|
|
|
|
|
|
|
|
|
0 |
0 |
0 |
|
|
|
|
|
|
|
|
|
|
100 |
3,861676 |
0,197603 |
|
|
|
|
|
|
|
|
|
|
200 |
7,317577 |
0,626418 |
|
|
|
|
|
|
|
|
|
|
300 |
9,324503 |
1,509201 |
|
|
|
|
|
|
|
|
|
|
400 |
9,345662 |
2,71921 |
|
|
|
|
|
|
|
|
|
|
500 |
9,347763 |
3,919028 |
|
|
|
|
|
|
|
|
|
|
600 |
9,350527 |
5,10254 |
|
|
|
|
|
|
|
|
|
|
700 |
9,35456 |
6,263083 |
|
|
|
|
|
|
|
|
|
|
800 |
9,361382 |
7,386741 |
|
|
|
|
|
|
|
|
|
|
900 |
9,376054 |
8,436975 |
|
|
|
|
|
|
|
|
|
|
1000 |
9,421958 |
9,273752 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
3) Режим аварийного слива
|
K1-3 |
0 |
|
|
|
|
|
|
|
|
|
|
|||||
|
x0 |
y0(1) |
y0(2) |
|
|
|
|
|
|
|
|
||||||
|
0 |
5 |
5 |
|
|
|
|
|
|
|
|
||||||
|
100 |
4,846792 |
3,746256 |
|
|
|
|
|
|
|
|
||||||
|
200 |
4,608491 |
2,714879 |
|
|
|
|
|
|
|
|
||||||
|
300 |
4,337529 |
1,818757 |
|
|
|
|
|
|
|
|
||||||
|
400 |
4,049469 |
1,023702 |
|
|
|
|
|
|
|
|
||||||
|
500 |
3,751773 |
0,310429 |
|
|
|
|
|
|
|
|
||||||
|
600 |
3,448661 |
|
|
|
|
|
|
|
|
|
||||||
|
700 |
3,142777 |
|
|
|
|
|
|
|
|
|
||||||
|
800 |
2,835903 |
|
|
|
|
|
|
|
|
|
||||||
|
900 |
2,529316 |
|
|
|
|
|
|
|
|
|
||||||
|
1000 |
2,223969 |
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
||||||
