Уравнения связного графа. Решение системы

В этом видео получена система. Решим эти уравнения связного графа, применив Python, в частности библиотеки numpy, scipy и matplotlib. Построим графики решения.

Скачать только скрипт.

Скачать целиком виртуальное окружение Python (Jupiter Notebook и все библиотеки).

Имя архива виртуального окружения: HydraulicEquationsSolve01.zip
Размер: 195841054 байтов (186 MiB)
CRC32: 836DFC91

Исходный код в виртуальном окружении приложен на Яндекс-диске. Это сделано для удобства. Можно просто скачать, запустить виртуальное окружение и редактировать код в Jupyter Notebook.

О решении системы уравнений связного графа

На всякий случай, напомним. Система состоит из сопротивлений, двух баков и одной длинной трубы, у которой нельзя пренебречь инерцией.

О неточностях. У трубы ‘lf’ должно быть сопротивление. В задаче учитывается только инерция. Более в задаче трубопроводов нет. Соединения между баками и дросселями — идеальные. Повторю здесь картинку со схемой из видео.

Система решается путем создания функции, вычисляющей производные. То есть, во-первых надо записать уравнения системы соответственно, а во-вторых, надо сделать функцию, которая их вычисляет и скормить эту функцию решателю вместе с начальными условиями. Да, и не следует забывать и о том, что одно из уравнений алгебраическое.

Кстати, решатель меняется очень легко. Попробуйте поменять BDF на RK45 в вызове solve_ivp( ). Получится шум в давлении. Это из-за того, что давление определяется дифференцированием момента Г. Говорят, что это т.н. жёсткая система. С решателем BDF картинка более гладкая.

Код — Расчет констант

Здесь, согласно параметрам трубопровода, баков и дросселей мы считаем сопротивления, емкости и индуктивность. В соответствии с этой статьей.

Значения выбраны так, чтобы расход через второй дроссель был минимальным. Его диаметр 1 мм. Вся рабочая жидкость уходит в баки. Это позволяет провести небольшую проверку. Сумма объемов в баках должна быть равна объему подачи от источника. Давление в каждом баке определяется высотой столба жидкости. Давление в первом падает на первом дросселе R1 но не всё, а до уровня давления во втором баке. На дросселе R2 должно падать все давление из второго бака.

import numpy as np
import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import csv

# ------------------ Define constants --------------------
mu = 1e-3      # Pa·s (dynamic viscosity of water)
rho = 1000     # kg/m^3 (density of water)
g = 9.81       # m/s^2 (gravitational acceleration)

# Tank, pipe, throttles dimensions
D1 = 0.03      # m (diameter of the long pipe)
Dthr1 = 0.03   # m (diameter of throttle R1)
Dthr2 = 0.001  # m (diameter of throttle R2)
Dt1 = 0.5      # m (C1 tank diameter)
Dt2 = 0.5      # m (C2 tank diameter)
L = 200        # m (length of the long pipe)

# Flow rate in m^3/s from the source.
q = 100 / 60 / 1000             # 100 liters/min converted to m^3/s

# Cross-sectional areas of the tanks.
A_tank_1 = np.pi * Dt1**2 / 4   # Area of C1
A_tank_2 = np.pi * Dt2**2 / 4   # Area of C2
print('{0}, {1}'.format(A_tank_1, A_tank_2))

# Calculate orifice resistances.
R1 = (128 * mu * L) / (np.pi * Dthr1**4)
R2 = (128 * mu * L) / (np.pi * Dthr2**4)

# Calculate capacities of the tanks.
C1 = A_tank_1 / (rho * g)
C2 = A_tank_2 / (rho * g)

# Calculate hydraulic inertia of the liquid in the pipe.
ALP = np.pi * D1**2 / 4
lf = rho * L / ALP

Код — Решение системы уравнений

Система уравнений решается численно. Используется встроенная функция библиотеки SciPy, которая называется solve_ivp( ). Подход здесь стандартный для библиотечных численных методов. Система выражается в виде функции, подходящей для алгоритма. А именно, функция дает значения производных в зависимости от значений искомых величин.

Величин три. Объем в первом баке, во втором и гидравлический момент в длинной трубке, которая имеет инерцию.

# Create Pandas dataframe to store results
data = pd.DataFrame(0, index=range(1), columns=['Time', 'PC1', 'PR1', 'QR1', 'PC2', 'Qlf', 'PR2'])

# Define the system of DAEs
def hydraulic_system(t, y):
    global data
    V1, V2, G = y  # V1 is volume in C1, V2 is volume in C2, G is Gamma

    # Ptimary conditions (PC) in junctions.
    PC1 = (V1 / C1)             # PC in 0-junction for C1.
    PR1 = (V1 / C1) - (V2 / C2) # Algebraic equation for PR1 (over throttle R1).
    QR1 = PR1 / R1              # PC in 1-junction for R2.
    PC2 = (V2 / C2)             # PC in 0-junction for C2.
    Qlf = G / lf                # PC in 1-junction for R2 and pipe.
    
    PR2 = Qlf * R2              # Calculation for PR2 drop.

    
    # Save the calculated PC and PR1 to the dataframe.
    new_row = pd.DataFrame({'Time': [t], 'PC1': [PC1], 'PR1': [PR1], 
                            'QR1': [QR1], 'PC2': [PC2], 'Qlf': [Qlf], 'PR2': [PR2]})
    if new_row.notna().all(axis=None):
        data = pd.concat([data, new_row], ignore_index=True)
    
    # Differential equations
    V1_dot = q - PR1 / R1
    V2_dot = PR1 / R1 - G / lf
    G_dot = V2 / C2 - G * R2 / lf
    
    return [V1_dot, V2_dot, G_dot]

# Initial conditions: (initial volumes V1, V2, and G)
V1_0 = 0  # initial volume in tank C1
V2_0 = 0  # initial volume in tank C2
G_0 = 0   # initial Gamma (pressure in pipe)

# Time span for the solution
process_time = 60  # seconds
t_span = (0, process_time)  # time span in seconds
y0 = [V1_0, V2_0, G_0]

# Solve the system using solve_ivp
sol = solve_ivp(hydraulic_system, t_span, y0, method='BDF', dense_output=True)

Постпроцессинг некоторых величин

# Time points for plotting
t = np.linspace(0, process_time, 1000)
y = sol.sol(t)

# Postprocessing Gamma
# Calculate Gamma' (numerical derivative of Gamma)
Gamma = y[2]
Gamma_prime = np.gradient(Gamma, t)  # Numerical derivative of Gamma
# Calculate pressure difference
pressure_difference = Gamma_prime

Графики результата

# Create a figure and a set of subplots
fig, axes = plt.subplots(7, 1, figsize=(10, 25))

# Plot V1 (Volume in C1)
axes[0].plot(t, y[0], label='V1 (Volume in C1)', color='b')
axes[0].set_ylabel('Volume')
axes[0].set_title('V1 (Volume in C1)')
axes[0].legend()
axes[0].grid(True)

# Plot V2 (Volume in C2)
axes[1].plot(t, y[1], label='V2 (Volume in C2)', color='g')
axes[1].set_ylabel('Volume')
axes[1].set_title('V2 (Volume in C2)')
axes[1].legend()
axes[1].grid(True)

# Plot Gamma (Pressure in pipe)
axes[2].plot(t, pressure_difference, label='dPi (Pressure diff in pipe)', color='r')
axes[2].set_ylabel('Pressure')
axes[2].set_title('dPi (Pressure diff in pipe)')
axes[2].legend()
axes[2].grid(True)

# Plot calculated parameter PR1.
axes[3].plot(data['Time'], data['PR1'], label='PR1 (delta over R1 throttle)', color='b')
axes[3].set_ylabel('Pressure drop')
axes[3].set_title('Pressure drop PR1')
axes[3].legend()
axes[3].grid(True)

# Plot calculated parameter PR2.
axes[4].plot(data['Time'], data['PR2'], label='PR2 (delta over R2 throttle)', color='r')
axes[4].set_ylabel('Pressure drop')
axes[4].set_title('Pressure drop PR2')
axes[4].legend()
axes[4].grid(True)

# Plot QR1
axes[5].plot(data['Time'], data['QR1'], label='QR1', color='b')
axes[5].set_ylabel('QR1')
axes[5].set_title('QR1 flow rate')
axes[5].legend()
axes[5].grid(True)

# Plot Qlf
axes[6].plot(data['Time'], data['Qlf'], label='Qlf', color='g')
axes[6].set_ylabel('Qlf')
axes[6].set_title('Qlf flow rate through pipe and R2')
axes[6].legend()
axes[6].grid(True)

Проверки

# Print check results
print('Check:')
print('Total volume by source = {0}(m^3).'.format(q*process_time))
print('V1 = {0}(m^3) V2 = {1}(m^3), sum = {2}'.format(y[0][-1], y[1][-1], y[0][-1] + y[1][-1]))
print('Pressure in the C2: from liquid level {0}, pressure in the pipe {1}'.format(rho*g*(y[1][-1]/A_tank_2), data['PC2'].iloc[-1]))

Повторюсь, в самом начале, мы установили диаметр дросселя R2 1 мм. Это сделано, чтобы вся РЖ пошла в баки.
Поэтому, по сумме объемов в баке и имея заданный расход источника, можно проверить расходы. Проверка по конечному значению (на 60-й секунде). Давление проверяется иначе, нужно считать гидростатическое давление баков. Оно падает на дросселях и трубе.

Давление в первом баке определяется высотой столба РЖ. Для приведенных графиков (t=60c),

$V_1 = 0,08; A_1 = 0,196$ , тогда $h_1 = 0,4$ и эта высота столба соответствует давлению $\rho*g*h = 4000$ Па.

Аналогично из уровня во втором баке можно получить высоту 0,061 м, которая соответствует 600 Па. Примерно эти цифры и падают на графиках Pressure Drop PR1 и Pressure Drop PR2, там они рассчитаны из расходов.