После того, как сформирована модель пчелы и сота (рамки), следующим шагом представляется модель улочки. Что такое улочка - это пространство между двумя сотами (или сотом и стенкой улья, колоды), в котором на одном и другом соте (стороне улочки) могут находиться пчелы, а между ними - воздушный зазор, ширина которого зависит от ширины улочки. Каждая пчела вырабатывает тепло (модель термогенеза у нас уже есть) и обменивается им с ячейками сот, с другими пчелами и с воздухом в улочке. Но математической модели этого обмена - пока нет, а без него, очевидно, модели улочки не построить.
Поэтому модель улочки начинаем разрабатывать с модели пчелы на соте.
Поскольку пчела будет нагревать сот, а статус ячеек под пчелой может быть любым, в том числе там может находиться расплод или пчела, постольку нужна тепловая модель сота (она у нас также есть), и модель пчелы на соте явится по сути ее развитием и продолжением.
Напомню, в качестве геометрической модели пчелы были приняты три соосных цилиндра диаметром 4 мм и высотой, соответственно, 1.5 (голова), 4.5 (торакс) и 6 мм (брюшко) с небольшими зазорами, зависящими от температуры окружающего воздуха. При расположении на соте пчела может характеризоваться своим положением и ориентацией. Ограничим множество положений пчелы множеством совпадений центра торакса с центром ячейки, а возможных ориентаций оставим 6 - по направлениям, перпендикулярным боковым стенкам ячеек. Полагаю, что для тепловой модели этого множества будет достаточно. Таким образом, если торакс находится в центре ячейки, то площадь его проекции (18 мм²) будет перекрывать 75% площади входа в ячейку (24 мм²), над соседней ячейкой голова будет перекрывать 25%, а над другой соседней ячейкой брюшко будет перекрывать практически 100% площади входа. (Пишу "практически", т.к. диаметр описанной окружности ячейки немного больше диаметра пчелы, но брюшко пчелы "выглядывает" за пределы ячейки и перекрывает небольшую часть четвертой ячейки, поэтому интегральный эффект экранирования будет примерно, как при перекрытии одной ячейки на 100%, т.к. площадь проекции брюшка равна площади входа в ячейку - 24 мм²). Эти проценты будут учтены при расчете теплообмена.
Рассмотрим небольшой кластер (7×7 ячеек) с пчелой в центральной ячейке (Y, 0, 0), как это уже делали ранее, но теперь разместим на стороне Y еще одну пчелу в строке №1 на соте. Адрес ее торакса будет (Y, 1, 0), головы - (Y, 1, -1), брюшка - (Y, 1, 1). Нижняя (вентральная) поверхность пчелы обменивается теплом с ячейками, а верхняя (дорсальная, спина) - с улочкой. Тогда модель теплообмена торакса можно представить так:
C_th · dT_th/dt = P_met(T_sense) ← генерация
− (T_th − T_hd)/R_th_hd ← теплопроводность → голова
− (T_th − T_ab)/R_th_ab ← теплопроводность → брюшко
− (T_th − T_amb)/R_th_amb_surf ← дорсальная конвекция
− (T_th − T_air)/R_th_cell_surf ← вентральная конвекция
+ P_rad_th_dorsal(T_th, T_amb) ← излучение в улочку
+ P_rad_th_ventral(T_th, T_wall) ← излучение в ячейки
Модели для головы и брюшка устроены аналогично, но без P_met и с добавлением испарения q_evap = ke0·T_th.
Мощность метаболизма зависит от температуры, которую пчела измеряет усиками, т.е. от температуры в районе головы. Для пчелы в ячейке это была температура воздуха в этой ячейке, а для пчелы на соте, очевидно, это будет температура в ячейке под головой, в нашем случае (Y, 1, -1).
Полный текст Python-скрипта - в спойлере внизу заметки. Если его запустить, то получим такие результаты:
Как видим на верхнем графике, картина температур в центральной ячейке с пчелой практически не изменилась - там она определяется прежде всего именно пчелой в ячейке, а пчела на соте слабо на нее влияет. На нижнем графике показано, как меняются температуры пчелы на соте. Температуры пчелы в ячейке и пчелы на соте существенно различаются: пчела в ячейке холоднее, т.к. больше отдает тепла сотам, стенки ячейки работают как радиатор охлаждения, или иными словами - пчела в ячейке греет сот эффективнее, чем пчела на соте.
При этом, т.к. температуры в ячейке с пчелой и в ячейке под головой пчелы на соте разные, то и мощность термогенеза у них разная:
у пчелы в ячейке она составляет 19.28 мВт, а на соте - 19.86 мВт.
На следующем рисунке показано, как эта мощность распределяется между нагревом сота и нагревом воздуха в улочке:
Видно, что в улочку отдается примерно вдвое больше, чем соту - 11.36 мВт против 6.71 мВт. Составляющая dU/dt показывает мощность, идущую на нагрев пчелы пока режим неустановившийся; видно, что в установившемся режиме она стремится к нулю. Баланс мощности сходится с точностью ~10⁻¹⁵ мВт (машинная точность): P_met = Q_улочка + Q_соты + Q_испарение + dU/dt.
Воздух в ячейках под пчелой на соте заметно прогревается, это видно на тепловой карте кластера:
если в соседних с центральной ячейкой без пчелы на соте температура составляет около 12.45°С, то тоже в соседней с центральной ячейкой, но находящейся под тораксом пчелы на соте она достигает 12.71°С, а под головой - 12.55°С.
Вывод: модель пчелы на соте есть, можем двигаться дальше - строить модель улочки.
🐝 Python-скрипт
"""
=============================================================================
Программа к заметке "Пчела на соте: теплообмен на границе кластера"
=============================================================================
Модель теплообмена пчелы, сидящей НА ПОВЕРХНОСТИ сота (в составе кластера).
Развитие модели пчелы В ЯЧЕЙКЕ (fast0.py).
Идея: пчела на поверхности — граница между тёплым нутром кластера и холодной
средой. Её тело (3 соосных цилиндра: голова, торакс, брюшко) лежит над тремя
соседними ячейками и обменивается теплом с РАЗНОЙ средой сверху и снизу:
улочка (T_ambient)
═══════════════════════════════════════════
┌────┐ ┌──────────┐ ┌────────────┐
│Голова│ │ Торакс │ │ Брюшко │ ← пчела (вид сверху)
│1.5×4│ │ 4.5×4 │ │ 6×4 мм │
└──┬─┘ └────┬─────┘ └─────┬──────┘
┌─────▼───┐ ┌────▼─────┐ ┌─────▼──────┐
│(Y,1,−1) │ │ (Y,1,0) │ │ (Y,1,+1) │ ← ячейки сота
└─────────┘ └──────────┘ └────────────┘
Ключевые особенности модели:
* Дорсально-вентральная асимметрия теплообмена:
- ДОРСАЛЬНАЯ (спинная, верхняя) сторона → улочка (T_ambient);
- ВЕНТРАЛЬНАЯ (брюшная, нижняя) сторона → ячейки сота (T_air, T_wall).
* T_sense (температура, измеряемая усиками) = T воздуха в ячейке под ГОЛОВОЙ,
т.е. в (Y,1,−1). Усики направлены вперёд-вниз, поэтому «термометр» пчелы
отделён от источника тепла (торакса) — разнесённая обратная связь.
* Доли перекрытия входа ячейки выводятся из геометрии проекции цилиндра:
cover = (d_bee · h_segment) / A_entrance
голова ≈ 25%, торакс ≈ 75%, брюшко ≈ 100%.
* Вентральный теплообмен идёт ЧЕРЕЗ ПЕРЕКРЫТУЮ ЧАСТЬ ВХОДА (площадь
cover·A_entrance), а НЕ через кольцевой зазор вокруг тела.
Баланс мощности поверхностной пчелы (первый закон термодинамики):
P_met = Q_улочка + Q_соты + Q_испарение + dU/dt
где dU/dt — скорость накопления тепла телом (в стационаре → 0).
Архитектура: содержимое ячейки (States) и перекрытие входа (Cover) НЕЗАВИСИМЫ.
Ячейка под пчелой может содержать что угодно (пустая/расплод/мёд/перга/пчела) —
поверхностная пчела лишь добавляет перекрытие входа поверх содержимого.
=============================================================================
"""
import time
import math
import numpy as np
from numba import njit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# ==========================================================
# 1. ФИЗИЧЕСКИЕ ПАРАМЕТРЫ И ГЕОМЕТРИЯ
# ==========================================================
# --- Сетка и условия ---
# cols, rows = 11, 11 # Размер сетки (для расчёта кластера: 76, 57)
cols, rows = 7, 7
variant = 2 # 1: Тонкий сот, 2: Коричневый сот, 3: Толстый сот
T_ambient = 12.0 # Температура окружающего воздуха / улочки (°C)
# --- Пчела: три соосных цилиндра ---
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 # Коэффициент испарения (Вт/°C)
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 = 3500 # Удельная теплоёмкость тела (Дж/кг·К)
C_hd, C_th, C_ab = m_head * c, m_thorax * c, m_ab * c # Теплоёмкости (Дж/К)
# --- Воздух ---
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 # Перевод °C → К
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)
# Режим "печки" (0.0 — обычный; 45e-3 — нагрев центральной пчелы)
P_warm_up = 0.0
# ==========================================================
# 1b. ПЧЕЛА НА ПОВЕРХНОСТИ: АСИММЕТРИЯ ТЕПЛООБМЕНА
# ==========================================================
# Тело пчелы на поверхности мысленно делится по оси на две половины:
# ДОРСАЛЬНАЯ (верхняя) — открыта в улочку;
# ВЕНТРАЛЬНАЯ (нижняя) — обращена ко входу ячейки.
# Сверху и снизу среда разная, поэтому и теплообмен разный.
# --- ДОРСАЛЬНАЯ сторона (вверх, в улочку) ---
# Свободная конвекция над горизонтальной поверхностью (Вт/м²·К)
h_amb_surface = 8.0
# Дорсальные площади = верхняя половина боковой поверхности + торец(ы)
A_dorsal_th = 0.5 * A_th + s_r
A_dorsal_hd = 0.5 * (A_hd - s_r) + s_r
A_dorsal_ab = 0.5 * (A_ab - s_r) + s_r
# Сопротивления дорсальной конвекции (К/Вт)
R_th_amb_surf = 1.0 / (h_amb_surface * A_dorsal_th)
R_hd_amb_surf = 1.0 / (h_amb_surface * A_dorsal_hd)
R_ab_amb_surf = 1.0 / (h_amb_surface * A_dorsal_ab)
# --- Доли перекрытия входа ячейки (из геометрии проекции цилиндра) ---
# Проекция сегмента на вход = d_bee · h_segment; делим на площадь входа.
cover_head = min((d_bee * h_head) / A_entrance, 1.0) # ≈ 0.25
cover_thorax = min((d_bee * h_thorax) / A_entrance, 1.0) # ≈ 0.75
cover_ab = min((d_bee * h_ab) / A_entrance, 1.0) # ≈ 1.00
print(f"Доли перекрытия входа: голова={cover_head*100:.1f}%, "
f"торакс={cover_thorax*100:.1f}%, брюшко={cover_ab*100:.1f}%")
# --- ВЕНТРАЛЬНАЯ сторона (вниз, в ячейку) ---
# Пчела сидит сверху как «крышка». Тепло идёт ЧЕРЕЗ ПЕРЕКРЫТУЮ ЧАСТЬ ВХОДА
# (площадь cover·A_entrance), а не через кольцевой зазор вокруг тела.
# Площадь вентрального теплообмена = проекция сегмента:
A_cov_hd = cover_head * A_entrance
A_cov_th = cover_thorax * A_entrance
A_cov_ab = cover_ab * A_entrance
# Вентральная конвекция: тот же коэффициент h_entrance, что и для открытой
# части входа (вход — два параллельных пути: открытый и закрытый пчелой).
G_cell_hd = h_entrance * A_cov_hd # проводимость (Вт/К)
G_cell_th = h_entrance * A_cov_th
G_cell_ab = h_entrance * A_cov_ab
# Вентральное излучение: тело → стенки ячейки, через перекрытую площадь
se_cell_hd = sigma * eps_wall * A_cov_hd
se_cell_th = sigma * eps_wall * A_cov_th
se_cell_ab = sigma * eps_wall * A_cov_ab
# ==========================================================
# 2. ГЕНЕРАЦИЯ ТОПОЛОГИИ И СЕТКИ
# ==========================================================
# Типы СОДЕРЖИМОГО ячейки (не зависят от перекрытия входом!)
CELL_EMPTY = 0
CELL_BEE = 5
# Индексы в массиве состояний [13 элементов]:
# [0]=T_hd, [1]=T_th, [2]=T_ab, [3]=T_air, [4..12]=T_wall (6 бок. + 3 дна)
I_HD, I_TH, I_AB, I_AI = 0, 1, 2, 3
# Массивы сопротивлений/теплоёмкостей для Numba
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)
# --- Адреса ячеек (две стороны сота: Y и ^) ---
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):
"""Адрес соседа для грани face_id (0-5: боковые, 6-8: дно)."""
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
# --- Содержимое ячеек (States). Одна пчела в ячейке (центральная). ---
States = np.full(N_CELLS, CELL_EMPTY, dtype=np.int32)
bee_cells = [(0, 'Y', 0, 0)]
# 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
# print(f">>> НАША: cols={cols}, rows={rows}, пчёл в ячейках={len(bee_cells)}")
# --- Пчела на поверхности: какие ячейки перекрывает ---
# Торакс над (Y,1,0), голова над (Y,1,−1), брюшко над (Y,1,+1)
surf_bee_thorax_cell = (0, 'Y', 0, 1)
surf_bee_head_cell = (0, 'Y', -1, 1)
surf_bee_ab_cell = (0, 'Y', 1, 1)
surf_bee_thorax_idx = cell_to_index[surf_bee_thorax_cell]
surf_bee_head_idx = cell_to_index[surf_bee_head_cell]
surf_bee_ab_idx = cell_to_index[surf_bee_ab_cell]
# --- Перекрытие входа (Cover) — НЕЗАВИСИМО от содержимого ячейки ---
# В будущем ячейка может содержать расплод/мёд/пергу/пчелу — Cover просто
# добавляет перекрытие входа поверх любого содержимого.
Cover = np.zeros(N_CELLS, dtype=np.float64)
Cover[surf_bee_head_idx] = cover_head
Cover[surf_bee_thorax_idx] = cover_thorax
Cover[surf_bee_ab_idx] = cover_ab
print(f"Ячейки под поверхностной пчелой (содержимое не меняется, добавлено перекрытие):")
print(f" Голова: {surf_bee_head_cell} cover={Cover[surf_bee_head_idx]:.2f}")
print(f" Торакс: {surf_bee_thorax_cell} cover={Cover[surf_bee_thorax_idx]:.2f}")
print(f" Брюшко: {surf_bee_ab_cell} cover={Cover[surf_bee_ab_idx]:.2f}")
# Cover = np.zeros(N_CELLS, dtype=np.float64)
# ==========================================================
# 3. РАСЧЁТ НАЧАЛЬНЫХ ТЕМПЕРАТУР
# ==========================================================
def metabolic_power(T):
"""Мощность термогенеза пчелы (Вт) в зависимости от T_sense."""
return max((29.133 - 0.739 * T) * 1e-3, Pmin)
def calculate_initial_temperatures(T_air_val):
"""Стационарные температуры пчелы В ЯЧЕЙКЕ при заданной T воздуха."""
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 = metabolic_power(T_air_val)
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)
T_bee_avg_init = (m_head * Thd_init + m_thorax * Tth_init + m_ab * Tab_init) / (m_head + m_thorax + m_ab)
# Начальные температуры поверхностной пчелы = те же, что у пчелы в ячейке:
# при T_sense = T_ambient = 12°C метаболизм одинаков. Разница в теплопотерях
# (дорсальная конвекция в улочку) проявится уже в динамике.
Thd_surf_init, Tth_surf_init, Tab_surf_init = Thd_init, Tth_init, Tab_init
print(f"Начальные T (обе пчелы, T_sense={T_ambient}°C): "
f"Hd={Thd_init:.2f}, Th={Tth_init:.2f}, Ab={Tab_init:.2f}")
# ==========================================================
# 4. ПОДГОТОВКА МАТРИЦ ДЛЯ 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
# Состояние поверхностной пчелы: [T_hd, T_th, T_ab]
T_surf_bee = np.array([Thd_surf_init, Tth_surf_init, Tab_surf_init], dtype=np.float64)
# ==========================================================
# 5. ЯДРО СИМУЛЯЦИИ (NUMBA JIT)
# ==========================================================
@njit(fastmath=True, cache=True)
def do_cluster_step_surf(T_cur, T_next, Topology, Neighbor_Faces, States, Cover, C_grid,
dt, T_ambient_val, heater_idx,
surf_thorax_idx, surf_head_idx, surf_ab_idx,
T_surf_hd, T_surf_th, T_surf_ab):
"""
Один шаг тепловой сети кластера.
Вход ячейки может быть перекрыт поверхностной пчелой (Cover[i] > 0)
НЕЗАВИСИМО от содержимого ячейки (States[i]).
Модель входа для перекрытой ячейки (содержимое — не пчела):
открытая часть (1−cover)·A_entrance : конвекция воздух ↔ T_ambient (h_entrance)
закрытая часть cover·A_entrance : конвекция воздух ↔ T_тела (h_entrance)
Возвращает: (Q_entrance_Y, Q_entrance_X) — теплопотери через входы (Вт).
"""
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]
cover = Cover[i]
Thd = T_cur[i, I_HD]
Tth = T_cur[i, I_TH]
Tab = T_cur[i, I_AB]
# --- Температура воздуха в ячейке (Tai) ---
if state == CELL_BEE:
# Пчела внутри ячейки: воздух греется от частей тела,
# охлаждается через торец брюшка (R_ab_)
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_
elif cover > 0.0:
# Вход частично перекрыт поверхностной пчелой
if i == surf_thorax_idx:
T_part = T_surf_th
elif i == surf_head_idx:
T_part = T_surf_hd
else:
T_part = T_surf_ab
open_area = (1.0 - cover) * A_entrance
covered_area = cover * A_entrance
# Два параллельных пути через вход
num = T_ambient_val * h_entrance * open_area + T_part * h_entrance * covered_area
den = h_entrance * A_entrance
else:
# Пустая ячейка, вход полностью открыт
num = T_ambient_val * hA
den = hA
# Вклад стенок (9 граней)
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)
elif cover > 0.0:
# В улочку теряет только ОТКРЫТАЯ часть входа
# (закрытая обменивается с пчелой — учтено в Q_cells поверхностной пчелы)
Q_cell = h_entrance * (1.0 - cover) * A_entrance * (Tai - T_ambient_val)
if is_side_Y:
Q_entrance_Y += Q_cell
else:
Q_entrance_X += Q_cell
else:
Q_cell = hA * (Tai - T_ambient_val)
if is_side_Y:
Q_entrance_Y += Q_cell
else:
Q_entrance_X += Q_cell
# --- Температуры стенок (9 граней) ---
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
# +++ Излучение ПОВЕРХНОСТНОЙ пчелы к стенкам покрытой ячейки
if cover > 0.0 and state != CELL_BEE:
if i == surf_thorax_idx:
T_part_rad = T_surf_th
elif i == surf_head_idx:
T_part_rad = T_surf_hd
else:
T_part_rad = T_surf_ab
p_surf_to_wall = sigma * eps_wall * (cover * A_entrance) * \
((T_part_rad - t0) ** 4 - (T_wall_avg - t0) ** 4)
else:
p_surf_to_wall = 0.0
# Что греет стенки: пчела внутри (p_s_total) ИЛИ поверхностная пчела (p_surf_to_wall)
p_rad_to_walls = p_s_total if state == CELL_BEE else p_surf_to_wall
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_to_wall_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_to_wall_k)
# --- Температуры пчелы В ЯЧЕЙКЕ ---
if state == CELL_BEE:
# Метаболизм зависит от T воздуха в ячейке (= T_sense для ячеечной пчелы)
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
@njit(fastmath=True, cache=True)
def do_surface_bee_step(T_surf, T_cur, dt, T_ambient_val,
surf_head_idx, surf_thorax_idx, surf_ab_idx):
"""
Расчёт температур пчелы НА ПОВЕРХНОСТИ + все тепловые потоки НАПРЯМУЮ.
Дорсально (в улочку): свободная конвекция (h_amb) + излучение → T_ambient,
площадь = дорсальная поверхность тела (A_dorsal).
Вентрально (в соты): конвекция (h_entrance) + излучение → T_air/T_wall,
площадь = перекрытая часть входа (A_cov = cover·A_entrance).
Возвращает: (P_met, Q_street, Q_cells, Q_evap, dU_dt)
P_met — мощность термогенеза (Вт)
Q_street — поток в улочку: дорсальная конвекция + излучение (Вт)
Q_cells — поток в соты: вентральная конвекция + излучение (Вт)
Q_evap — поток на испарение (Вт)
dU_dt — скорость накопления тепла телом (Вт)
Проверка (первый закон термодинамики):
P_met − Q_street − Q_cells − Q_evap − dU_dt == 0 (с машинной точностью)
В стационаре dU_dt → 0, тогда P_met = Q_street + Q_cells + Q_evap.
"""
T_hd = T_surf[0]
T_th = T_surf[1]
T_ab = T_surf[2]
# T_sense — температура, измеряемая усиками (в области головы)
T_sense = T_cur[surf_head_idx, I_AI]
# Средние температуры стенок ячеек под частями тела
T_wall_hd = 0.0
for f in range(9):
T_wall_hd += T_cur[surf_head_idx, 4 + f]
T_wall_hd /= 9.0
T_wall_th = 0.0
for f in range(9):
T_wall_th += T_cur[surf_thorax_idx, 4 + f]
T_wall_th /= 9.0
T_wall_ab = 0.0
for f in range(9):
T_wall_ab += T_cur[surf_ab_idx, 4 + f]
T_wall_ab /= 9.0
# Температуры воздуха в ячейках под частями тела
T_air_hd = T_cur[surf_head_idx, I_AI]
T_air_th = T_cur[surf_thorax_idx, I_AI]
T_air_ab = T_cur[surf_ab_idx, I_AI]
# Метаболизм (зависит от T_sense)
P_met = (29.133 - 0.739 * T_sense) * 1e-3
if P_met < Pmin:
P_met = Pmin
# Динамические сопротивления теплопроводности перешейков
l_dyn = (1.3493 * np.exp(-0.058 * T_sense)) * 1e-3
R_th_hd_dyn = l_dyn / (k_hd * s_r)
R_th_ab_dyn = l_dyn / (k_ab * s_r)
# --- Излучение дорсальное (→ T_ambient, площадь тела A_dorsal) ---
p_rad_hd_d = -sigma * eps * ((T_hd - t0)**4 - (T_ambient_val - t0)**4) * A_dorsal_hd
p_rad_th_d = -sigma * eps * ((T_th - t0)**4 - (T_ambient_val - t0)**4) * A_dorsal_th
p_rad_ab_d = -sigma * eps * ((T_ab - t0)**4 - (T_ambient_val - t0)**4) * A_dorsal_ab
# --- Излучение вентральное (→ T_wall, площадь перекрытого входа A_cov) ---
p_rad_hd_v = -se_cell_hd * ((T_hd - t0)**4 - (T_wall_hd - t0)**4)
p_rad_th_v = -se_cell_th * ((T_th - t0)**4 - (T_wall_th - t0)**4)
p_rad_ab_v = -se_cell_ab * ((T_ab - t0)**4 - (T_wall_ab - t0)**4)
# Испарение с одной поверхности (головы или брюшка)
q_evap = ke0 * T_th
# ==========================================================
# ПРЯМЫЕ ТЕПЛОВЫЕ ПОТОКИ (положительные = от пчелы наружу)
# ==========================================================
# В УЛОЧКУ (дорсально): конвекция + излучение
Q_street_conv = ((T_hd - T_ambient_val) / R_hd_amb_surf +
(T_th - T_ambient_val) / R_th_amb_surf +
(T_ab - T_ambient_val) / R_ab_amb_surf)
Q_street_rad = -(p_rad_hd_d + p_rad_th_d + p_rad_ab_d)
Q_street = Q_street_conv + Q_street_rad
# В СОТЫ (вентрально): конвекция через перекрытый вход + излучение
Q_cells_conv = (G_cell_hd * (T_hd - T_air_hd) +
G_cell_th * (T_th - T_air_th) +
G_cell_ab * (T_ab - T_air_ab))
Q_cells_rad = -(p_rad_hd_v + p_rad_th_v + p_rad_ab_v)
Q_cells = Q_cells_conv + Q_cells_rad
# ИСПАРЕНИЕ (с головы и брюшка)
Q_evap = 2.0 * q_evap
# ==========================================================
# УРАВНЕНИЯ БАЛАНСА (для интегрирования температур)
# ==========================================================
# Голова: теплопроводность от торакса = дорс.конв + дорс.изл + вентр.конв + вентр.изл + испарение
dT_hd = ((T_th - T_hd) / R_th_hd_dyn
- (T_hd - T_ambient_val) / R_hd_amb_surf
- G_cell_hd * (T_hd - T_air_hd)
+ p_rad_hd_d + p_rad_hd_v - q_evap)
# Торакс: метаболизм = теплопроводность + дорс.конв + дорс.изл + вентр.конв + вентр.изл
dT_th = (P_met
- (T_th - T_hd) / R_th_hd_dyn - (T_th - T_ab) / R_th_ab_dyn
- (T_th - T_ambient_val) / R_th_amb_surf
- G_cell_th * (T_th - T_air_th)
+ p_rad_th_d + p_rad_th_v)
# Брюшко: теплопроводность от торакса = дорс.конв + дорс.изл + вентр.конв + вентр.изл + испарение
dT_ab = ((T_th - T_ab) / R_th_ab_dyn
- (T_ab - T_ambient_val) / R_ab_amb_surf
- G_cell_ab * (T_ab - T_air_ab)
+ p_rad_ab_d + p_rad_ab_v - q_evap)
# dU/dt — мощность, идущая на изменение внутренней энергии тела
# (dT_hd, dT_th, dT_ab имеют размерность Вт; C·dT/dt = dT_здесь)
dU_dt = dT_hd + dT_th + dT_ab
# Интегрирование (явный Эйлер)
T_surf[0] = T_hd + (dt / C_hd) * dT_hd
T_surf[1] = T_th + (dt / C_th) * dT_th
T_surf[2] = T_ab + (dt / C_ab) * dT_ab
return P_met, Q_street, Q_cells, Q_evap, dU_dt
# ==========================================================
# 6. ЦИКЛ СИМУЛЯЦИИ
# ==========================================================
dt_sim = 0.1 # Шаг интегрирования (с)
time_total = 600.0 # Длительность симуляции (с)
steps = int(time_total / dt_sim)
# --- Предаллокация логов ---
log_time = []
# Пчела в ячейке:
log_T_hd, log_T_th, log_T_ab = [], [], []
log_T_ai_bee, log_T_bee_avg = [], []
log_Power = []
# Пчела на поверхности:
log_T_hd_surf, log_T_th_surf, log_T_ab_surf = [], [], []
log_T_sense_surf = []
log_Power_surf = []
# Потоки поверхностной пчелы (все — прямые расчёты):
log_Q_street_surf, log_Q_cells_surf, log_Q_evap_surf = [], [], []
log_dU_dt_surf = [] # накопление тепла (в стационаре → 0)
log_balance_residual = [] # невязка баланса (для контроля, ~машинная точность)
# Теплопотери кластера:
log_Q_Y, log_Q_X = [], []
print(f"\nСетка {cols}x{rows}. Запуск...")
t_idx = bee_idx[0]
# Прогрев Numba (первый вызов — компиляция)
do_cluster_step_surf(T_grid_A, T_grid_B, Topology, Neighbor_Faces, States, Cover, C_grid,
dt_sim, T_ambient, bee_idx[0],
surf_bee_thorax_idx, surf_bee_head_idx, surf_bee_ab_idx,
T_surf_bee[0], T_surf_bee[1], T_surf_bee[2])
do_surface_bee_step(T_surf_bee, T_grid_A, dt_sim, T_ambient,
surf_bee_head_idx, surf_bee_thorax_idx, surf_bee_ab_idx)
# --- Основной цикл ---
start_time = time.time()
T_cur, T_next = T_grid_A, T_grid_B
for step in range(steps):
# Шаг 1: тепловая сеть кластера
Q_Y, Q_X = do_cluster_step_surf(
T_cur, T_next, Topology, Neighbor_Faces, States, Cover, C_grid,
dt_sim, T_ambient, bee_idx[0],
surf_bee_thorax_idx, surf_bee_head_idx, surf_bee_ab_idx,
T_surf_bee[0], T_surf_bee[1], T_surf_bee[2])
# Swap буферов
T_cur, T_next = T_next, T_cur
# Шаг 2: пчела на поверхности (все потоки напрямую)
P_surf, Q_street, Q_cells, Q_evap, dU_dt = do_surface_bee_step(
T_surf_bee, T_cur, dt_sim, T_ambient,
surf_bee_head_idx, surf_bee_thorax_idx, surf_bee_ab_idx)
# Шаг 3: дорсальный поток поверхностной пчелы — в теплопотери стороны Y
Q_Y += Q_street
# --- Логирование (каждый 10-й шаг) ---
if step % 10 == 0:
t = step * dt_sim
log_time.append(t)
# Пчела в ячейке
log_T_hd.append(T_cur[t_idx, I_HD])
log_T_th.append(T_cur[t_idx, I_TH])
log_T_ab.append(T_cur[t_idx, I_AB])
log_T_ai_bee.append(T_cur[t_idx, I_AI])
T_bee_avg = (m_head * T_cur[t_idx, I_HD] + m_thorax * T_cur[t_idx, I_TH] +
m_ab * T_cur[t_idx, I_AB]) / (m_head + m_thorax + m_ab)
log_T_bee_avg.append(T_bee_avg)
P_rest = max((29.133 - 0.739 * T_cur[t_idx, I_AI]) * 1e-3, Pmin)
log_Power.append((P_rest + P_warm_up) * 1000)
# Пчела на поверхности
log_T_hd_surf.append(T_surf_bee[0])
log_T_th_surf.append(T_surf_bee[1])
log_T_ab_surf.append(T_surf_bee[2])
log_T_sense_surf.append(T_cur[surf_bee_head_idx, I_AI])
log_Power_surf.append(P_surf * 1000)
# Потоки (мВт)
log_Q_street_surf.append(Q_street * 1000)
log_Q_cells_surf.append(Q_cells * 1000)
log_Q_evap_surf.append(Q_evap * 1000)
log_dU_dt_surf.append(dU_dt * 1000)
# Невязка баланса (контроль корректности кода)
residual = (P_surf - Q_street - Q_cells - Q_evap - dU_dt) * 1000
log_balance_residual.append(residual)
# Теплопотери кластера
log_Q_Y.append(Q_Y * 1000)
log_Q_X.append(Q_X * 1000)
time_calc = time.time() - start_time
print(f"Общее время: {time_calc:.3f} c ({steps / time_calc:.0f} шагов/сек)")
# ==========================================================
# 7. ВИЗУАЛИЗАЦИЯ
# ==========================================================
# --- График 1: Динамика температур ---
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 8), sharex=True)
# Пчела в ячейке
ax1.scatter(log_time, log_T_th, s=2, c='r', label='Торакс', zorder=5)
ax1.plot(log_time, log_T_hd, 'b-', label='Голова', linewidth=1.2)
ax1.plot(log_time, log_T_ab, 'g-', label='Брюшко', linewidth=1.2)
ax1.plot(log_time, log_T_ai_bee, 'c--', label='Воздух в ячейке', linewidth=1.2)
ax1.plot(log_time, log_T_bee_avg, 'k:', label='Средняя T', linewidth=1.5)
ax1.axhline(T_ambient, color='gray', linestyle=':', label=f'T_amb={T_ambient}°C')
ax1.set_ylabel('Температура (°C)')
ax1.set_title(f'Пчела В ЯЧЕЙКЕ (центральная). Вариант {variant}')
ax1.legend(fontsize=8); ax1.grid(True, linestyle=':')
# Пчела на поверхности
ax2.plot(log_time, log_T_th_surf, 'r-', label='Торакс', linewidth=1.5)
ax2.plot(log_time, log_T_hd_surf, 'b-', label='Голова', linewidth=1.2)
ax2.plot(log_time, log_T_ab_surf, 'g-', label='Брюшко', linewidth=1.2)
ax2.plot(log_time, log_T_sense_surf, 'c--', label='T_sense (воздух под головой)', linewidth=1.2)
ax2.axhline(T_ambient, color='gray', linestyle=':', label=f'T_amb={T_ambient}°C')
ax2.set_xlabel('Время (с)'); ax2.set_ylabel('Температура (°C)')
ax2.set_title('Пчела НА ПОВЕРХНОСТИ СОТА')
ax2.legend(fontsize=8); ax2.grid(True, linestyle=':')
plt.tight_layout(); plt.show()
# --- График 2: Мощность термогенеза (в ячейке vs на поверхности) ---
fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(log_time, log_Power, 'orange', linewidth=2.2, label='Пчела В ЯЧЕЙКЕ')
ax.plot(log_time, log_Power_surf, 'darkred', linewidth=2.2, linestyle='--',
label='Пчела НА ПОВЕРХНОСТИ')
ax.set_xlabel('Время (с)')
ax.set_ylabel('Мощность термогенеза P_met (мВт)')
ax.set_title('Мощность термогенеза:\n пчела в ячейке vs пчела на поверхности')
ax.legend(); ax.grid(True, linestyle=':')
plt.tight_layout(); plt.show()
# --- График 3: Баланс мощности поверхностной пчелы ---
# Все слагаемые первого закона термодинамики:
# P_met = Q_улочка + Q_соты + Q_испарение + dU/dt
fig, ax = plt.subplots(figsize=(6, 6))
ax.plot(log_time, log_Power_surf, 'k-', linewidth=2.5, label='P_met (генерация)')
ax.plot(log_time, log_Q_street_surf, 'r-', linewidth=1.8, label='Q → улочка (дорсально)')
ax.plot(log_time, log_Q_cells_surf, 'b-', linewidth=1.8, label='Q → соты (вентрально)')
ax.plot(log_time, log_Q_evap_surf, 'g--', linewidth=1.5, label='Q испарение')
ax.plot(log_time, log_dU_dt_surf, color='purple', linewidth=1.5, label='dU/dt (накопление)')
# Сумма потоков — должна совпадать с P_met (проверка баланса)
sum_flows = (np.array(log_Q_street_surf) + np.array(log_Q_cells_surf) +
np.array(log_Q_evap_surf) + np.array(log_dU_dt_surf))
ax.plot(log_time, sum_flows, color='gray', linestyle=':', linewidth=2.0,
label='Σ потоков (= P_met)')
ax.axhline(0, color='k', linestyle='-', linewidth=0.5)
ax.set_xlabel('Время (с)')
ax.set_ylabel('Мощность (мВт)')
ax.set_title('Баланс мощности поверхностной пчелы:\n'
'P_met = Q_улочка + Q_соты + Q_испарение + dU/dt')
ax.legend(fontsize=8); ax.grid(True, linestyle=':')
plt.tight_layout(); plt.show()
# --- График 4: 2D-карта температур сот + пчела (тепловизор) ---
temp_Y_clean, temp_X_clean = {}, {}
for i, cell in enumerate(cells):
_, side, col, row = cell
T_air_final = T_cur[i, I_AI]
if side == 'Y':
temp_Y_clean[(col, row)] = T_air_final
else:
temp_X_clean[(col, row)] = T_air_final
size = 1.0
width = np.sqrt(3) * size
dy_hex = 1.5 * size
def hex_corner(center, i):
angle_rad = np.pi / 180 * (60 * i - 30)
return (center[0] + size * np.cos(angle_rad), center[1] + size * np.sin(angle_rad))
def primal_to_pixel(col, row):
return (width * (col + 0.5 * (row & 1)), -dy_hex * row)
def dual_to_pixel(col, row):
return (width * (col + 0.5 * (1 - (row & 1))), -dy_hex * row - size / 2)
fig, ax = plt.subplots(figsize=(9, 7))
ax.set_aspect('equal')
# Диапазон шкалы: включаем температуры пчелы
all_temps = list(temp_Y_clean.values()) + list(temp_X_clean.values())
all_temps += [T_surf_bee[0], T_surf_bee[1], T_surf_bee[2]]
vmin, vmax = min(all_temps), max(all_temps)
# Гексагоны стороны Y
for col in range(-int(cols / 2), int(cols / 2) + 1):
for row in range(-int(rows / 2), int(rows / 2) + 1):
if (col, row) in temp_Y_clean:
x, y = primal_to_pixel(col, row)
corners = [hex_corner((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
color = plt.cm.viridis((temp_Y_clean[(col, row)] - vmin) / max(vmax - vmin, 1e-5))
ax.fill(hex_x, hex_y, color=color, edgecolor='steelblue', linewidth=1.2, alpha=0.8)
if (col, row) == (0, 0):
ax.fill(hex_x, hex_y, facecolor='none', edgecolor='black', linewidth=3)
ax.text(x, y, f"{temp_Y_clean[(col, row)]:.2f}", ha='center', va='center',
fontsize=12, color='b')
# Гексагоны стороны ^
for col in range(-int(cols / 2), int(cols / 2) + 1):
for row in range(-int(rows / 2), int(rows / 2) + 1):
if (col, row) in temp_X_clean:
x, y = dual_to_pixel(col, row)
corners = [hex_corner((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
ax.plot(hex_x + (hex_x[0],), hex_y + (hex_y[0],), color='crimson', linewidth=1.2)
ax.text(x, y, f"{temp_X_clean[(col, row)]:.2f}", ha='center', va='center',
fontsize=12, color='crimson')
# Пчела на поверхности: три прямоугольника (вид сверху = проекции сегментов)
scale = size / a # масштаб: map units на метр
w_bee_map = d_bee * scale # ширина тела (4 мм)
l_head_map = h_head * scale # длина головы (1.5 мм)
l_thorax_map = h_thorax * scale # длина торакса (4.5 мм)
l_ab_map = h_ab * scale # длина брюшка (6 мм)
# Зазор между сегментами (длина перешейка, зависит от T_sense)
T_sense_final = T_cur[surf_bee_head_idx, I_AI]
l_gap_map = (1.3493 * np.exp(-0.058 * T_sense_final)) * 1e-3 * scale
# Привязка: центр торакса = центр ячейки (Y,1,0)
cx_th, cy_th = primal_to_pixel(surf_bee_thorax_cell[2], surf_bee_thorax_cell[3])
cx_head = cx_th - l_thorax_map / 2 - l_gap_map - l_head_map / 2
cx_ab = cx_th + l_thorax_map / 2 + l_gap_map + l_ab_map / 2
cy_bee = cy_th
# Прямоугольники (полупрозрачные, без подписей T — видны ячейки под пчелой)
for cx_seg, l_seg, T_part in [(cx_head, l_head_map, T_surf_bee[0]),
(cx_th, l_thorax_map, T_surf_bee[1]),
(cx_ab, l_ab_map, T_surf_bee[2])]:
norm_T = (T_part - vmin) / max(vmax - vmin, 1e-5)
rect = mpatches.Rectangle((cx_seg - l_seg / 2, cy_bee - w_bee_map / 2), l_seg, w_bee_map,
facecolor=plt.cm.viridis(norm_T), edgecolor='black',
linewidth=1.5, alpha=0.45, zorder=10)
ax.add_patch(rect)
# Осевая линия тела
ax.plot([cx_head - l_head_map/2, cx_ab + l_ab_map/2], [cy_bee, cy_bee],
'k-', linewidth=0.8, alpha=0.4, zorder=9)
# ax.text(cx_th, cy_bee + w_bee_map / 2 + 0.2, 'Пчела на поверхности',
# ha='center', va='bottom', fontsize=9, fontweight='bold',
# bbox=dict(boxstyle='round,pad=0.2', facecolor='yellow', alpha=0.7))
ax.axis('off'); ax.autoscale_view()
sm = plt.cm.ScalarMappable(cmap=plt.cm.viridis, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array([])
plt.colorbar(sm, ax=ax, shrink=0.8).set_label('Температура (°C)', fontsize=10)
plt.title(f'Тепловизор: температурное поле сот + пчела на поверхности.\n'
f'Вариант {variant}, t = {time_total:.0f} с. '
f'(Синяя — сторона "Y", Красная — сторона "^")')
plt.tight_layout(); plt.show()
# ==========================================================
# 8. ФИНАЛЬНЫЙ ВЫВОД
# ==========================================================
print(f"\n{'='*60}")
print(f" РЕЗУЛЬТАТЫ (t = {time_total:.0f} с)")
print(f"{'='*60}")
print(f"\n Пчела В ЯЧЕЙКЕ (центральная):")
print(f" Торакс: {log_T_th[-1]:.2f} °C")
print(f" Голова: {log_T_hd[-1]:.2f} °C")
print(f" Брюшко: {log_T_ab[-1]:.2f} °C")
print(f" Воздух: {log_T_ai_bee[-1]:.2f} °C")
print(f" Мощность: {log_Power[-1]:.2f} мВт")
print(f"\n Пчела НА ПОВЕРХНОСТИ:")
print(f" Торакс: {log_T_th_surf[-1]:.2f} °C")
print(f" Голова: {log_T_hd_surf[-1]:.2f} °C")
print(f" Брюшко: {log_T_ab_surf[-1]:.2f} °C")
print(f" T_sense: {log_T_sense_surf[-1]:.2f} °C")
print(f"\n БАЛАНС МОЩНОСТИ поверхностной пчелы (прямые расчёты):")
print(f" P_met: {log_Power_surf[-1]:.3f} мВт (генерация)")
print(f" Q → улочка: {log_Q_street_surf[-1]:.3f} мВт (дорсально)")
print(f" Q → соты: {log_Q_cells_surf[-1]:.3f} мВт (вентрально)")
print(f" Q испарение: {log_Q_evap_surf[-1]:.3f} мВт")
print(f" dU/dt: {log_dU_dt_surf[-1]:.3f} мВт (→ 0 в стационаре)")
print(f"\n ВЕРИФИКАЦИЯ:")
print(f" Макс. невязка баланса: {max(abs(r) for r in log_balance_residual):.2e} мВт")
print(f" (должна быть ~1e-12 — машинная точность)")
print(f" Q_cells > 0? {'ДА ✓' if log_Q_cells_surf[-1] > 0 else 'НЕТ ✗'}")
print(f" dU/dt → 0? {'ДА ✓' if abs(log_dU_dt_surf[-1]) < 0.01 else 'ещё нет (сот прогревается)'}")
print(f"\n ТЕПЛОПОТЕРИ КЛАСТЕРА:")
print(f" Сторона Y: {log_Q_Y[-1]:.2f} мВт (включая дорсальный поток пчелы)")
print(f" Сторона ^: {log_Q_X[-1]:.2f} мВт")
print(f"{'='*60}")
Комментарии
Отправить комментарий