Введение
Жизнь пчелиной семьи - сложный процесс, меняющий свое содержание от сезона к сезону: она категорически отличается, например, зимой от лета. Поэтому и модели теплового обмена будут разные для каждого времени года. Один из важнейших периодов и, пожалуй, сложнейший с точки зрения тепловой модели - это период выращивания расплода, который продолжается практически полгода. В этой заметке фиксирую текущее состояние модели: как устроен обогрев расплода, пока без учета улочки, т.е. на одном соте.
Расплод семья располагает в центре сота - так эффективней обеспечить его теплом и пищей. Именно целый сот, а не его небольшой фрагмент нужен для того, чтобы промоделировать и понять, как происходит обогрев расплода.
Повышение быстродействия модели позволяет перейти от моделирования небольших фрагментов сота к любому стандартному или нестандартному размеру реальных сотов и рамок. Для определенности возьмем наиболее распространенный размер - Дадана. Рамка Дадана 435×300 мм (за вычетом толщины реек 415×268 мм) моделируется гексагональной сеткой 76×57 ячеек (шаг 5.4 мм).
1. Распределение расплода
Расплод разместим эллиптическим пятном в центре рамки. Полуоси эллипса: 30% ширины и 35% высоты рамки. Внутри эллипса каждая ячейка с вероятностью 92% содержит расплод; остальные 8% — пустые (пропуски, как в реальном соте).
Возраст. Для данной конкретной стороны рамки положим средний возраст расплода 300 часов (~12.5 суток), разброс ±18 ч. По радиусу пятна задан градиент: к центру расплод старше, к краю - моложе (матка движется по спирали от центра). Случайный шум ±6 ч имитирует нерегулярность засева.
Запечатывание. При возрасте 203 ч (~8.5 суток) ячейка запечатывается восковой крышечкой. В модели крышечка:
- добавляет теплоёмкость
C_cap = 10.35 мДж/Кк ячейке с расплодом; - открывает канал прямого контакта «торакс пчелы → ячейка с расплодом».
Масса, теплоемкость и термогенез эмбриона описаны в заметке Статус ячейки с расплодом 2
Напомню.
-
До запечатывания масса растёт экспоненциально:
m(t)=m_egg⋅exp(t/τ), m_egg=0.000289 мг, τ=15.5 ч
К моменту запечатывания(t = 203 ч): m = 132 мг. После запечатывания масса линейно уменьшается (куколка расходует запасы):
m(t)=132−0.12⋅(t−203) мг
К выходу взрослой пчелы(t ≈ 502 ч): m ≈ 96 мг. -
Для ячейки с расплодом возраста 300 ч:
P_brood(300) = 10⁻⁵·300² − 0.005019·300 + 0.877 ≈ 0.27 мВт
На всю сторону рамки (в пятне получилось 1313 ячеек расплода): суммарная мощность ≈ 0.35 Вт. Это существенный внутренний источник (350 мВт - это примерно 7 пчел-грелок) - расплод сам себя греет, но, как увидим далее, этого недостаточно. -
Теплоёмкость. В модели ячейка с расплодом имеет три тепловых узла:
Расплод (узел I_TH):C = m_embryo · c_brood + C_cap = 0.42 + 0.01 = 0.43 Дж/К
Стенки (9 узлов):C_wall = 6·C_lat + 3·C_bot = 0.125 Дж/К
Воздух:C_air = 0.00035 Дж/К(пренебрежимо мала)
Расплод и стенки - разные тела: расплод греет стенки конвекцией и излучением, стенки передают тепло соседям теплопроводностью. Крышечка (10 мДж/К) добавляется к расплоду, поскольку прижата к нему и имеет с ним хороший тепловой контакт.2. Пчёлы-грелки в ячейках
Пчёлы-грелки размещаются в пустых ячейках среди расплода - тех самых 8% пропусков. Для моделирования из пропусков отбирается 30% (минимум 2), равномерно по азимуту вокруг пятна расплода. В нашем случае: 33 грелки на 1313 ячеек расплода.
Пчела в ячейке моделируется, как и ранее, тремя соосными цилиндрами: голова(Ø4×1.5 мм), торакс(Ø4×4.5 мм), брюшко(Ø4×6 мм). Зазор между пчелой и стенками ячейки (~0.65 мм на сторону) заполнен воздухом; конвекция в зазоре описывается коэффициентом
h_in = k_air / ((d_cell − d_bee)/2) ≈ 40 Вт/(м²·К). Брюшко "смотрит" из ячейки и обменивается теплом с воздухом улочки.
Метаболизм. Базовая мощность определяется температурой, которую ощущают усики:
P_met = max(29.133−0.739⋅T_sense, 0.2), мВт
Для грелки в ячейкеT_sense = T_airэтой ячейки (усики направлены вперёд, вглубь ячейки). ПриT_air = 34°C: P_met ≈ 4 мВт. Этого достаточно, чтобы поддерживать торакс около ~35°C - на 1° выше среды. Но для обогрева расплода этого мало.
Режим печки. Пчела-грелка способна кратковременно увеличить тепловыделение, сокращая летательные мышцы торакса без движения крыльев. В модели это добавкаP_heater = 45 мВтк метаболизму торакса. Итого: ~49 мВт. По данным литературы (Stabentheiner, 1987; UMT thesis), торакс грелки при этом разогревается до 43–44°C, максимально до 47°C.
Термостат: один датчик - усики. Ключевое биологическое уточнение: пчела регулирует нагрев не по температуре собственного тела и не по температуре расплода, а по температуре воздуха, которую ощущают усики, т.е. обратная связь - косвенная.
Предположим, что логика термостатирования у пчелы такая:
измеренная температура ниже 35°С - включить режим печки;
измеренная температура выше 36°С - выключить режим печки;
температура головы выше 47°С - аварийное выключение режима печки.
Гистерезис 1°C предотвращает дребезг, аварийный потолок 47°C- денатурацию белков.
Длительность нагрева. Пчёлы-грелки остаются в пустых ячейках неподвижно до 30 минут [20, с.176 и др.]. В модели максимальную продолжительность включения режима печки зададим 1800 с (30 мин).Результат (без обсидчиков,
T_ambient = 34°C):- Воздух в ячейке грелки: 35–36°C (термостат держит в этом диапазоне)
- Брюшко грелки: 36–37.5°C (на рисунке — цвет кружка)
- Торакс грелки: 37°C (печка ВЫКЛ) → 37–43°C (печка ВКЛ), на рисунке не показан, т.к. экранируется брюшком
- Расплод в центре пятна: ~34.7°C (без обсидчиков до цели 35.5 не дотягивает)
- Грелки работают циклически: печка ВКЛ ~28% времени (в среднем по 33 грелкам). Период цикла определяется тепловой инерцией ячейки: нагрев воздуха до 36°C (выкл) → остывание до 35°C (вкл).
На рисунке пчелы показаны кружками Ø4 мм (диаметр пчелы), цвет - температура брюшка; обводка: красная - печка ВКЛ, оранжевая - ВЫКЛ.3. "Обсидчики" - так для краткости будем называть пчёл на поверхности сота
сидят сверху на соте, частично перекрывая входы ячеек своим телом. В модели 351 обсидчик, размещены эллиптическим пятном (чуть шире пятна расплода: полуоси 32%×37%).
Ориентация и упаковка. Каждая пчела задаётся парой
(thorax_cell, orientation), гдеorientation ∈ {0..5} —направление на голову (к одному из 6 боковых соседей, см. Пчела на соте). Брюшко смотрит в противоположного соседа. Три ячейки на пчелу (голова, торакс, брюшко) не должны пересекаться с ячейками других пчёл - генератор выборки перебирает кандидатов и проверяет уникальность.Дорсально-вентральная асимметрия. Среда сверху и снизу разная, поэтому теплообмен разный:
Направление Среда Механизм Коэффициент Дорсально (спина, вверх) Воздух улочки Свободная конвекция + излучение h_amb = 8 Вт/(м²·К) Вентрально (вниз) Ячейки сота Через перекрытую часть входа h_entrance = 11 Вт/(м²·К) Доли перекрытия входа. Проекция части тела на площадь входа (
A_entrance ≈ 24 мм²):Часть тела Экран Физический смысл Голова 25% 75% входа открыто Торакс 75% 25% входа открыто Брюшко 100% вход полностью закрыт Закрытая часть входа «видит» температуру тела пчелы (конвекция + излучение к стенкам ячейки). Открытая часть — воздух улочки.
Контакт через крышечку. Если торакс обсидчика находится над запечатанной ячейкой расплода, тепло идёт не через воздух, а напрямую через восковую крышечку.
Контакт в ~14 раз эффективнее конвекции. В режиме печки пчела прижимает торакс к крышечке (boost=1.5), что ещё улучшает теплоотдачу;g_contact = 150 Вт/(м²·К) —оценочное значение для прижима хитинового торакса к восковой мембране толщиной 0.15 мм. Реальное значение зависит от силы прижима и площади фактического контакта; требует экспериментальной проверки.Термостат обсидчика. Датчик - усики, направленные вперёд-вниз.
T_sense- температура воздуха в ячейке под головой (не под тораксом!) или температура крышечки. Это «разнесённая обратная связь»: источник тепла (торакс) и датчик (усики) разделены на несколько миллиметров. Пчела может иметь горячий торакс (43°C) и отдавать тепло в ячейку под тораксом, но измеряемая температура в соседней ячейке (крышечки соседней ячейки) будет, например, 35–36°C - и пчела будет продолжать греть.Логика та же, что для грелок в ячейках:
Условие Действие T_sense < 35°C ВКЛЮЧИТЬ печку T_sense ≥ 36°C ВЫКЛЮЧИТЬ T_hd ≥ 47°C АВАРИЙНОЕ ВЫКЛ Таймер ≥ 1800 с ВЫКЛЮЧИТЬ (пора отдохнуть) Баланс мощности обсидчика:
P_met + P_heater = Q_street + Q_cells + Q_evap + dU/dtгде
Q_street —дорсальные потери (конвекция + излучение в улочку),Q_cells —вентральный поток (конвекция + контакт + излучение в ячейки),Q_evap — испарение,dU/dt- мощность нагрева самой пчелы (до установившегося режима). Все потоки считаются напрямую из физики; баланс - проверка. Невязка: ~10⁻¹⁵ мВт (машинная точность).Численный пример (обсидчик в режиме печки):
Поток Знач. Куда P_met + P_heater ~49 мВт источник (торакс) Контакт торакс → крышечка ~30 мВт в расплод (основной канал) Конвекция + излучение дорсально ~3.5 мВт в улочку Конвекция головы + брюшка вентрально ~2 мВт в воздух ячеек Излучение вентрально (к стенкам) ~1 мВт в стенки Испарение ~1.5 мВт потеря Примерно две трети мощности идёт в расплод. Поэтому обсидчики на запечатанном расплоде - главный инструмент обогрева, т.к. их гораздо больше, чем пчел-грелок в пустых ячейках и, как увидим дальше, они чаще включают режим печки.
Видим, что расплод прогревается неравномерно, некоторые ячейки перегреты, но, во-первых, это не постоянно и не смертельно, а, во-вторых, мы пока не учитываем рост температуры в улочке и поэтому не видим эффект усреднения температуры в улье (увидим это в следующей заметке).
На следующем рисунке показано, как могла бы выглядеть картина в тепловизоре:Визуально заметно, что плотность "включенных" обсидчиков больше на периферии расплода.
4. Итоговая температура расплода
Для нормального развития расплода требуется 33–36°C [32]. Ниже 33°C - гибель, выше 36°C - дефекты/отклонения развития. Модель должна удерживать расплод в этом коридоре при
T_ambient = 34°C(воздух в гнезде рядом с расплодной зоной).Иерархия нагрева. Три сценария:
Сценарий T расплода (средн/мин/макс), °C ΔT от T_amb, °C Ячеек с Т>34.5°С Без пчёл (1313) 34.40/34.22/34.51 +0.40 56 Только грелки (1313+33) 34.72/34.29/35.11 +0.72 1160 Грелки + обсидчики (1313+33+351) 35.29/35.01/38.94 +1.29°C 1313 Расплод без пчёл греет себя недостаточно: 0.27 мВт на ячейку при тепловом сопротивлении дают в среднем ΔT ≈ 0.40°C.
Грелки в ячейках добавляют еще 0.32°С (+0.72).Обсидчики также включают режим печки по алгоритму локального термостатирования и доводят общий нагрев до 1.29°С, расплод нагревается в среднем до 35.29°С.
Как видим, основная роль в поддержании температуры расплода принадлежит обсидке.
В таблице даны установившиеся значения после моделирования часа (3600 секунд) реального времени. Однако представляет интерес и динамика термостатирования. На следующем рисунке представлено изменение во времени средней температуры расплода, пчел и как меняется доля пчел в ячейках и на соте, включивших режиме печки.
Видим, что в установившемся режиме 33 пчелы в ячейках (грелки) способны поддерживать среднюю температуру всего расплода, но для этого им приходится держать среднюю температуру торакса около 38°С, при этом в режиме печки должны находиться около четверти всех пчел грелок. При наличии обсидки (351 пчела на соте) грелки в ячейках вообще не включают режим печки, и средняя температура торакса у них составляет около 35.5°С. Около 10% обсидчиков обеспечивают поддержание температуры расплода, включая режим печки, остальные - стабилизируют температуру просто своим присутствием. В этом тоже есть своя логика природы: откуда-то же надо брать эти
T_ambient = 34°C.Выводы
Главный результат этой заметки - не конкретные температуры, а демонстрация эмерджентного эффекта. В модели нет ни пчелы-регулятора, ни глобального плана обогрева. Каждая пчела следует одному простому правилу - локальному термостатированию по температуре воздуха под усиками, без перемещения по соту, обмена сигналами с соседями или любого другого поведения. Рассмотрена всего одна сторона одного сота.
Тем не менее из множества независимых локальных решений самоорганизуется глобальный результат: температура расплода стабилизируется в узком коридоре, необходимом для его развития. Ни одна пчела «не знает» о состоянии расплода в целом, каждая реагирует лишь на воздух в своей ячейке. Согласованность возникает через общую среду: расплод греет стенки, стенки греют воздух, воздух включает и выключает печки соседних пчёл. Это косвенная, стигмергическая координация - тот же принцип, что лежит в основе постройки сота и терморегуляции живой семьи.
Именно поэтому модель при всей её простоте воспроизводит ключевое свойство семьи как сверхорганизма: надёжный гомеостаз расплода без какого-либо централизованного управления.
🐝 Python-скрипт
"""
Один сот (сторона Y), расплод + грелки + обсидчики.
"""
import time, math
import numpy as np
from numba import njit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize
from matplotlib.patches import Patch, Polygon, Circle
import matplotlib
matplotlib.rcParams['font.size'] = 11
# ==========================================================
# ПАРАМЕТРЫ (копия из рабочего скрипта)
# ==========================================================
cols, rows = 76, 57
variant = 2
T_ambient = 34.0
CELL_EMPTY, CELL_BEE, CELL_BROOD = 0, 5, 6
I_HD, I_TH, I_AB, I_AI = 0, 1, 2, 3
LF5 = 203.0; tau_h = 1.0/0.064506; m_egg_mg = 0.000289; p_spec = 8.0; c_brood_spec = 3500.0
kf_mg_h = 0.12
_hours1 = np.arange(1, LF5); _hours2 = np.arange(LF5, 503)
_m1 = m_egg_mg * np.exp(_hours1 / tau_h)
_m2 = _m1[-1] - kf_mg_h * (_hours2 - LF5)
_hours_all = np.hstack((_hours1, _hours2)); _mass_all_mg = np.hstack((_m1, _m2))
def brood_power_W(age_h):
if age_h <= LF5 - 1: return m_egg_mg*np.exp(age_h/tau_h)*p_spec*1e-6
return (1e-5*age_h**2 - 0.005019*age_h + 0.877249)*1e-3
def brood_mass_kg(age_h):
return np.interp(np.clip(age_h,1.0,502.0), _hours_all, _mass_all_mg)*1e-6
cap_thickness = 0.15e-3; cap_density = 960.0
P_heater = 45e-3; T_sense_off = 36.0; T_sense_on = 35.0; T_hd_emerg = 47.0
pulse_surf = 1800.0; contact_boost = 1.5; g_contact = 150.0
brood_params = dict(center_col_frac=0.5, center_row_frac=0.5, radius_col_frac=0.30,
radius_row_frac=0.35, fill_rate=0.92, heater_fraction=0.30,
min_heaters=2, age_mean=300.0, age_span=36.0, age_noise=6.0, seed=42)
surf_params = dict(center_col_frac=0.5, center_row_frac=0.5, radius_col_frac=0.32,
radius_row_frac=0.37, density_max=0.6, seed=7)
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
C_cap=(A_entrance*cap_thickness*cap_density)*c_wax
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)
A_brood=math.pi*d_cell*h_total+math.pi*(d_cell/2)**2
h_in_brood=k_air/0.3e-3; R_br_ai=1.0/(h_in_brood*A_brood); G_conv_brood=h_in_brood*A_brood
se_brood=sigma*eps_wall*A_brood
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)
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)
cover_head=min((d_bee*h_head)/A_entrance,1.0); cover_thorax=min((d_bee*h_thorax)/A_entrance,1.0); cover_ab=min((d_bee*h_ab)/A_entrance,1.0)
G_cell_hd=h_entrance*cover_head*A_entrance; G_cell_th=h_entrance*cover_thorax*A_entrance; G_cell_ab=h_entrance*cover_ab*A_entrance
se_cell_hd=sigma*eps_wall*cover_head*A_entrance; se_cell_th=sigma*eps_wall*cover_thorax*A_entrance; se_cell_ab=sigma*eps_wall*cover_ab*A_entrance
G_contact_th=g_contact*cover_thorax*A_entrance
# ==========================================================
# ГЕНЕРАЦИЯ + ТОПОЛОГИЯ + ЯДРО (копия из рабочего скрипта)
# ==========================================================
def generate_brood_and_heaters(cols,rows,params):
rng=np.random.RandomState(params['seed'])
cmin=-(cols//2); cmax=cmin+cols-1; rmin=-(rows//2); rmax=rmin+rows-1
cx=cmin+params['center_col_frac']*cols; cy=rmin+params['center_row_frac']*rows
rx=params['radius_col_frac']*cols; ry=params['radius_row_frac']*rows
brood,gap=[],[]
for col in range(cmin,cmax+1):
for row in range(rmin,rmax+1):
ds=((col-cx)/rx)**2+((row-cy)/ry)**2
if ds<=1.0:
if rng.rand()0 and len(gap)0:
ang=[math.atan2(r-cy,c-cx) for (c,r) in gap]; order=np.argsort(ang); step=len(order)/nh
heaters=[gap[order[int(k*step)%len(order)]] for k in range(nh)]
else: heaters=[]
return brood,heaters,gap
brood_cells,heater_cells,gap_cells=generate_brood_and_heaters(cols,rows,brood_params)
N_heater = len(heater_cells)
print('Количество грелок ', N_heater)
R_ai_wall_arr=np.array([L_wall/(k_air*A) for A in [A_lateral]*6+[A_bottom]*3])
R_wall_self=np.array([L_wall/(lam_wax*A_lateral)]*6+[L_bottom/(lam_wax*A_bottom)]*3)
C_wall_base=np.array([C_wax_lateral]*6+[C_wax_bottom]*3)
cmin=-(cols//2); cmax=cmin+cols-1; rmin=-(rows//2); rmax=rmin+rows-1
cells=[]
for side in ["Y","^"]:
for col in range(cmin,cmax+1):
for row in range(rmin,rmax+1): cells.append((0,side,col,row))
cell_to_index={c:i for i,c 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: _="" age_arr="" age_h="" allset:="" allset="set(cells)" bl="" break="" brood_cells:="" brood_idx="[]" brood_power_w="" c_brood_arr="np.zeros(N_CELLS)" capped="1" cell="" cell_to_index="" cells="" col-1="" col="" dc="" deltas="EVEN_ROW_DELTAS" dr="" else:="" else="" enumerate="" even_row_deltas="" fid="" for="" get_neighbor_for_face="" heater_on="np.zeros(N_CELLS,np.int32);" i="" idx="" if="" in="" int="" is_capped="np.zeros(N_CELLS,np.int32)" k="" na="" nc="" neighbor_faces="" nr="" ns="" odd_row_deltas="" p_brood_arr="" p_extra="np.zeros(N_CELLS);" range="" return="" rf="" row-1="" row="" side="" states="" topology="">=LF5 else 0; is_capped[idx]=capped
C_brood_arr[idx]=brood_mass_kg(age_h)*c_brood_spec+(C_cap if capped else 0.0)
brood_idx.append(idx)
heater_idx=[]
for (col,row) in heater_cells:
idx=cell_to_index[(0,'Y',col,row)]; States[idx]=CELL_BEE; P_extra[idx]=P_heater; heater_on[idx]=1; heater_idx.append(idx)
def generate_surface_bees(cols,rows,side,params,States,cell_to_index):
rng=np.random.RandomState(params['seed'])
cmin=-(cols//2); cmax=cmin+cols-1; rmin=-(rows//2); rmax=rmin+rows-1
cx=cmin+params['center_col_frac']*cols; cy=rmin+params['center_row_frac']*rows
rx=params['radius_col_frac']*cols; ry=params['radius_row_frac']*rows
occupied=set(); bees=[]; bid=0
for col in range(cmin,cmax+1):
for row in range(rmin,rmax+1):
ds=((col-cx)/rx)**2+((row-cy)/ry)**2
# if ds<=1.0 and rng.rand()0.0: G=h_entrance*of*A_entrance; num+=Te*G; den+=G
if cover>0.0:
part=Cover_part[i]; is_contact=(part==1 and state==CELL_BROOD and is_capped[i]==1)
if not is_contact:
T_part=T_surf_bees[Cover_bee[i],part]; G=h_entrance*cover*A_entrance; num+=T_part*G; den+=G
if state==CELL_BROOD: num+=Tth/R_br_ai; den+=1.0/R_br_ai
Tai=num/den; T_next[i,I_AI]=Tai
isY=(i0.0:
part=Cover_part[i]; is_contact=(part==1 and state==CELL_BROOD and is_capped[i]==1)
if not is_contact:
T_part=T_surf_bees[Cover_bee[i],part]
p_rad+=sigma*eps_wall*(cover*A_entrance)*((T_part-t0)**4-(Tw-t0)**4)
for f in range(9):
Qa=(Tai-T_cur[i,4+f])/R_ai_wall_arr[f]; ni=Topology[i,f]
if ni!=-1: Qn=(T_cur[ni,4+Neighbor_Faces[i,f]]-T_cur[i,4+f])/R_wall_self[f]
else: Qn=(T_amb-T_cur[i,4+f])/R_wall_self[f]
psk=((A_lateral if f<6 a="" a_bottom="" cur="" dt="" else="" f="" grid="" i="" if="" min:="" n="" p_extra="" p_rad="" pr="Pmin" psk="" state="=CELL_BEE:" t_next="" wall_total="">0.0:
Tsense=Tai; on=heater_on[i]==1
woff=(Tsense>=T_sense_off) or (Thd>=T_hd_emerg)
won=(Tsense0.0 and Cover_part[i]==1 and is_capped[i]==1:
bee=Cover_bee[i]; Gc=G_contact_th*(contact_boost if surf_heater_on[bee]==1 else 1.0); Tc=T_surf_bees[bee,1]
Tnew=((Cbr/dt)*Tth+P_brood_arr[i]+G_conv_brood*Tai+Gc*Tc-Qrad)/(Cbr/dt+G_conv_brood+Gc)
T_next[i,I_TH]=Tnew; T_next[i,I_HD]=Tnew; T_next[i,I_AB]=Tnew
return Q_Y,Q_X
@njit(fastmath=True,cache=True)
def do_surface_bees_step(T_surf,T_cur,dt,T_amb,head_idx,th_idx,ab_idx,N_surf,States,is_capped,surf_heater_on,surf_heater_timer):
Qst=0.0; Qce=0.0; Qev=0.0; dU=0.0; maxres=0.0; Qst_arr=np.zeros(N_surf)
for k in range(N_surf):
hi=head_idx[k]; ti=th_idx[k]; ai=ab_idx[k]
Thd=T_surf[k,0]; Tth=T_surf[k,1]; Tab=T_surf[k,2]
Tsense=T_cur[hi,I_AI]
Twh=0.0; Twt=0.0; Twb=0.0
for f in range(9): Twh+=T_cur[hi,4+f]; Twt+=T_cur[ti,4+f]; Twb+=T_cur[ai,4+f]
Twh/=9.0; Twt/=9.0; Twb/=9.0
Tah=T_cur[hi,I_AI]; Tat=T_cur[ti,I_AI]; Tab_a=T_cur[ai,I_AI]
Pmet=(29.133-0.739*Tsense)*1e-3
if Pmet=T_sense_off) or (Thd>=T_hd_emerg)
won=(Tsensemaxres: maxres=res
T_surf[k,0]=Thd+(dt/C_hd)*dThd; T_surf[k,1]=Tth+(dt/C_th)*dTth; T_surf[k,2]=Tab+(dt/C_ab)*dTab
Qst+=Qst_c; Qce+=Qce_c; Qev+=Qev_c; dU+=dUk; Qst_arr[k]=Qst_c
return Qst,Qce,Qev,dU,maxres,Qst_arr
# ==========================================================
# ПРОГОН
# ==========================================================
dt_sim=0.1; time_total=3600.0
def make_buffers():
TA=np.full((N_CELLS,13),T_ambient); TB=np.full((N_CELLS,13),T_ambient)
for idx in heater_idx:
for buf in (TA,TB): buf[idx,I_HD]=Thd_i; buf[idx,I_TH]=Tth_i; buf[idx,I_AB]=Tab_i
hon=heater_on.copy()
Ts=np.tile(np.array([Thd_i,Tth_i,Tab_i]),(max(N_surf,1),1)).astype(np.float64)
sho=surf_heater_on.copy(); sht=surf_heater_timer.copy()
return TA,TB,hon,Ts,sho,sht
def run_simulation(use_surf,time_total,log=False):
TA,TB,hon,Ts,sho,sht=make_buffers()
Cf,Cb,Cp=(Cover_frac,Cover_bee,Cover_part) if use_surf else (Cover_frac0,Cover_bee0,Cover_part0)
steps=int(time_total/dt_sim); T_cur,T_next=TA,TB
log_t,log_Tbr,log_Tht,log_hon,log_TthS,log_QceS,log_surf_on=[],[],[],[],[],[],[]
maxres=0.0
for step in range(steps):
Q_Y,Q_X=do_cluster_step_surf(T_cur,T_next,Topology,Neighbor_Faces,States,C_grid,P_extra,
P_brood_arr,C_brood_arr,hon,is_capped,Cf,Cb,Cp,Ts,sho,T_entrance,dt_sim,T_ambient)
T_cur,T_next=T_next,T_cur
if use_surf and N_surf:
Qst,Qce,Qev,dU,mr,_=do_surface_bees_step(Ts,T_cur,dt_sim,T_ambient,surf_head_idx,surf_thorax_idx,surf_ab_idx,N_surf,States,is_capped,sho,sht)
Q_Y+=Qst
if mr>maxres: maxres=mr
if log and step%60==0:
log_t.append(step*dt_sim)
log_Tbr.append(T_cur[brood_idx, I_TH].mean() if brood_idx else T_ambient)
log_Tht.append(T_cur[heater_idx, I_TH].mean() if heater_idx else T_ambient)
log_hon.append(hon[heater_idx].mean() if heater_idx else 0.0)
log_TthS.append(Ts[:,1].mean() if N_surf else T_ambient)
log_QceS.append((Qce/N_surf*1000) if (use_surf and N_surf) else 0.0)
log_surf_on.append(sho[:N_surf].mean() if N_surf else 0.0)
T_air=np.full((cols,rows),T_ambient); T_br=np.full((cols,rows),np.nan)
for i,cell in enumerate(cells):
_,side,col,row=cell
if side!="Y": continue
ci,ri=col-cmin,row-rmin; T_air[ci,ri]=T_cur[i,I_AI]
if States[i]==CELL_BROOD: T_br[ci,ri]=T_cur[i,I_TH]
return dict(T_cur=T_cur,Ts=Ts,hon=hon,sho=sho,T_air=T_air,T_br=T_br,
log_t=log_t,log_Tbr=log_Tbr,log_Tht=log_Tht,log_hon=log_hon,
log_TthS=log_TthS,log_QceS=log_QceS,log_surf_on=log_surf_on,maxres=maxres)
# компиляция + прогоны
TA0,TB0,hon0,Ts0,sho0,sht0=make_buffers()
do_cluster_step_surf(TA0,TB0,Topology,Neighbor_Faces,States,C_grid,P_extra,P_brood_arr,C_brood_arr,
hon0,is_capped,Cover_frac,Cover_bee,Cover_part,Ts0,sho0,T_entrance,dt_sim,T_ambient)
if N_surf: do_surface_bees_step(Ts0,TA0,dt_sim,T_ambient,surf_head_idx,surf_thorax_idx,surf_ab_idx,N_surf,States,is_capped,sho0,sht0)
print("Baseline..."); res_base = run_simulation(use_surf=False, time_total=time_total, log=True)
duty_cycle = np.mean(res_base['log_hon']) * 100
print(f" Доля времени грелок в печке: {duty_cycle:.0f}%")
print("Основной..."); res_surf=run_simulation(True,time_total,True)
T_cur=res_surf['T_cur']; T_surf_bees=res_surf['Ts']; heater_on=res_surf['hon']; surf_heater_on=res_surf['sho']
# --- Прогон: только расплод (без грелок, без обсидчиков) ---
print("Только расплод (без пчёл)...")
States_save2=States.copy(); P_extra_save2=P_extra.copy(); heater_on_save2=heater_on.copy()
for idx in heater_idx:
States[idx]=CELL_EMPTY; P_extra[idx]=0; heater_on[idx]=0
res_brood_only=run_simulation(use_surf=False, time_total=time_total, log=False)
States[:]=States_save2; P_extra[:]=P_extra_save2; heater_on[:]=heater_on_save2
T_br_only=res_brood_only['T_br']
print(f" T расплода без пчёл: {np.nanmin(T_br_only):.3f}..{np.nanmax(T_br_only):.3f}°C")
# ==========================================================
# ВИЗУАЛИЗАЦИЯ (в столбик, для блога)
# ==========================================================
FS=13; size=1.0; width_hex=math.sqrt(3)*size; dy_hex=1.5*size
n_cols,n_rows=cols,rows
def hex_xy(ci,ri):
col=ci+cmin; row=ri+rmin; return width_hex*(col+0.5*(row&1)), -dy_hex*row
def extent():
xs=[]; ys=[]
for ci in range(n_cols):
for ri in range(n_rows): x,y=hex_xy(ci,ri); xs.append(x); ys.append(y)
return min(xs)-size,max(xs)+size,min(ys)-size,max(ys)+size
def cbar_inset(ax,fig,cmap,vmin,vmax,x0,x1,y0,y1,lab):
sx=x1-x0; gap=sx*0.02; cw=sx*0.035
cax=ax.inset_axes([x1+gap,y0,cw,y1-y0],transform=ax.transData)
sm=ScalarMappable(cmap=cmap,norm=Normalize(vmin,vmax)); sm.set_array([])
cb=fig.colorbar(sm,cax=cax); cb.ax.tick_params(labelsize=FS-1); cb.set_label(lab,fontsize=FS); return cw,gap
def draw_hex(ax,fig,field,cmap,vmin,vmax,title='',cb='',brood_field=None):
co=plt.get_cmap(cmap); ax.set_aspect('equal')
for ci in range(n_cols):
for ri in range(n_rows):
x,y=hex_xy(ci,ri); cx_,cy_=[],[]
for kk in range(6):
ang=math.pi/180*(60*kk-30); cx_.append(x+size*math.cos(ang)); cy_.append(y+size*math.sin(ang))
cx_.append(cx_[0]); cy_.append(cy_[0]); val=field[ci,ri]
if np.isnan(val):
ax.fill(cx_,cy_,color='#f0f0f0',edgecolor='#cccccc',lw=0.3); continue
norm=np.clip((val-vmin)/max(vmax-vmin,1e-5),0,1); ec,ew='#888888',0.3
if brood_field is not None and not np.isnan(brood_field[ci,ri]): ec,ew='#228B22',1.2
ax.fill(cx_,cy_,color=co(norm),edgecolor=ec,lw=ew)
x0,x1,y0,y1=extent(); cw,gap=cbar_inset(ax,fig,co,vmin,vmax,x0,x1,y0,y1,cb)
ax.set_xlim(x0-0.3,x1+gap+cw+0.3); ax.set_ylim(y0-0.3,y1+0.3)
ax.set_title(title,fontsize=FS,pad=8); ax.axis('off')
def draw_heater_circles(ax,cmap_obj,norm_obj):
r_circle=(d_bee/thread)*size
for idx in heater_idx:
_,_,col,row=cells[idx]; ci,ri=col-cmin,row-rmin; x,y=hex_xy(ci,ri)
Tab_val=T_cur[idx,I_AB]; on=heater_on[idx]==1
circ=Circle((x,y),radius=r_circle,facecolor=cmap_obj(norm_obj(Tab_val)),
edgecolor='red' if on else 'orange',linewidth=1.2,zorder=8,alpha=0.95)
ax.add_patch(circ)
def draw_bee_segments(ax,cmap_obj,norm_obj):
s_map=width_hex/5.4; L_hd,L_th,L_ab=1.5*s_map,4.5*s_map,6.0*s_map
W=d_bee*1000*s_map; gap=1.3*s_map
for k in range(N_surf):
_,_,tc,tr=cells[surf_thorax_idx[k]]; _,_,hc,hr=cells[surf_head_idx[k]]; _,_,ac,ar=cells[surf_ab_idx[k]]
px_th=hex_xy(tc-cmin,tr-rmin); px_hd=hex_xy(hc-cmin,hr-rmin); px_ab=hex_xy(ac-cmin,ar-rmin)
dxu=px_hd[0]-px_ab[0]; dyu=px_hd[1]-px_ab[1]; ln=math.hypot(dxu,dyu)
if ln<1e-9: ---="" -0.01="" 0.04="" 0.5="" 0="" 1.2:="" 1.3:="" 100="" 12="" 1:="" 1="" 35.5="" 3:="" 3="" 4:="" 5:="" 9="" _="" a1.axhline="" a1.grid="" a1.legend="" a1.plot="" a1.set_title="" a1.set_ylabel="" a1="" a2.axhline="" a2.grid="" a2.legend="" a2.plot="" a2.set_title="" a2.set_ylabel="" a2="" a3.grid="" a3.legend="" a3.plot="" a3.set_title="" a3.set_xlabel="" a3.set_ylabel="" a3="" ab="" age_arr="" age_y="" air="" allv="np.concatenate(vals_all);" allv_base.max="" allv_base="np.concatenate(vals_base)" alpha="0.7," ambient="" ax.add_patch="" ax.legend="" ax1="" ax2="" ax3="" ax="" bbox_inches="tight" bbox_to_anchor="(0.5,-0.01)," bmax="" bmin="" br="" br_max="" br_min="" br_y="" brood_arr="" brood_field="T_br_Y)" brood_idx="" cb="°C" cell="" cells="" ci="" circ="" cmap_g="" cmap_obj="" col="" color="gray" continue="" cur="" cx="" cy="" def="" dmax="" dpi="150," draw_bee_segments="" draw_heater_circles="" draw_hex="" dt="" dt_br="" duty_heater_base="np.mean(res_base[" duty_heater_surf="np.mean(res_surf[" duty_surf="np.mean(res_surf[" edgecolor="black" else="" enumerate="" ext="extent();" f="" facecolor="none" fh2="figsize_for(*ext," fh="" fig1_brood_age_power_t.png="" fig3_brood_dt.png="" fig4_dynamics.png="" fig5_dynamics.png="" fig="" fig_heaters_and_surface.png="" fig_heaters_only.png="" figsize="(9," float="" fontsize="FS-1);" for="" frameon="False)" from="" fw2="" fw="10;" g--="" g-="" gap="" gmax="" gmin="" handles="legend_full()," hd="" hdr="f" heater_idx:="" heater_idx="" heater_on="" hon="" hot="" i="" i_ab="" i_th="" idx="" if="" import="" in="" inferno="" k="" label="обсидчики" legend_full="" legend_heaters="" len="" line2d="" linestyle="none" linewidth="1.2," loc="upper center" log="" log_b="" log_hon="" log_s="" log_surf_on="" log_t="" log_tbr="" log_tht="" log_tths="" lorrd="" ls=":" lw="2," m-="" magma="" marker="o" markeredgecolor="orange" markeredgewidth="1.5," markerfacecolor="#bdbdbd" markersize="9," matplotlib.lines="" min="" n="" n_surf:="" ncol="4," norm_b="" norm_g="" norm_obj="" np.any="" np.array="" np.isnan="" np.nanmax="" np.nanmin="" on="" orange="" p1="(cx+ux*L/2+vx*W/2,cy+uy*L/2+vy*W/2);" p2="(cx+ux*L/2-vx*W/2,cy+uy*L/2-vy*W/2)" p3="(cx-ux*L/2-vx*W/2,cy-uy*L/2-vy*W/2);" p4="(cx-ux*L/2+vx*W/2,cy-uy*L/2+vy*W/2)" p="" p_br_y="" patch="" plasma="" plt.savefig="" plt.show="" plt.subplots="" plt.tight_layout="" poly="" print="" pulse_surf="" px_th="" r--="" r-="" r_circle="(d_bee/thread)*size" radius="r_circle," ravel="" rect="[0,0.06,1,1])" res_base="" return="" ri="col-cmin,row-rmin;" row="cells[idx];" segs:="" segs="[(px_th,L_th,T_surf_bees[k,1])," sharex="True)" side="" states="" surf="" surf_bees.ravel="" surf_bees="" t="" t_abdomen.="" t_air_y="" t_amb="{T_ambient}°C" t_ambient="" t_br_base="res_base[" t_br_only="" t_br_y="res_surf[" t_brood="" t_cur="" tab_val="Tcur[idx,I_AB];" tb_base="res_base[" tb_only="res_brood_only[" tb_surf="res_surf[" tcur="" th="" tight="" title="Обогрев обсидчиками (dT от базового)" ux="" uy="" vals_all.append="" vals_all="[T_air_Y.ravel()]" vals_base.append="" vals_base="[res_base[" vx="" vy="-uy,ux" x="" y="hex_xy(ci,ri)" zorder="8,">6} {'мин':>6} {'макс':>6} {'ΔT':>6} {'≥34.5':>6} {'>38':>4}"
print(hdr); print(f" {'-'*66}")
for name, Tb in [('Без пчёл', Tb_only),
('Только грелки ', Tb_base),
('Грелки+обсидчики ', Tb_surf)]:
print(f" {name:<28 b.mean="">6.2f} {Tb.min():>6.2f} {Tb.max():>6.2f} "
f"{Tb.mean()-T_ambient:>+6.2f} {int((Tb>=34.5).sum()):>6} {int((Tb>38).sum()):>4}")
print(f" {'-'*66}")
print(f" Доля времени в режиме печки:")
print(f" грелки (без обсидчиков): {duty_heater_base:.0f}%")
print(f" грелки (с обсидчиками): {duty_heater_surf:.0f}%")
print(f" обсидчики: {duty_surf:.0f}%")
print(f"{'='*72}")
Комментарии
Отправить комментарий