РХТУ им. Д.И. Менделеева

Кафедра информатики и компьютерного моделирования

Практическая работа №1

по курсу

«Компьютерное моделирование химических процессов»

Компьютерное моделирование простой гидравлической системы в стационарном режиме

Группа: О-41

Студент: Faz Birmutt :)

Вариант: 3

Преподаватель: Тамбовцев И.И.

Москва 2006

Схема гидравлической системы:

Известные величины:

Давления на входе и выходе системы: P₁ - P₆;

Коэффициенты сужающих устройств: k₁ - k₇;

Геометрические высоты емкостей: H₁, H₂.

Неизвестные величины:

Высоты столбов жидкости в емкостях: h₁, h₂;

Давления жидкости на дно сосуда: P₇, P₈;

Давления газа над жидкостью: P₉, P₁₀;

Расходы жидкости v₁ - v₇.

Основные допущения, накладываемые на модель:

1) Система работает в стационарном режиме;

2) Рассматривается только гидравлика процесса;

3) Сопротивление трения потоку пренебрежимо мало;

4) Жидкость несжимаема и имеет постоянную плотность;

5) Газ - идеальный pV=nRT;

6) Режим - изотермический T=const;

7) Трубопроводы находятся на одном уровне;

8) Цилиндрическая форма закрытой емкости с площадью поперечного сечения

S и геометрической высотой H;

9) Одинаковое давление газа P_н в емкостях, не заполненных жидкостью;

Система уравнений математического описания:

A) Уравнения расходов жидкости через клапаны:

1. v₁ = k₁ (p₁ - p₇)^1/2

2. v₂ = k₂ (p₂ - p₇) ^1/2

3. v₃ = k₃ (p₃ - p₇) ^1/2

4. v₄ = k₄ (p₈ - p₄) ^1/2

5. v₅ = k₅ (p₈ - p₅) ^1/2

6. v₆ = k₆ (p₈ - p₆) ^1/2

7. v₇ = k₇ (p₇ - p₈) ^1/2

B) Уравнения балансов емкостей:

8. v₁ + v₂ + v₃ - v₇ = 0

9. v₇ - v₄ - v₅ - v₆ = 0

C) Определение давлений жидкости и газа в емкостях:

10. p₇ = p₉ + ρ g h₁

11. p₉ = P_н H₁ / (H₁ - h₁)

12. p₈ = p₁₀ + ρ g h₂

13. p₁₀ = P_н H₂ / (H₂ - h₂)

Информационная матрица системы уравнений:

N п/п	v1	v2	v3	v4	v5	v6	v7	p7	p8	p9	p10	h1	h2
1
2
3
4
5
6
7
8
9
10
11
12
13

Блок-схема алгоритма расчета:

Текст программы:

Option Explicit

Option Base 1

Const np% = 10, nk% = 7, nv% = 7

Dim vm!(nk), v!(nk), ak!(nk), p!(np), hg!(2), h!(2)

Dim a!, b!, c!, e!, ro!, pn!, g!, x!

Dim i%, kl%, ipr%

Dim bu As Boolean

Public Sub stat()

ipr = 1

With Worksheets("TASK")

hg(1) = .Cells(4, 5): hg(2) = .Cells(5, 5) ' Высота емкостей (1-2) м

ro = .Cells(6, 5) ' Плотность кг/м3

pn = .Cells(6, 9) ' Начальное давление МПа

For i = 1 To 6: p(i) = .Cells(8, i + 4): Next i ' Давление (1-6) МПа

For i = 1 To nk: ak(i) = .Cells(9, i + 4): Next i ' Коэф. пpопускной способности (1-7)

e = .Cells(11, 6) ' Относительная локальная погpешность %

kl = .Cells(10, 6) ' вывод промежеточных результатов 0-нет, 1-частичный, 2-полный

End With

Worksheets("RESULTS").Activate

Cells.Select

Selection.Clear

Range("a1").Select

If kl = 2 Then

Cells(ipr, 5) = "Промежуточный вывод": ipr = ipr + 1

Cells(ipr, 5) = "h": Cells(ipr, 6) = "p(5-7)": Cells(ipr, 7) = "vm": ipr = ipr + 1

End If

g = 9.815: e = e / 100: a = 0: b = hg(1) * (1 - e)

Call MPD(a, b, e, bu, x)

With Worksheets("RESULTS")

If bu Then

' здесь подставляем свои коэффициенты

a = ro * g * 0.000001: b = p(8) + ro * g * hg(2) * 0.000001: c = (p(8) - pn) * hg(2) ''11'

h(2) = (b - Sqr(b * b - 4 * a * c)) / 2 / a

p(10) = pn * hg(2) / (hg(2) - h(2)) ''10

For i = 1 To nk: vm(i) = v(i) * ro: Next i

.Cells(1, 1) = "РЕЗУЛЬТАТ"

.Cells(2, 1) = "h": .Cells(2, 2) = "p(7-10)": .Cells(2, 3) = "vm(1-7)"

.Cells(3, 1) = h(1): .Cells(4, 1) = h(2)

For i = 1 To 4: .Cells(2 + i, 2) = p(6 + i): Next i

For i = 1 To nv: .Cells(2 + i, 3) = vm(i): Next i

'MsgBox ("V1 + V2 + V3 = " + Str(vm(1) + vm(2) + vm(3)) + Chr(13) + "V4 + V5 + V6 = " + Str(vm(4) + vm(5) + vm(6)) + Chr(13) + "V7 = " + Str(vm(7))) ' проверка на равенство

Else

kl = 2

.Cells(1, 1) = "РЕШЕНИЯ НЕТ"

.Cells(2, 1) = "a": .Cells(2, 2) = "f(a)": .Cells(2, 3) = "b": .Cells(2, 4) = "f(b)"

.Cells(3, 1) = a

.Cells(ipr, 5) = "Промежуточный вывод a": ipr = ipr + 1

.Cells(ipr, 5) = "h": .Cells(ipr, 6) = "p(5-7)": .Cells(ipr, 7) = "vm": ipr = ipr + 1

.Cells(3, 2) = FUNC(a)

.Cells(3, 3) = b

.Cells(1, 5) = "Промежуточный вывод b": ipr = ipr + 1

.Cells(2, 5) = "h": .Cells(2, 6) = "p(5-7)": .Cells(2, 7) = "vm": ipr = ipr + 1

.Cells(3, 4) = FUNC(b)

End If

End With

End Sub

Function FUNC(x!) As Single

Dim vm!(7), fx!

' здесь подставляем свои уравнения

h(1) = x

p(9) = pn * hg(1) / (hg(1) - h(1)) ''8

p(7) = p(9) + ro * g * h(1) * 0.000001 ''9

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(7) = v(1) + v(2) + v(3) ''12

p(8) = p(7) - Sgn(v(7)) * (v(7) / ak(7)) ^ 2 ''7

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

fx = (v(7) - v(4) - v(5) - v(6)) * ro ''13

For i = 1 To nk: vm(i) = v(i) * ro: Next i

If kl = 0 Then GoTo 400

If kl = 1 Then GoTo 300

Cells(ipr, 5) = h(1): Cells(ipr, 6) = p(5): Cells(ipr, 7) = vm(1): ipr = ipr + 1

Cells(ipr, 6) = p(6): Cells(ipr, 7) = vm(2): ipr = ipr + 1

Cells(ipr, 6) = p(7): Cells(ipr, 7) = vm(3): ipr = ipr + 1

Cells(ipr, 7) = vm(4): ipr = ipr + 1

Cells(ipr, 7) = vm(5): ipr = ipr + 1

300: Cells(ipr, 5) = "x = ": Cells(ipr, 6) = x: Cells(ipr, 7) = " fx = ": Cells(ipr, 8) = fx: ipr = ipr + 1

400: FUNC = fx

End Function

Sub MPD(a!, b!, eps!, bu As Boolean, xcon!)

Dim fa!, fb!, x!, fx!

fa = FUNC(a): fb = FUNC(b)

If fa * fb > 0 Then: bu = False: GoTo 100

Do

x = (a + b) / 2: fx = FUNC(x)

If fx * fa < 0 Then b = x Else a = x

Loop While Abs(a - b) > eps

xcon = Abs(a + b) / 2: bu = True

100:

End Sub

Sub auto_open()

Worksheets("TASK").Activate

End Sub

Результаты моделирования (h1-2, p7-10, v1-7):

1) Чувствительность к изменению высоты емкостей

H1	5
4,710637	1,774169	-5,236114
9,038226	1,128456	4,752169
	1,727934	8,519572
	1,039745	-3,486321
		7,927519
		3,584071
		8,035628
H1	10
9,405434	1,774212	-5,236530
9,038432	1,128681	4,751710
	1,681898	8,519316
	1,039968	-3,483083
		7,928942
		3,587218
		8,034495
H1	15
14,083122	1,774213	-5,236532
9,038433	1,128682	4,751708
	1,635987	8,519315
	1,039969	-3,483070
		7,928949
		3,587231
		8,034491