Постановка задачи
Ранее (см., например, в заметках Пчела в ячейке, Пчела на соте, Тепловое взаимодействие ячеек в соте и др.) мы моделировали теплообмен отдельной пчелы или нескольких пчел с сотом и воздухом. При этом воздух улочки (среды между сотами) задавался как бесконечный резервуар с постоянной температурой T_ambient, например, 12°С.
При моделировании улочки это допущение уже не подойдет: воздух улочки имеет конечный объём и теплоемкость; его температура зависит от общего теплообмена. Он нагревается от сот и пчел в улочке и охлаждается за счет вентиляции (конвекции) и теплопроводности (от стенок улья).
В данной заметке мы сделаем лишь первый шаг: расчёт стационарного поля температур в зазоре между двумя сотами без пчёл в улочке и без вертикальной конвекции воздуха. Источники тепла: только входы ячеек обоих сотов, нагретые пчёлами внутри ячеек. Понятно, что такой ситуации не бывает, но полученный "демон" позволит нам увидеть условия для вертикальной конвекции в возникшем поле температур, и в дальнейшем построить модель этой конвекции.
Поскольку улочка образуется двумя сотами (стороной Y одного сота и стороной ^ - другого), мы сгенерируем фрагмент одного сота 11×11 ячеек, а второй зададим точно такой же. Тогда вторую сторону (сторону ^ другого сота) можно взять из сгенерированного первого, а не считать отдельно - для существа дела это ничего не портит.
Сетка узлов улочки
Ширина улочки составляет 8-12 мм (а при наличии пчел будет еще меньше), поэтому для описания состояния (в частности - температуры в разных местах улочки) ограничимся одной плоскостью, проходящей вертикально точно по середине улочки параллельно сотам, т.е. это будет двумерная сетка.
При моделировании сота мы ограничились гексагональной сеткой, соответствующей центрам входов в ячейки - нам известны температуры в каждой из этих точек, как на стороне Y, так и на стороне ^. Естественным представляется и для сетки узлов улочки оставить тот же вид (топологию): каждой точке на соте будет противостоять узел улочки. Но как же тогда для противоположной стороны ^ - ведь там, как мы помним, ячейки с соответствующими индексами смещены относительно стороны Y на пол ячейки?
Поскольку период узлов по вертикали и горизонтали у них одинаковый, то можно пренебречь этим расхождением, и считать, например, что узлы сетки улочки находятся посредине отрезка соединяющего центры ячеек с одинаковыми индексами строки и столбца на противоположных сторонах - тепловая картина в улочке от этого не изменится, но это избавит нас от необходимости создавать дополнительные сетки узлов в улочке. Вообще говоря, для тепловой модели узел - это не геометрическая точка, а дискрета, в пределах которой мы считаем все параметры одинаковыми. В соте, например, для пустой ячейки - это воздух в ячейке, а также любая из девяти граней ячейки - 10 узлов; при наличии пчелы в ячейке добавится еще три узла - голова, торакс и брюшко. В рассматриваемой пустой улочке - это один узел: воздух, объемом = (площадь шестиугольника ячейки × ширину улочки). Это лишь подчеркивает допустимость принятого подхода.
Физика узла
Считаем, что каждый узел улочки взаимодействует (обменивается теплом) только с ближайшими узлами сот (с такими же номерами строки и столбца) и шестью соседними узлами улочки (по гексагональной структуре). С ячейками сот происходит конвективный обмен (как это мы учитывали и раньше с коэффициентом конвекции h_entr = 11 Вт/м²K):
$$ Q_{entr}=h_{entr}\cdot A_{entr} \cdot(T_{air}^Y−T_{street})+h_{entr}\cdot A_{entr} \cdot(T_{air}^∧−T_{street}) $$ а с соседними узлами - теплопроводностью, т.к. конвекцию в улочке пока не учитываем:
$$ Q_{lat}=∑_{n=1}^6G_{lat}(T_{street}^{(n)}−T_{street}) $$ Уравнение теплового баланса для узла:
$$ C_{node}\frac{dT_{street}}{dt}=Q_{entr}+Q_{lat} $$ Примечание: строго говоря, обмен теплопроводностью с ячейками сот также есть, но он значительно меньше конвективной составляющей (перемешивания воздуха на входе в ячейку).
Разместим по четыре пчелы на каждом соте, как в заметке Тепловое взаимодействие ячеек в соте . Программа для такой модели прикреплена в конце заметки. Если ее запустить, то увидим вот такое распределение температур в узлах улочки при ширине улочки 8, 10 и 12 мм:
Видим, что действительно воздух в узлах улочки прогрелся и прогрелся неравномерно. Это создает условия для возникновения конвекции.
Данную заметку публикую отдельно, чтобы зафиксировать именно топологию узлов сетки улочки. Она далась не сразу, т.к. хотелось по аналогии с пакетами вычислительной гидродинамики задать прямоугольную, а не гексагональную сетку, т.к. конвекция возникает именно вертикальная. Предполагалось также, в этом случае задать три слоя сетки - посредине улочки (как сейчас), и по краям - ближе к сотам. Они, естественно, были бы смещены друг относительно друга по горизонтали и по вертикали. В общем, получалась довольно сложная структура, которая была бы трудна для восприятия и громоздка для последующих вычислений. А так - получился еще один ряд "виртуальных" гексагональных ячеек, с которыми уже привычно работать. И, как видите, там можно рассчитывать температуру в улочке.
🐝 Python-скрипт
"""
fast_street_empty.py — Модель улочки БЕЗ пчёл в улочке.
Улочка = зазор между двумя одинаковыми сотами:
Сот 1: сторона Y (к улочке) — пчёлы в ячейках греют
Сот 2: сторона ^ (к улочке) — пчёл нет, греется через дно
Задача: рассчитать поле T_street(col, row) — температуру воздуха в улочке.
Физика: 2D теплопроводность воздуха + обмен с входами ячеек обоих сотов.
Параметр: d_street (ширина улочки, 8–12 мм).
Верификация:
1. Нет пчёл в сотах → T_street = T_ambient везде.
2. G_lat → ∞ → T_street = среднее(T_air_Y, T_air_X) везде.
3. d_street → 0 → каждый узел независимо: T_street = (T_air_Y + T_air_X)/2.
"""
import time
import math
import numpy as np
from numba import njit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
# ==========================================================
# 1. ПАРАМЕТРЫ (из fast_surface_bee.py, без поверхностной пчелы)
# ==========================================================
cols, rows = 11, 11
variant = 2
T_ambient = 12.0
# --- Пчела ---
d_bee = 0.004
s_r = math.pi * (d_bee / 2) ** 2
h_head, h_thorax, h_ab = 0.0015, 0.0045, 0.006
h_total = h_head + h_thorax + h_ab
A_hd = math.pi * d_bee * h_head + s_r
A_th = math.pi * d_bee * h_thorax
A_ab = math.pi * d_bee * h_ab + s_r
V_bee = h_total * s_r
ke0 = 3e-5
Pmin = 0.2e-3
k_hd, k_ab = 0.065, 0.042
m_head, m_thorax, m_ab = 10.25e-6, 32.5e-6, 57e-6
c_bee = 3500
C_hd, C_th, C_ab = m_head * c_bee, m_thorax * c_bee, m_ab * c_bee
# --- Воздух ---
rho_air, k_air, c_air = 1.205, 0.025, 1005
# --- Воск ---
L_w = [0.0001, 0.00014, 0.0002]
L_b = [0.00013, 0.00055, 0.0012]
c_w = [3300., 3000., 2700.]
lam_w = [0.235, 0.1608, 0.055]
eps_w = [0.22, 0.5, 0.9]
L_wall = L_w[variant - 1]
L_bottom = L_b[variant - 1]
c_wax = c_w[variant - 1]
lam_wax = lam_w[variant - 1]
rho_wax = 960
# --- Геометрия ячейки ---
thread = 5.4e-3
d_cell = thread - L_wall
h_cell = 12e-3
dh = 1e-3
a = d_cell / 2. / np.cos(np.radians(30))
A_lateral = a * (h_cell - dh / 2.)
d1, d2 = d_cell, np.sqrt(a ** 2 + dh ** 2)
A_bottom = d1 * d2 / 2.
C_wax_lateral = rho_wax * A_lateral * L_wall * c_wax
C_wax_bottom = rho_wax * A_bottom * L_bottom * c_wax
A_wall_total = 6 * A_lateral + 3 * A_bottom
area_hex = (3 * math.sqrt(3) / 2) * a ** 2
A_entrance = area_hex
V_cell = area_hex * h_total
C_air_cell = rho_air * V_cell * c_air
C_air_bee = rho_air * (V_cell - V_bee) * c_air
# --- Излучение ---
sigma = 5.67e-8
t0 = -273.15
eps = 0.22
eps_wall = eps_w[variant - 1]
se = sigma * eps_wall
se_th, se_hd, se_ab_ = se * A_th, se * (A_hd - s_r), se * (A_ab - s_r)
# --- Коэффициенты конвекции ---
h_in = k_air / ((d_cell - d_bee) / 2)
h_entrance = 11.0
# --- Сопротивления (пчела в ячейке) ---
R_th_ai = 1 / (h_in * A_th)
R_hd_ai = 1 / (h_in * A_hd)
R_ab_ai = 1 / (h_in * (A_ab - s_r))
R_ab_ = 1 / (h_entrance * s_r)
hA_ = h_entrance * (A_entrance - s_r)
hA = h_entrance * A_entrance
R_head_air_open = 1 / (h_entrance * A_hd)
R_thorax_air_open = 1 / (h_entrance * A_th)
R_ab_air_open = 1 / (h_entrance * A_ab)
P_warm_up = 0.0
# ==========================================================
# 2. ТОПОЛОГИЯ СОТА (без поверхностной пчелы)
# ==========================================================
CELL_EMPTY = 0
CELL_BEE = 5
I_HD, I_TH, I_AB, I_AI = 0, 1, 2, 3
R_ai_wall_arr = np.array([L_wall / (k_air * A) for A in [A_lateral] * 6 + [A_bottom] * 3], dtype=np.float64)
R_wall_self = np.array([L_wall / (lam_wax * A_lateral)] * 6 + [L_bottom / (lam_wax * A_bottom)] * 3, dtype=np.float64)
C_wall_base = np.array([C_wax_lateral] * 6 + [C_wax_bottom] * 3, dtype=np.float64)
cells = []
for side in ["Y", "^"]:
for col in range(-int(cols / 2), int(cols / 2) + 1):
for row in range(-int(rows / 2), int(rows / 2) + 1):
cells.append((0, side, col, row))
cell_to_index = {cell: i for i, cell in enumerate(cells)}
N_CELLS = len(cells)
EVEN_ROW_DELTAS = [(-1, 0), (-1, -1), (0, -1), (1, 0), (0, 1), (-1, 1)]
ODD_ROW_DELTAS = [(1, 0), (1, -1), (0, -1), (-1, 0), (0, 1), (1, 1)]
def get_neighbor_for_face(side, col, row, face_id):
if face_id < 6:
deltas = ODD_ROW_DELTAS if (row % 2 == 1) else EVEN_ROW_DELTAS
if side == "^":
deltas = EVEN_ROW_DELTAS if (row % 2 == 1) else ODD_ROW_DELTAS
dc, dr = deltas[face_id]
return (0, side, int(col + dc), int(row + dr))
else:
k = face_id - 6
if side == "Y":
bottom_list = [(col, row - 1), (col, row), (col - 1, row)] if row % 2 == 0 else \
[(col, row - 1), (col, row), (col + 1, row)]
n_side = "^"
else:
bottom_list = [(col, row), (col, row + 1), (col + 1, row)] if row % 2 == 0 else \
[(col, row), (col, row + 1), (col - 1, row)]
n_side = "Y"
n_col, n_row = bottom_list[k]
return (0, n_side, n_col, n_row)
Topology = np.full((N_CELLS, 9), -1, dtype=np.int32)
Neighbor_Faces = np.full((N_CELLS, 9), -1, dtype=np.int32)
all_cells_set = set(cells)
for i, cell in enumerate(cells):
_, side, col, row = cell
for face_id in range(9):
n_addr = get_neighbor_for_face(side, col, row, face_id)
if n_addr in all_cells_set:
Topology[i, face_id] = cell_to_index[n_addr]
_, n_side, n_col, n_row = n_addr
for rev_f in range(9):
if get_neighbor_for_face(n_side, n_col, n_row, rev_f) == (0, side, col, row):
Neighbor_Faces[i, face_id] = rev_f
break
# --- Пчёлы в ячейках (только сторона Y) ---
States = np.full(N_CELLS, CELL_EMPTY, dtype=np.int32)
bee_cells = [(0, 'Y', 0, 0), (0, 'Y', 0, 1), (0, 'Y', 0, 2), (0, 'Y', 2, -1)]
bee_idx = [cell_to_index[cell] for cell in bee_cells]
for idx in bee_idx:
States[idx] = CELL_BEE
# Cover = 0 везде (нет поверхностной пчелы)
Cover = np.zeros(N_CELLS, dtype=np.float64)
print(f">>> Сот: cols={cols}, rows={rows}, пчёл в ячейках={len(bee_cells)}, поверхностных=0")
# ==========================================================
# 3. РАСЧЁТ СОСТОЯНИЯ СОТА (600 с, без поверхностной пчелы)
# ==========================================================
def calculate_initial_temperatures(T_air_val):
l_val = (1.3493 * np.exp(-0.058 * T_air_val)) * 1e-3
R_th_hd_val = l_val / (k_hd * s_r)
R_th_ab_val = l_val / (k_ab * s_r)
def equations(vars_val):
Thd_v, Tth_v, Tab_v = vars_val
q_evap = ke0 * Tth_v
P_rest = max((29.133 - 0.739 * T_air_val) * 1e-3, Pmin)
p_s_hd = -eps * sigma * ((Thd_v - t0) ** 4 - (T_air_val - t0) ** 4) * A_hd
p_s_th = -eps * sigma * ((Tth_v - t0) ** 4 - (T_air_val - t0) ** 4) * A_th
p_s_ab = -eps * sigma * ((Tab_v - t0) ** 4 - (T_air_val - t0) ** 4) * A_ab
eq1 = p_s_hd + (Tth_v - Thd_v) / R_th_hd_val - (Thd_v - T_air_val) / R_head_air_open - q_evap
eq2 = P_rest + p_s_th - (Tth_v - Thd_v) / R_th_hd_val - (Tth_v - Tab_v) / R_th_ab_val - \
(Tth_v - T_air_val) / R_thorax_air_open
eq3 = p_s_ab + (Tth_v - Tab_v) / R_th_ab_val - (Tab_v - T_air_val) / R_ab_air_open - q_evap
return [eq1, eq2, eq3]
return fsolve(equations, [T_air_val + 10, T_air_val + 15, T_air_val + 10], xtol=1e-8)
Thd_init, Tth_init, Tab_init = calculate_initial_temperatures(T_ambient)
print(f"Начальные T пчелы: Hd={Thd_init:.2f}, Th={Tth_init:.2f}, Ab={Tab_init:.2f}")
# --- Матрицы для Numba ---
C_grid = np.zeros((N_CELLS, 13), dtype=np.float64)
for i in range(N_CELLS):
C_grid[i, 4:13] = C_wall_base
if States[i] == CELL_BEE:
C_grid[i, I_HD] = C_hd
C_grid[i, I_TH] = C_th
C_grid[i, I_AB] = C_ab
C_grid[i, I_AI] = C_air_bee
else:
C_grid[i, I_AI] = C_air_cell
T_grid_A = np.full((N_CELLS, 13), T_ambient, dtype=np.float64)
T_grid_B = np.full((N_CELLS, 13), T_ambient, dtype=np.float64)
for t_idx_b in bee_idx:
T_grid_A[t_idx_b, I_HD] = Thd_init
T_grid_A[t_idx_b, I_TH] = Tth_init
T_grid_A[t_idx_b, I_AB] = Tab_init
T_grid_B[t_idx_b, I_HD] = Thd_init
T_grid_B[t_idx_b, I_TH] = Tth_init
T_grid_B[t_idx_b, I_AB] = Tab_init
# --- Ядро симуляции сота (без поверхностной пчелы) ---
@njit(fastmath=True, cache=True)
def do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid,
dt, T_ambient_val, heater_idx):
"""Один шаг тепловой сети сота. Без поверхностной пчелы (Cover=0)."""
Q_entrance_Y = 0.0
Q_entrance_X = 0.0
half_cells = len(T_cur) // 2
for i in range(len(T_cur)):
state = States[i]
Thd = T_cur[i, I_HD]
Tth = T_cur[i, I_TH]
Tab = T_cur[i, I_AB]
# --- T воздуха в ячейке ---
if state == CELL_BEE:
num = Tth / R_th_ai + Thd / R_hd_ai + Tab / R_ab_ai + T_ambient_val / R_ab_
den = 1.0 / R_th_ai + 1.0 / R_hd_ai + 1.0 / R_ab_ai + 1.0 / R_ab_
else:
num = T_ambient_val * hA
den = hA
for f in range(9):
num += T_cur[i, 4 + f] / R_ai_wall_arr[f]
den += 1.0 / R_ai_wall_arr[f]
Tai = num / den
T_next[i, I_AI] = Tai
is_side_Y = (i < half_cells)
# --- Теплопотери через вход ---
if state == CELL_BEE:
Q_cell = hA_ * (Tai - T_ambient_val)
Q_ab_loss = (Tab - T_ambient_val) / R_ab_
if is_side_Y:
Q_entrance_Y += (Q_cell + Q_ab_loss)
else:
Q_entrance_X += (Q_cell + Q_ab_loss)
else:
Q_cell = hA * (Tai - T_ambient_val)
if is_side_Y:
Q_entrance_Y += Q_cell
else:
Q_entrance_X += Q_cell
# --- Стенки ---
T_wall_avg = 0.0
for f in range(9):
T_wall_avg += T_cur[i, 4 + f]
T_wall_avg /= 9.0
p_s_th = -se_th * ((Tth - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_hd = -se_hd * ((Thd - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_ab = -se_ab_ * ((Tab - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_total = p_s_th + p_s_hd + p_s_ab
p_rad_to_walls = p_s_total if state == CELL_BEE else 0.0
for f in range(9):
Q_from_air = (Tai - T_cur[i, 4 + f]) / R_ai_wall_arr[f]
neighbor_idx = Topology[i, f]
if neighbor_idx != -1:
n_face = Neighbor_Faces[i, f]
Q_from_neighbor = (T_cur[neighbor_idx, 4 + n_face] - T_cur[i, 4 + f]) / R_wall_self[f]
else:
Q_from_neighbor = (T_ambient_val - T_cur[i, 4 + f]) / R_wall_self[f]
p_s_k = ((A_lateral if f < 6 else A_bottom) / A_wall_total) * p_rad_to_walls
T_next[i, 4 + f] = T_cur[i, 4 + f] + (dt / C_grid[i, 4 + f]) * (
Q_from_air + Q_from_neighbor + p_s_k)
# --- Пчела в ячейке ---
if state == CELL_BEE:
P_rest = (29.133 - 0.739 * Tai) * 1e-3
if P_rest < Pmin:
P_rest = Pmin
P_local = P_warm_up if i == heater_idx else 0.0
l_dyn = (1.3493 * np.exp(-0.058 * Tai)) * 1e-3
R_th_hd_dyn = l_dyn / (k_hd * s_r)
R_th_ab_dyn = l_dyn / (k_ab * s_r)
T_next[i, I_TH] = Tth + (dt / C_grid[i, I_TH]) * (
P_rest + P_local + p_s_th
- (Tth - Thd) / R_th_hd_dyn - (Tth - Tab) / R_th_ab_dyn - (Tth - Tai) / R_th_ai)
T_next[i, I_HD] = Thd + (dt / C_grid[i, I_HD]) * (
(Tth - Thd) / R_th_hd_dyn - (Thd - Tai) / R_hd_ai - ke0 * Tth)
T_next[i, I_AB] = Tab + (dt / C_grid[i, I_AB]) * (
p_s_ab + (Tth - Tab) / R_th_ab_dyn - (Tab - Tai) / R_ab_ai
- (Tab - T_ambient_val) / R_ab_ - ke0 * Tth)
return Q_entrance_Y, Q_entrance_X
# --- Запуск расчёта сота ---
dt_sim = 0.1
time_total = 600.0
steps_comb = int(time_total / dt_sim)
# Прогрев Numba
do_cluster_step(T_grid_A, T_grid_B, Topology, Neighbor_Faces, States, C_grid,
dt_sim, T_ambient, bee_idx[0])
print(f"Расчёт сота ({steps_comb} шагов)...")
start_time = time.time()
T_cur, T_next = T_grid_A, T_grid_B
for step in range(steps_comb):
Q_Y, Q_X = do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid,
dt_sim, T_ambient, bee_idx[0])
T_cur, T_next = T_next, T_cur
print(f" Сот рассчитан за {time.time() - start_time:.2f} с. "
f"Q_Y={Q_Y * 1000:.3f} мВт, Q_X={Q_X * 1000:.3f} мВт")
# ==========================================================
# 4. ИЗВЛЕЧЕНИЕ T_air ДЛЯ УЛОЧКИ
# ==========================================================
# Сот 1 (сторона Y → к улочке): T_air_Y[col, row]
# Сот 2 (сторона ^ → к улочке): T_air_X[col, row]
col_range = range(-int(cols / 2), int(cols / 2) + 1)
row_range = range(-int(rows / 2), int(rows / 2) + 1)
n_cols = len(col_range)
n_rows = len(row_range)
T_air_Y = np.full((n_cols, n_rows), T_ambient)
T_air_X = np.full((n_cols, n_rows), T_ambient)
for i, cell in enumerate(cells):
_, side, col, row = cell
ci = col + int(cols / 2)
ri = row + int(rows / 2)
if side == "Y":
T_air_Y[ci, ri] = T_cur[i, I_AI]
else:
T_air_X[ci, ri] = T_cur[i, I_AI]
print(f"\nT_air_Y (сторона Y, к улочке):")
print(f" min={T_air_Y.min():.3f}, max={T_air_Y.max():.3f}, mean={T_air_Y.mean():.3f} °C")
print(f"T_air_X (сторона ^, к улочке):")
print(f" min={T_air_X.min():.3f}, max={T_air_X.max():.3f}, mean={T_air_X.mean():.3f} °C")
# ==========================================================
# 5. МОДЕЛЬ УЛОЧКИ: 2D ТЕПЛОПРОВОДНОСТЬ
# ==========================================================
# --- Параметры улочки ---
d_street_values = [8e-3, 10e-3, 12e-3] # Ширина улочки (м)
# Гексагональная сетка: те же смещения, что у ячеек
# Для стороны Y: чётные ряды → EVEN_ROW_DELTAS, нечётные → ODD_ROW_DELTAS
# (уличка наследует топологию стороны Y сота 1)
STREET_EVEN_DELTAS = [(-1, 0), (-1, -1), (0, -1), (1, 0), (0, 1), (-1, 1)]
STREET_ODD_DELTAS = [(1, 0), (1, -1), (0, -1), (-1, 0), (0, 1), (1, 1)]
def get_street_neighbors(ci, ri, n_cols, n_rows):
"""Соседи узла (ci, ri) в гексагональной сетке улочки.
ci, ri — индексы в массиве [0..n_cols-1, 0..n_rows-1].
row = ri - int(rows/2) — для определения чётности."""
row = ri - int(rows / 2)
deltas = STREET_ODD_DELTAS if (row % 2 == 1) else STREET_EVEN_DELTAS
neighbors = []
for dc, dr in deltas:
ni, nj = ci + dc, ri + dr
if 0 ≤ ni < n_cols and 0 ≤ nj < n_rows:
neighbors.append((ni, nj))
return neighbors
# Предвычисление соседей для всех узлов
street_neighbors = [[None] * n_rows for _ in range(n_cols)]
for ci in range(n_cols):
for ri in range(n_rows):
street_neighbors[ci][ri] = get_street_neighbors(ci, ri, n_cols, n_rows)
# --- Граничные условия: периметр = T_ambient (Дирихле) ---
is_boundary = np.zeros((n_cols, n_rows), dtype=bool)
for ci in range(n_cols):
for ri in range(n_rows):
if ci == 0 or ci == n_cols - 1 or ri == 0 or ri == n_rows - 1:
is_boundary[ci, ri] = True
print(f"\nУличка: сетка {n_cols}×{n_rows}, "
f"внутренних узлов={np.sum(~is_boundary)}, граничных={np.sum(is_boundary)}")
# --- Решатель улочки (явный Эйлер, 2D теплопроводность) ---
def solve_street(T_air_Y, T_air_X, d_street, dt_street=0.1, n_steps=200):
"""
Решить 2D теплопроводность в улочке.
Уравнение узла (col, row):
C_node * dT/dt = G_entr*(T_air_Y - T) + G_entr*(T_air_X - T)
+ Σ G_lat*(T_neighbor - T)
Граничные условия: T = T_ambient на периметре.
"""
# Параметры узла
A_node = (math.sqrt(3) / 2) * thread ** 2 # Площадь контрольного объёма (м²)
V_node = A_node * d_street # Объём узла (м³)
C_node = rho_air * V_node * c_air # Теплоёмкость узла (Дж/К)
G_entr = h_entrance * A_entrance # Проводимость входа (Вт/К)
G_lat = k_air * d_street / math.sqrt(3) # Проводимость между соседями (Вт/К)
# Начальное условие: T_ambient везде
T_street = np.full((n_cols, n_rows), T_ambient)
# Логирование
log_T_center = [] # T в центре (над пчёлами)
log_T_mean = [] # Средняя T по внутренним узлам
log_steps = []
for step in range(n_steps):
T_new = T_street.copy()
for ci in range(n_cols):
for ri in range(n_rows):
if is_boundary[ci, ri]:
T_new[ci, ri] = T_ambient # Дирихле
continue
T_here = T_street[ci, ri]
# Источник: обмен с входами ячеек обоих сотов
Q_source = G_entr * (T_air_Y[ci, ri] - T_here) + \
G_entr * (T_air_X[ci, ri] - T_here)
# Латеральная теплопроводность (6 соседей)
Q_lat = 0.0
for (ni, nj) in street_neighbors[ci][ri]:
Q_lat += G_lat * (T_street[ni, nj] - T_here)
# Явный Эйлер
T_new[ci, ri] = T_here + (dt_street / C_node) * (Q_source + Q_lat)
T_street = T_new
# Логирование
if step % 10 == 0 or step == n_steps - 1:
log_steps.append(step)
# Центр сетки (col=0, row=0) → индекс [3, 3] для 7×7
ci_c, ri_c = int(cols / 2), int(rows / 2)
log_T_center.append(T_street[ci_c, ri_c])
# Средняя по внутренним узлам
interior_mask = ~is_boundary
log_T_mean.append(T_street[interior_mask].mean())
return T_street, log_steps, log_T_center, log_T_mean, C_node, G_entr, G_lat
# --- Расчёт для трёх значений d_street ---
results = {}
for d_street in d_street_values:
d_mm = d_street * 1000
print(f"\n d_street = {d_mm:.0f} мм...")
T_street, log_steps, log_T_center, log_T_mean, C_node, G_entr, G_lat = \
solve_street(T_air_Y, T_air_X, d_street, dt_street=0.05, n_steps=400)
results[d_mm] = {
'T_street': T_street,
'log_steps': log_steps,
'log_T_center': log_T_center,
'log_T_mean': log_T_mean,
'C_node': C_node,
'G_entr': G_entr,
'G_lat': G_lat,
}
print(f" C_node={C_node * 1e4:.3f}×10⁻⁴ Дж/К, "
f"G_entr={G_entr * 1e4:.3f}×10⁻⁴ Вт/К, "
f"G_lat={G_lat * 1e4:.3f}×10⁻⁴ Вт/К")
print(f" T_center={T_street[int(cols / 2), int(rows / 2)]:.3f}°C, "
f"T_mean_interior={T_street[~is_boundary].mean():.3f}°C")
# ==========================================================
# 6. ВЕРИФИКАЦИЯ
# ==========================================================
print(f"\n{'=' * 60}")
print(" ВЕРИФИКАЦИЯ")
print(f"{'=' * 60}")
# Проверка 1: стационар (dT/dt ≈ 0)
for d_mm, res in results.items():
T = res['T_street']
# Ещё 10 шагов — изменение должно быть < 1e-6
T_check = T.copy()
for _ in range(10):
T_new = T_check.copy()
for ci in range(n_cols):
for ri in range(n_rows):
if is_boundary[ci, ri]:
continue
T_here = T_check[ci, ri]
Q = G_entr * (T_air_Y[ci, ri] - T_here) + G_entr * (T_air_X[ci, ri] - T_here)
for (ni, nj) in street_neighbors[ci][ri]:
Q += res['G_lat'] * (T_check[ni, nj] - T_here)
T_new[ci, ri] = T_here + (0.05 / res['C_node']) * Q
T_check = T_new
max_change = np.max(np.abs(T_check - T))
print(f" d={d_mm:.0f}мм: макс. изменение за 10 доп. шагов = {max_change:.2e} °C "
f"{'✓ стационар' if max_change < 1e-6 else '✗ НЕ стационар'}")
# Проверка 2: без пчёл → T_street = T_ambient
T_air_Y_cold = np.full((n_cols, n_rows), T_ambient)
T_air_X_cold = np.full((n_cols, n_rows), T_ambient)
T_cold, _, _, _, _, _, _ = solve_street(T_air_Y_cold, T_air_X_cold, 10e-3, n_steps=50)
max_dev_cold = np.max(np.abs(T_cold - T_ambient))
print(f" Без пчёл: макс. отклонение от T_amb = {max_dev_cold:.2e} °C "
f"{'✓' if max_dev_cold < 1e-10 else '✗'}")
# Проверка 3: G_lat → ∞ → T_street = (T_air_Y + T_air_X)/2 (для внутренних)
# (При бесконечной латеральной проводимости все узлы выравниваются)
T_avg_expected = (T_air_Y + T_air_X) / 2
T_avg_interior = T_avg_expected[~is_boundary].mean()
print(f" G_lat→∞: ожидаемая T (среднее по внутренним) = {T_avg_interior:.4f} °C")
# ==========================================================
# 7. ВИЗУАЛИЗАЦИЯ
# ==========================================================
# --- График 1: Динамика установления T_street ---
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)
colors = {8: 'blue', 10: 'green', 12: 'red'}
for d_mm, res in results.items():
t_arr = np.array(res['log_steps']) * 0.05 # время в секундах
ax1.plot(t_arr, res['log_T_center'], color=colors[d_mm], linewidth=2,
label=f'd={d_mm:.0f} мм (центр)')
ax2.plot(t_arr, res['log_T_mean'], color=colors[d_mm], linewidth=2,
label=f'd={d_mm:.0f} мм (среднее)')
ax1.axhline(T_ambient, color='gray', linestyle=':', label=f'T_amb={T_ambient}°C')
ax1.set_ylabel('T_street (°C)')
ax1.set_title('Установление T_street: центр сетки (над пчёлами)')
ax1.legend(fontsize=9);
ax1.grid(True, linestyle=':')
ax2.axhline(T_ambient, color='gray', linestyle=':')
ax2.set_xlabel('Время (с)')
ax2.set_ylabel('T_street (°C)')
ax2.set_title('Установление T_street: среднее по внутренним узлам')
ax2.legend(fontsize=9);
ax2.grid(True, linestyle=':')
plt.tight_layout();
plt.show()
# --- График 2: 2D поля T_street для трёх d_street ---
fig, axes = plt.subplots(1, 4, figsize=(12, 4), gridspec_kw={'width_ratios': [1, 1, 1, 0.08]})
ax_cbar = axes[3]
# Общий диапазон температур
all_T = np.concatenate([res['T_street'][~is_boundary].ravel() for res in results.values()])
vmin, vmax = all_T.min() - 0.1, all_T.max() + 0.1
for ax, (d_mm, res) in zip(axes, results.items()):
T = res['T_street']
# Гексагональная отрисовка
size = 1.0
width_hex = math.sqrt(3) * size
dy_hex = 1.5 * size
im_data = []
for ci in range(n_cols):
for ri in range(n_rows):
col = ci - int(cols / 2)
row = ri - int(rows / 2)
# Позиция центра гексагона (гексагональная сетка)
x = width_hex * (col + 0.5 * (row & 1))
y = -dy_hex * row
# Гексагон
corners_x, corners_y = [], []
for k in range(6):
angle = math.pi / 180 * (60 * k - 30)
corners_x.append(x + size * math.cos(angle))
corners_y.append(y + size * math.sin(angle))
corners_x.append(corners_x[0])
corners_y.append(corners_y[0])
norm_T = (T[ci, ri] - vmin) / max(vmax - vmin, 1e-5)
color = plt.cm.inferno(norm_T)
if is_boundary[ci, ri]:
ax.fill(corners_x, corners_y, color='lightgray',
edgecolor='gray', linewidth=0.8, alpha=0.5)
else:
ax.fill(corners_x, corners_y, color=color,
edgecolor='black', linewidth=0.8)
ax.text(x, y, f"{T[ci, ri]:.2f}", ha='center', va='center',
fontsize=7, color='white' if norm_T > 0.5 else 'black')
ax.set_title(f'd_street = {d_mm:.0f} мм\n'
f'T_center={T[int(cols / 2), int(rows / 2)]:.2f}°C, '
f'T_mean={T[~is_boundary].mean():.2f}°C', fontsize=10)
ax.set_aspect('equal')
ax.axis('off')
sm = plt.cm.ScalarMappable(cmap=plt.cm.inferno, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array([])
sm = plt.cm.ScalarMappable(cmap=plt.cm.inferno,
norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array([])
cb = fig.colorbar(sm, cax=ax_cbar)
cb.set_label('T_street (°C)')
axes[3].set_title('')
plt.suptitle('Поле температур в улочке без пчёл\n', fontsize=11)
plt.tight_layout()
plt.show()
# --- График 3: Сечение через центр (row=0) ---
fig, ax = plt.subplots(figsize=(8, 5))
row_center = int(rows / 2) # row=0
col_arr = np.array([c - int(cols / 2) for c in range(n_cols)])
for d_mm, res in results.items():
T_row = res['T_street'][:, row_center]
ax.plot(col_arr, T_row, 'o-', color=colors[d_mm], linewidth=2,
markersize=6, label=f'd={d_mm:.0f} мм')
# T_air источников для сравнения
ax.plot(col_arr, T_air_Y[:, row_center], 's--', color='orange', linewidth=1.5,
markersize=5, alpha=0.7, label='T_air_Y (сот 1, сторона Y)')
ax.plot(col_arr, T_air_X[:, row_center], '^--', color='purple', linewidth=1.5,
markersize=5, alpha=0.7, label='T_air_X (сот 2, сторона ^)')
ax.axhline(T_ambient, color='gray', linestyle=':', label=f'T_amb={T_ambient}°C')
ax.set_xlabel('col (положение в соте)')
ax.set_ylabel('Температура (°C)')
ax.set_title(f'Сечение улочки через центр (row=0).\n'
f'Источники: T_air_Y (пчёлы) и T_air_X (пусто)')
ax.legend(fontsize=9);
ax.grid(True, linestyle=':')
plt.tight_layout();
plt.show()
# ==========================================================
# 8. ИТОГОВЫЙ ВЫВОД
# ==========================================================
print(f"\n{'=' * 60}")
print(f" РЕЗУЛЬТАТЫ: УЛОЧКА БЕЗ ПЧЁЛ (d_street = 8, 10, 12 мм)")
print(f"{'=' * 60}")
print(f"\n Источники тепла (входы ячеек):")
print(f" Сот 1 (Y, пчёлы): T_air = {T_air_Y.min():.2f}...{T_air_Y.max():.2f}°C")
print(f" Сот 2 (^, пусто): T_air = {T_air_X.min():.2f}...{T_air_X.max():.2f}°C")
for d_mm, res in results.items():
T = res['T_street']
interior = T[~is_boundary]
print(f"\n d_street = {d_mm:.0f} мм:")
print(f" C_node = {res['C_node'] * 1e4:.3f}×10⁻⁴ Дж/К")
print(f" G_entr = {res['G_entr'] * 1e4:.3f}×10⁻⁴ Вт/К (один вход)")
print(f" G_lat = {res['G_lat'] * 1e4:.3f}×10⁻⁴ Вт/К (один сосед)")
print(f" T_street: min={interior.min():.3f}, max={interior.max():.3f}, "
f"mean={interior.mean():.3f}°C")
print(f" Градиент (центр→край): "
f"{(T[int(cols / 2), int(rows / 2)] - T[1, int(rows / 2)]) / ((int(cols / 2) - 1) * thread * 100):.2f} °C/см")
print(f"\n ФИЗИКА:")
print(f" Улочка асимметрична: сот 1 (Y, пчёлы) греет сильнее,")
print(f" сот 2 (^, пусто) — слабее (греется через дно от стороны Y).")
print(f" Теплопроводность воздуха (k=0.025 Вт/м·К) сглаживает градиент.")
print(f" Ширина улочки влияет на G_lat ∝ d_street:")
print(f" шире → сильнее латеральная связь → поле более однородное.")
print(f"{'=' * 60}")
Комментарии
Отправить комментарий