Моделирование системы водоснабжения

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

# Откроем модель
modelName = "WaterSupply_New"
if !(modelName in [m.name for m in engee.get_all_models()]) engee.load( "$(@__DIR__)/$(modelName).engee"); end;

Рассчитаем параметры системы в зависимости от числа домохозяйств

# =================================================================
# 1. КОНФИГУРАТОР ПОТРЕБИТЕЛЕЙ (ДОМОХОЗЯЙСТВ)
# =================================================================
# ВВЕДИТЕ КОЛИЧЕСТВО ДОМОХОЗЯЙСТВ:
num_households = 175  # При 150-200 домохозяйствах параметры совпадут с прошлым расчетом

# Эмпирические зависимости проектирования по СНиП / СП 31.13330:
# Требуемый рабочий объем башни (примерно 7-10% от суточного потребления поселка)
V_tower = round(num_households * 0.085, digits=1) 

# Необходимая высота башни для обеспечения напора у потребителей
# (с ростом поселка и длины труб требуемая высота башни увеличивается)
H_tower = 10.0 + (num_households - 175) * 0.05 
Reservoir_IN_H = H_tower + 2 # Точка излива напорной трубы насоса

# Глубина всасывания (из скважины/резервуара до насоса)
H_suction = 5.0  

println("--- Параметры системы для $num_households домохозяйств ---")
println("Расчетный объем водонапорной башни:      ", V_tower, " м³")
println("Геометрическая высота башни (нагнетание):", round(H_tower, digits=1), " м")
println("Глубина всасывания (динамический уровень):", H_suction, " м")

Зададим параметры насоса

# =================================================================
# 2. РАСЧЕТ И ПОДБОР ПАРАМЕТРОВ НАСОСА
# =================================================================
# Статический напор системы (полный геометрический перепад высот)
H_static = H_tower + H_suction

# Динамические потери в трубах (примем стандартные ~33% от статики при номинальном расходе)
H_losses = H_static * 0.3333
H_nom = H_static + H_losses # Полный требуемый номинальный напор насоса

# Номинальная подача насоса (рассчитывается так, чтобы заполнить башню за ~15 минут)
Q_nom_m3h = round((V_tower / 15.0) * 60.0, digits=1)

# Формирование граничных точек для Engee (параболическая кривая Q-H)
H_shutoff = round(H_nom * 1.25, digits=1)       # Напор при закрытой задвижке (+50%)
Q_max_m3h = round(Q_nom_m3h / sqrt(1.0 - (H_nom / H_shutoff)), digits=2)  # Подача при нулевом напоре (+41.6%)

# Перевод подачи насоса в систему СИ (м³/с)
Q_nom = Q_nom_m3h / 3600.0
Q_max = Q_max_m3h / 3600.0

Подберем электродвигатель

# =================================================================
# 3. РАСЧЕТ И ПОДБОР МОЩНОСТИ ДВИГАТЕЛЯ ДЛЯ НАСОСА
# =================================================================
# Гидравлические константы и КПД
eta_pump  = 0.72      # Номинальный КПД насоса (72%)
rho       = 1000.0    # Плотность перекачиваемой жидкости (вода), кг/м³
g         = 9.81      # Ускорение свободного падения, м/с²

P_hydr = rho * g * Q_nom * H_nom 
P_shaft = P_hydr / eta_pump
k_zap = 1.15
P_calc_kw = (P_shaft * k_zap) / 1000.0

# Автоматический выбор ближайшего стандартного двигателя
standard_powers = [0.55, 0.75, 1.1, 1.5, 2.2, 3.0, 4.0, 5.5, 7.5, 11.0, 15.0, 18.5, 22.0, 30.0]
P_n = minimum(filter(x -> x >= P_calc_kw, standard_powers)) * 1000.0

# =================================================================
# 3. ВХОДНЫЕ КАТАЛОЖНЫЕ ДАННЫЕ ДЛЯ ВЫБРАННОГО ДВИГАТЕЛЯ (5.5 кВт)
# =================================================================
U_line  = 400.0; f = 50.0; n_n = 2920.0; eta_mot = 0.857; cos_phi = 0.88
k_i     = 7.0; p = 1; J_phys = 0.014; F_phys = 0.0011

# =================================================================
# 4. ВЫЧИСЛЕНИЕ БАЗОВЫХ ЭЛЕКТРИЧЕСКИХ ВЕЛИЧИН И СХЕМЫ ЗАМЕЩЕНИЯ
# =================================================================
U_ph   = U_line / sqrt(3)
S_n    = P_n / (eta_mot * cos_phi)
I_n    = S_n / (sqrt(3) * U_line)
Z_base = U_ph / I_n
omega_base = 2 * pi * f
L_base = Z_base / omega_base

R_s = 0.0435 * Z_base
n_synch = (60.0 * f) / p
s_n     = (n_synch - n_n) / n_synch
R_r_prime = s_n * Z_base * cos_phi
X_k = Z_base / k_i; L_k = X_k / omega_base
L_ls = 0.4 * L_k; L_lr_prime = 0.6 * L_k
I_mu = 0.35 * I_n; X_m = (U_ph / I_mu) - (omega_base * L_ls); L_m = X_m / omega_base

# Перевод в p.u. (относительные единицы)
R_s_pu = R_s / Z_base; R_r_pu = R_r_prime / Z_base
L_ls_pu = L_ls / L_base; L_lr_pu = L_lr_prime / L_base; L_m_pu = L_m / L_base; L_0_pu = L_ls_pu
omega_synch = (2 * pi * n_synch) / 60.0
H_inertia = (J_phys * omega_synch^2) / (2 * P_n)
M_base = P_n / (omega_synch * (1 - s_n))
F_pu = (F_phys * omega_synch) / M_base

println("--- Параметры схемы замещения  блока электродвигателя (p.u.) ---")
println("Rs: ", round(R_s_pu, digits=5), " | Rr': ", round(R_r_pu, digits=5))
println("Lls: ", round(L_ls_pu, digits=5), " | Llr': ", round(L_lr_pu, digits=5), " | Lm: ", round(L_m_pu, digits=5))
println("Inertia H: ", round(H_inertia, digits=5), " s")

Построим характеристику H-Q насоса

# =================================================================
# 4. МАТЕМАТИЧЕСКАЯ АППРОКСИМАЦИЯ КРИВОЙ Q-H ДЛЯ ГРАФИКА
# =================================================================
# Вычисляем коэффициент 'a' строго по точке нулевого напора (на излив)
a_coeff = H_shutoff / (Q_max_m3h^2)

# Функция напора от расхода
H_pump(Q) = max(0.0, H_shutoff - a_coeff * Q^2) # max не дает напору уйти в минус

# Генерация массива точек для графика
Q_vector = collect(0.0:0.5:Q_max_m3h)
H_vector = H_pump.(Q_vector)


# =================================================================
# 5. АВТОМАТИЧЕСКОЕ ПОСТРОЕНИЕ ГРАФИКА Q-H
# =================================================================
plot(Q_vector, H_vector, 
     label="Кривая Q-H насоса", 
     linewidth=3, 
     color=:blue,
     title="Напорно-расходная характеристика насоса",
     xlabel="Подача Q, м³/ч", 
     ylabel="Напор H, м",
     grid=true,
     xlims=(0, Q_max_m3h + 5),
     ylims=(0, H_shutoff + 5))

# Добавление номинальной рабочей точки (маркер на графике)
scatter!([Q_nom_m3h], [H_nom], 
         color=:red, 
         markersize=7, 
         label="Рабочая точка")

# Добавление точки закрытой задвижки
scatter!([0.0], [H_shutoff], 
         color=:darkgreen, 
         markersize=5, 
         label="Закрытая задвижка")

# Добавление точки работы на излив
scatter!([Q_max_m3h], [0.0], 
         color=:orange, 
         markersize=5, 
         label="Точка на излив")

Смасштабируем объем башни для анализа суточного потребления за 240 секунд

using Plots

println("=================================================================")
println("  МАСШТАБИРОВАНИЕ ОБЪЕМА БАШНИ ДЛЯ УСКОРЕНИЯ КОЛЕБАНИЙ (240 с)   ")
println("=================================================================\n")

# --- Исходная конфигурация поселка --- (175 домов по умолчанию)


# Среднесуточное потребление на 1 дом ≈ 900 литров (0.9 м³)
V_day_total_liters = num_households * 900.0 

# Коэффициент сжатия времени (1440 минут суток / 240 секунд симуляции)
time_compression = 1440.0 / 240.0 # = 6.0

# =================================================================
# 6. УПРАВЛЯЕМЫЙ ФЛАГ МАСШТАБИРОВАНИЯ И ЕДИНЫЕ ИМЕНА ПЕРЕМЕННЫХ БАКА
# =================================================================
# true -> для теста суточных колебаний за 240с (масштабированные параметры)
# false -> реальный физический номинал (реальные геометрические параметры)
apply_scaling = true  

# Расчет масштабированной геометрии
V_tower_scaled  = round(V_tower / time_compression, digits=2)
Area_real       = V_tower / H_tower
Area_scaled     = V_tower_scaled / H_tower

# --- ЕДИНЫЕ НАИМЕНОВАНИЯ ПЕРЕМЕННЫХ ДЛЯ МАСКИ МОДЕЛИ ---
# (Прописанные параметры в блоке бака)
V_tank    = apply_scaling ? V_tower_scaled : V_tower
Area_tank = apply_scaling ? Area_scaled : Area_real

println("--- Результаты масштабирования накопительной емкости ---")
println("Текущий режим: ", apply_scaling ? "МАСШТАБИРОВАННЫЙ (Симуляция 240с)" : "ФИЗИЧЕСКИЙ (Номинал)")
println("Активное значение объема V_tank:        ", V_tank, " м³")
println("Активное значение площади Area_tank:   ", round(Area_tank, digits=3), " м²")
println("Высота башни H_tower (константа):        ", H_tower, " м")

# =================================================================
# 6. ГЕНЕРАЦИЯ РЕАЛЬНОГО СУТОЧНОГО ГРАФИКА (1440 точек)
# =================================================================
V_day_total_liters = num_households * 900.0 
hourly_coeffs = [
    0.35, 0.25, 0.20, 0.15, 0.25, 0.50, # 00:00 - 06:00
    1.20, 2.10, 1.80, 1.20, 1.00, 0.95, # 06:00 - 12:00 (Утренний пик)
    1.10, 1.30, 1.00, 0.90, 1.10, 1.70, # 12:00 - 18:00
    2.30, 2.00, 1.50, 1.10, 0.70, 0.45  # 18:00 - 24:00 (Вечерний пик)
]
hourly_coeffs = hourly_coeffs ./ sum(hourly_coeffs)

Q_24h_lmin = zeros(1440)
for m in 0:1439
    hour_idx = div(m, 60) + 1
    Q_24h_lmin[m+1] = (V_day_total_liters * hourly_coeffs[hour_idx]) / 60.0
end

# =================================================================
# 7. МАСШТАБИРОВАНИЕ ВРЕМЕНИ НА ШКАЛУ 240 С С ШАГОМ 1 C
# =================================================================
time_sim_sec = collect(0.0:1.0:240.0) 
Q_sim_lmin   = zeros(241)

for t in 0:240
    corresponding_minute = t * time_compression
    
    # Безопасное определение индексов с ограничением диапазона [1, 1440]
    idx1 = clamp(floor(Int, corresponding_minute) + 1, 1, 1440)
    idx2 = clamp(ceil(Int, corresponding_minute) + 1, 1, 1440)
    
    weight = corresponding_minute - floor(corresponding_minute)
    
    # Сохраняем честные физические расходы потребителей (л/мин)
    Q_sim_lmin[t+1] = (1.0 - weight) * Q_24h_lmin[idx1] + weight * Q_24h_lmin[idx2]
end


println("\n--- РЕЗУЛЬТАТ МАСШТАБИРОВАНИЯ ВЕКТОРОВ ДЛЯ LOOKUP TABLE ---")
println("Длина вектора времени 'time_sim_sec': ", length(time_sim_sec), " элементов")
println("Длина вектора расходов 'Q_sim_lmin':   ", length(Q_sim_lmin), " элементов")
println("Пиковое потребление на графике:        ", round(maximum(Q_sim_lmin), digits=1), " л/мин")

Запуск модели

# =================================================================
# 8. Запуск модели
# =================================================================

data = engee.run( modelName )

Проверка выходных параметров

# Так как время симуляции достаточно большое (240 с), строим графики с прореживанием

# 1. Задаем параметры времени
step_sec = 1.0 
t_max = 240.0

# Массив секунд для оси X (0, 1, 2 ... 240)
target_times = 0.0:step_sec:t_max

# 2. Находим общее количество строк в исходных данных 
# Загружаем один массив, чтобы узнать, сколько всего шагов сделала симуляция
raw_throttle = collect(WorkspaceArray{Float64}("$(modelName)/Throttle"))
total_rows = size(raw_throttle, 1)

# Вычисляем индексы строк, которые равномерно распределены от 1 до конца таблицы на 241 точку
indices = round.(Int, range(1, total_rows, length=length(target_times)))

# 3. Загружаем данные из WorkspaceArray и фильтруем их по индексам
df_throttle = raw_throttle[indices, :]
df_power    = collect(WorkspaceArray{Float64}("$(modelName)/Активная мощность (о.е.)"))[indices, :]
df_speed    = collect(WorkspaceArray{Float64}("$(modelName)/Скорость ротора (о.е.)"))[indices, :]
df_q1       = collect(WorkspaceArray{Float64}("$(modelName)/Расходомер-1/Q1"))[indices, :]
df_q2       = collect(WorkspaceArray{Float64}("$(modelName)/Расходомер/Q2"))[indices, :]
df_fluid    = collect(WorkspaceArray{Float64}("$(modelName)/FluidLevel"))[indices, :]
df_current = collect(WorkspaceArray{Vector{Float64}}("$(modelName)/iABC/Ток по фазам"))

raw_current_time = df_current.time


# 4. Строим 4 графика с привязкой к оси времени в секундах
p1 = plot( target_times, [df_throttle.value .*1e6], title="Площадь отверстия истечения (открытие задвижки) (м²)", ylims=(0, 1000))

p2 = plot( raw_current_time, reduce(hcat, df_current.value)', title="Ток в обмотках (А)", xlims=(160.9, 161.0), ylims=(-20, 20))

p3 = plot( target_times, [df_power.value])
     plot!( target_times, [df_speed.value], title="Активная мощность и скорость ротора (отн.ед.)" )

p4 = plot( target_times, [df_q1.value])
     plot!( target_times, [df_q2.value], title="Расход потребителя и расход насоса (л/мин)", ylims=(0.0,1000) )

p5 = plot( target_times, [df_fluid.value], title="Уровень воды в башне", xlabel="Время, с" )

# 5. Сборка финального полотна
plot( p1, p2, p3, p4, p5, layout=(5,1), leg=false, size=(1200,600), titlefont=font(10) )

Заключение

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

Система водоснабжения с водонапорной башней

Моделирование системы водоснабжения

Рассчитаем параметры системы в зависимости от числа домохозяйств

Зададим параметры насоса

Подберем электродвигатель

Построим характеристику H-Q насоса

Смасштабируем объем башни для анализа суточного потребления за 240 секунд

Запуск модели

Проверка выходных параметров

Заключение

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник