Введение
Улочка - это пространство между двумя соседними сотами (или сотом и стенкой улья) шириной 8–12 мм (в модели задано d_street = 12 мм, как параметр для варьирования). Воздух в улочке переносит тепло между сотами, а при наличии градиента температуры возникает конвекция, но здесь ее пока не учитываем - это следующий шаг развития модели.
Собственно подход к формализации улочки для модели описан в заметке Улочка без пчел и без конвекции, здесь кратко повторим.
Улочка моделируется гексагональной сеткой узлов, по одному на каждую пару ячеек "смотрящих" друг на друга (76×57 = 4332 узла для рамки Дадана). Каждый узел имеет 8 связей:
- 2 ячейки (входы): ячейка
(0,'^',c,r)сота 0 и ячейка(1,'Y',c,r)сота 1 - обмен теплом через открытые входы (проводимостьG_entr) - 6 соседей (гексагональная сетка) - теплопроводность воздуха (проводимость
G_lat)
Параметры узла:
| Параметр | Формула | Значение |
|---|---|---|
A_node (площадь шестиугольника) |
(√3/2)·thread² |
25.3 мм² |
V_node (объем узла) |
A_node·d_street |
0.30 мм³ |
C_node (теплоемкость узла) |
ρ_air·V_node·c_air |
3.7×10⁻⁴ Дж/К |
G_entr(тепловая проводимость к ячейкам) |
h_entr·A_entr |
2.6×10⁻⁴ Вт/К |
G_lat(тепловая проводимость к соседям) |
k_air·d_street/√3 |
1.7×10⁻⁴ Вт/К |
C_node мала (воздух лёгкий), G_entr и G_lat одного порядка. Это значит, что узел улочки быстро реагирует на изменение потоков (малая инерция), а тепло распределяется между соседями сопоставимо с притоком от ячеек.
Периметр сетки (крайние узлы) фиксирован на T_ambient (условие Дирихле) - по краям улочка сообщается с общим объёмом улья, в котором пока температуру считаем постоянной. В дальнейшем (при моделировании улья целиком) этот параметр также будет меняться.
Уравнение теплообмена узла (явный Эйлер):
$$ C_{node}\frac{dT_s}{dt}=G_{entr}(T_{air,0}−T_s)+G_{entr}(T_{air,1}−T_s)+∑_{j=1}^6G_{lat}(T_{s,j}−T_s)$$
Текст программы, как обычно, в спойлере в конце заметки.
1. Эффект улочки
До сих пор модель считала (все результаты в заметке Как пчелы греют расплод), что по ту сторону входа ячейки находится бесконечный резервуар с постоянной температурой T_ambient = 34°C. Это верно для одиночного сота вне улья, но не для гнезда, где соты стоят плотно, через 8-12 мм друг от друга. Воздух улочки нагревается от обоих образующих её сотов, и его температура переменная, а не константа. Причем она не одинакова в разных узлах улочки, т.е. это массив температур T_street[col, row]:
Полный адрес ячейки сота состоит из номера сота, стороны , столбца и строки: frame, side, col, row, поэтому в расчете принимают участие N_CELLS = 4 × 76 × 57 = 17 328 ячеек сот.
Правило включения/выключения режима печки оставляем прежним :
- измеренная температура ниже 35°С - включить режим печки;
- измеренная температура выше 36°С - выключить режим печки;
- температура головы выше 47°С - аварийное выключение режима печки.
Порядок вычислений на каждом шаге:
do_cluster_step()— обновление всех 17 328 ячеек (T_entranceберётся изT_streetдля внутренних сторон, а изT_ambient- для внешних). Возвращает массивQ_entranceдля каждой ячейки.do_street_step()— обновление 4332 узлов улочки (используетQ_entranceот двух ячеек + теплопроводность к соседям).
Первый результат (расплод на одной стороне 1313 зародышей и пустой сот напротив через улочку 12 мм):
Видим, что температура в улочке напротив расплода повысилась примерно на пол градуса, а воздух в ячейках противоположной стороны - примерно на 0.2°С. Но, самое интересное, примерно на столько же (0.19°С) повысилась и средняя температура расплода: было (среднее/минимальное/максимальное) 34.40/34.22/34.51 (см. предыдущую заметку) стало 34.59/34.33/34.77. Т.е. простая постановка напротив расплода пустой рамки дает существенный эффект - эффект улочки.
Второй результат (расплод в 1313 ячейках и грелки в 33 - на одной стороне напротив пустого сота):
- тепловая карта
- динамика переходного процесса
Видим, что стационарный режим в улочке устанавливается примерно в течение часа (3600с), причем в начале этого часа все пчелы в ячейках вынуждены включать режим печки - это к вопросу о стрессе при осмотре пчеловодом рамок с расплодом. Температура расплода 34.87/34.40/35.16 (было 34.72/34.29/35.11 см. предыдущую заметку), т.е. средняя температура расплода выше, чем при отсутствии пустого сота на 0.15°С.
В центральной части (напротив расплода) температура воздуха в улочке близка к температуре расплода, поэтому и "потери" в улочку составляют всего 1.1 мВт.
В установившемся режиме (после 3600 секунд) включают режим печки в среднем около 10% пчел со средней продолжительностью 42 секунды. На одну пчелу в среднем за два часа пришлось 10 включений.
В гнезде не бывает, чтобы напротив сотов с расплодом был пустой сот, это может произойти только при вмешательстве пчеловода. Поэтому предыдущие оценки нужны были лишь для того, чтобы "пощупать", почувствовать, как возникает "эффект улочки". Далее мы поставим не пустой сот, а сот с расплодом и пчелами.
Рассмотрим два сценария: 1) расплод на двух сторонах плюс пчелы-грелки в пустых ячейках среди расплода 2) к первому сценарию добавим пчел-обсидчиков.
2. Расплод и пчелы в ячейках
Запускаем модель по первому сценарию: расплод на двух сторонах и пчелы-грелки в пустых ячейках среди расплода
Температура расплода составила 34.96/34.51/35.14. Средняя температура расплода оказалась опять близка к температуре улочки в центральной части, поэтому тепловой обмен между сотами и улочкой остался на небольшом уровне 1.9 мВт. Доля "включенных" грелок сократилась примерно вдвое (по сравнению со случаем одного пустого сота), среднее время работы осталось на уровне 42 секунд, но при этом на одну пчелу в среднем пришлось 5.2 включения на 2 часа.
3. Расплод на двух сторонах + пчелы-грелки в ячейках + обсидчики
Запускаем модель по второму сценарию - добавляем обсидку: 339 + 338 = 677 пчел.
Видим, что средняя температура расплода возросла примерно на пол градуса и составила 35.49/35.02/35.74. При этом грелки в ячейках вообще не включаются - нет необходимости, а, точнее, температура в районе усиков никогда не опускается ниже 35°С. Доля обсидчиков, находящихся в режиме печки, в среднем составляет около 1%. Средняя продолжительность включения составляет примерно 14 секунд, а за два часа пчела включается всего 3.2 раза.
Но обсидка создает еще и "эффект улочки 2": суммарный поток тепла через входы ячеек меняет знак - тепло идет не из ячеек, а в ячейки, и его величина значительна -212 мВт. Теплое одеяло обсидки греет расплод, как четыре-пять включенных грелок, но при этом гораздо равномернее и "мягче".
Выводы
Улочка не является пассивным пустым пространством между сотами. Это активная термодинамическая среда, которая связывает соты в единую систему и существенно влияет на тепловой баланс расплода. В предыдущей модели воздух за входом ячейки считался бесконечным резервуаром с постоянной температурой; теперь его температура переменная, и это меняет всю картину.
Эффект улочки проявляется в трёх ипостасях, по нарастающей:
- Тепловой буфер. Даже простой пустой сот напротив расплода повышает его среднюю температуру на 0.19°C (34.40 → 34.59°C): улочка становится прослойкой прогретого воздуха вместо холодного стока.
- Взаимный обогрев сотов. Два сота с расплодом греют друг друга через улочку. При том же правиле обогрева доля включённых грелок вдвое ниже, чем в случае пустого сота напротив - соты «помогают» друг другу.
- Активный обогрев обсидкой. Это главный эффект. Обсидчики, сидящие на сотах, даже без включения режима печки прогревают воздух улочки дорсально, и она превращается в «тепловое одеяло». Суммарный поток через входы ячеек меняет знак: тепло идёт не из ячеек в улочку, а из улочки в ячейки -212 мВт, что эквивалентно примерно пяти включённым грелкам, но распределённым равномерно, без локальных перегревов. Средняя температура расплода выходит на цель 35.5°C, а максимум не превышает 35.7°C, т.е. температура в установившемся режиме поддерживается очень точно практически без колебаний.
Как уже было отмечено в предыдущей заметке, здесь нет ни пчелы-регулятора, ни глобального плана обогрева. Каждая пчела следует одному локальному правилу - термостатированию по температуре воздуха под усиками. Но теперь координация происходит через более протяжённую общую среду: пчела греет воздух улочки, этот воздух доходит до соседних ячеек и узлов улочки, что влияет на решения других пчёл.
Показательна и «экономическая» сторона: при наличии обсидки грелки в ячейках вообще не включают режим печки - воздух в районе их усиков не опускается ниже 35°C. Система самоорганизуется так, что основную работу берёт на себя тот слой пчёл, который выполняет её эффективнее (обсидчики, греющие через крышечки), а грелки в пустых ячейках становятся резервом. Более того, как замечено в [32], молодые пчелы (до 2 дневного возраста), которые еще не способны работать печками, часто занимают пустые ячейки - видимо, приходят погреться в яслях и поспать. Оптимальное распределение ролей возникает само, без какого-либо управления.
🐝 Python-скрипт
"""
Архитектура:
Frame 0: (0,'Y') внешняя, (0,'^') внутренняя (к улочке)
Frame 1: (1,'Y') внутренняя (к улочке), (1,'^') внешняя
Узлы улочки: T_street[cols,rows], 8 связей (2 ячейки + 6 соседей).
Граница улочки: T_ambient (Дирихле).
Термостат грелок: по усикам (T_sense = воздух).
Масса эмбриона: формула (экспонента + линейное снижение).
Конфигурация запуска задается ключами BROOD_SIDES, ENABLE_HEATERS, ENABLE_OBSIDCHIKI
"""
import time
import 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
# ==========================================================
# 1. ПАРАМЕТРЫ
# ==========================================================
cols, rows = 76, 57
# cols, rows = 7, 7 # отладка
variant = 2
T_ambient = 35.0
CELL_EMPTY, CELL_BEE, CELL_BROOD = 0, 5, 6
I_HD, I_TH, I_AB, I_AI = 0, 1, 2, 3
# --- Улочка ---
d_street = 12e-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))
# --- конфигурация шага ---
BROOD_SIDES = [(0,'^'), (1,'Y')] # обе внутренние стороны
# BROOD_SIDES = [(0,'^')]
# ENABLE_HEATERS = False # Шаг 1: без грелок
ENABLE_HEATERS = True
ENABLE_OBSIDCHIKI = True
T_sense_off = 36.0; T_sense_on = 35.0; T_hd_emerg = 47.0
print(f">>> Конфигурация: расплод на {len(BROOD_SIDES)} сторонах {BROOD_SIDES}")
print(f" грелки={'ВКЛ' if ENABLE_HEATERS else 'выкл'}, "
f"обсидчики={'ВКЛ' if ENABLE_OBSIDCHIKI else 'выкл'}")
print(f" термостат: вкл<{T_sense_on}°C, выкл≥{T_sense_off}°C, авария {T_hd_emerg}°C")
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
# --- Расплод (параметры генерации) ---
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)
# --- Пчела ---
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
pulse_surf = 1800.0; contact_boost = 1.; g_contact = 100.0
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)
# Дорсальные площади и сопротивления (конвекция в воздух улочки)
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
# --- Параметры улочки ---
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)
print(f">>> T_amb={T_ambient}°C, рамка {cols}×{rows}, d_street={d_street*1e3:.0f} мм")
print(f" C_node={C_node*1e4:.3f}×10⁻⁴ Дж/К, G_entr={G_entr*1e4:.3f}×10⁻⁴ Вт/К, G_lat={G_lat*1e4:.3f}×10⁻⁴ Вт/К")
print(f" термостат: T_sense_on={T_sense_on}, T_sense_off={T_sense_off}, T_hd_emerg={T_hd_emerg}")
# ==========================================================
# 2. ГЕНЕРАЦИЯ РАСПЛОДА + ГРЕЛОК (на стороне (0,'^'))
# ==========================================================
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
# ==========================================================
# 3. ТОПОЛОГИЯ (4 стороны: frame 0 Y,^; frame 1 Y,^)
# ==========================================================
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 frame in [0, 1]:
for side in ["Y","^"]:
for col in range(cmin,cmax+1):
for row in range(rmin,rmax+1):
cells.append((frame,side,col,row))
cell_to_index={c:i for i,c in enumerate(cells)}; N_CELLS=len(cells)
n_cols, n_rows = cols, rows
N_street = n_cols * n_rows
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(frame, side, col, row, face_id):
if face_id<6: ---="" 0="" 42="" age_arr="" age_h="" allset:="" allset="set(cells)" bc:="" bc="" bl="" break="" brood_heaters_per_side="" brood_idx="[]" brood_params.get="" brood_power_w="" brood_sides:="" 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="" enable_heaters:="" enumerate="" even_row_deltas="" fid="" for="" fr="" frame="" gc="generate_brood_and_heaters(cols," get_neighbor_for_face="" hc="brood_heaters_per_side[(fr,side)]" heater_on="np.zeros(N_CELLS,np.int32);" i="" idx="" if="" in="" int="" is_capped="np.zeros(N_CELLS,np.int32)" k="" min_heaters="" na="" nc="" neighbor_faces="" nf="" not="" nr="" ns="" odd_row_deltas="" p_brood_arr="" p_extra="np.zeros(N_CELLS);" params_side="" range="" return="" rf="" row-1="" row="" rows="" seed="" 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=[]
if ENABLE_HEATERS:
for (fr,side) in BROOD_SIDES:
bc, hc = brood_heaters_per_side[(fr,side)]
for (col,row) in hc:
idx=cell_to_index[(fr,side,col,row)]; States[idx]=CELL_BEE
P_extra[idx]=P_heater; heater_on[idx]=1; heater_idx.append(idx)
# --- разбивка по сторонам ---
brood_per_side = []
heater_per_side = []
for (fr, side) in BROOD_SIDES:
bc, hc = brood_heaters_per_side[(fr, side)]
brood_per_side.append(f"({fr},'{side}')={len(bc)}")
heater_per_side.append(f"({fr},'{side}')={len(hc) if ENABLE_HEATERS else 0}")
n_capped = int(is_capped[brood_idx].sum()) if brood_idx else 0
print(f">>> Расплод: {', '.join(brood_per_side)}, всего={len(brood_idx)} (запечатано={n_capped})")
print(f">>> Грелки: {', '.join(heater_per_side)}, всего={len(heater_idx)}, "
f"ENABLE_HEATERS={ENABLE_HEATERS}")
def generate_surface_bees(cols, rows, frame, side, params, bee_cell_set_BEE):
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,'^')={len(surf_bees_0)}, сторона (1,'Y')={len(surf_bees_1)}, всего={N_surf}")
# ==========================================================
# 4. УЗЛЫ УЛОЧКИ
# ==========================================================
# Маппинг: ячейка -> индекс узла улочки
# Внутренние стороны: (0,'^') и (1,'Y') -> street_idx >= 0
# Внешние стороны: (0,'Y') и (1,'^') -> street_idx = -1
street_idx = np.full(N_CELLS, -1, np.int32)
street_cell0 = np.zeros(N_street, np.int32) # индекс ячейки (0,'^',c,r)
street_cell1 = np.zeros(N_street, np.int32) # индекс ячейки (1,'Y',c,r)
for ci in range(n_cols):
for ri in range(n_rows):
col = ci + cmin; row = ri + rmin
si = ci * n_rows + ri
idx0 = cell_to_index[(0,'^',col,row)]
idx1 = cell_to_index[(1,'Y',col,row)]
street_cell0[si] = idx0
street_cell1[si] = idx1
street_idx[idx0] = si
street_idx[idx1] = si
# Соседи узлов (гексагональная сетка, топология стороны Y)
street_neighbors = np.full((N_street, 6), -1, np.int32)
is_boundary = np.zeros(N_street, dtype=np.bool_)
for ci in range(n_cols):
for ri in range(n_rows):
si = ci * n_rows + ri
row = ri + rmin
deltas = ODD_ROW_DELTAS if (row%2==1) else EVEN_ROW_DELTAS
for j,(dc,dr) in enumerate(deltas):
ni, nj = ci+dc, ri+dr
if 0<=ni>> Улочка: {n_cols}×{n_rows} узлов, внутренних={n_interior}, граничных={int(np.sum(is_boundary))}")
# ==========================================================
# 5. НАЧАЛЬНЫЕ T И МАТРИЦЫ
# ==========================================================
def calc_init(T_air_val, P_extra_val=0.0):
l=(1.3493*np.exp(-0.058*T_air_val))*1e-3; Rthd=l/(k_hd*s_r); Rtab=l/(k_ab*s_r)
def eq(v):
h,t,b=v; q=ke0*t; Pr=max((29.133-0.739*T_air_val)*1e-3,Pmin)
ph=-eps*sigma*((h-t0)**4-(T_air_val-t0)**4)*A_hd; pt=-eps*sigma*((t-t0)**4-(T_air_val-t0)**4)*A_th
pb=-eps*sigma*((b-t0)**4-(T_air_val-t0)**4)*A_ab
return [ph+(t-h)/Rthd-(h-T_air_val)/R_head_air_open-q,
Pr+P_extra_val+pt-(t-h)/Rthd-(t-b)/Rtab-(t-T_air_val)/R_thorax_air_open,
pb+(t-b)/Rtab-(b-T_air_val)/R_ab_air_open-q]
return fsolve(eq,[T_air_val+1,T_air_val+2,T_air_val+1],xtol=1e-8)
Thd_i,Tth_i,Tab_i=calc_init(T_ambient,0.0)
C_grid=np.zeros((N_CELLS,13))
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
elif States[i]==CELL_BROOD: C_grid[i,I_TH]=C_brood_arr[i]; C_grid[i,I_AI]=C_air_cell
else: C_grid[i,I_AI]=C_air_cell
if N_surf > 0:
surf_thorax_idx = np.array([cell_to_index[b[1]] for b in surf_bees], np.int32)
surf_head_idx = np.array([cell_to_index[b[2]] for b in surf_bees], np.int32)
surf_ab_idx = np.array([cell_to_index[b[3]] for b in surf_bees], np.int32)
surf_node_head = street_idx[surf_head_idx]
surf_node_thorax = street_idx[surf_thorax_idx]
surf_node_ab = street_idx[surf_ab_idx]
else:
surf_thorax_idx = surf_head_idx = surf_ab_idx = np.zeros(0, np.int32)
surf_node_head = surf_node_thorax = surf_node_ab = np.zeros(0, np.int32)
# Cover-массивы: какая часть тела какого обсидчика перекрывает ячейку
Cover_frac = np.zeros(N_CELLS)
Cover_frac0=np.zeros(N_CELLS)
Cover_bee0 = np.full(N_CELLS, -1, np.int32)
Cover_part0 = np.zeros(N_CELLS, np.int32)
Cover_bee = np.full(N_CELLS, -1, np.int32)
Cover_part = np.zeros(N_CELLS, np.int32)
for k in range(N_surf):
Cover_frac[surf_head_idx[k]]=cover_head; Cover_bee[surf_head_idx[k]]=k; Cover_part[surf_head_idx[k]]=0
Cover_frac[surf_thorax_idx[k]]=cover_thorax; Cover_bee[surf_thorax_idx[k]]=k; Cover_part[surf_thorax_idx[k]]=1
Cover_frac[surf_ab_idx[k]]=cover_ab; Cover_bee[surf_ab_idx[k]]=k; Cover_part[surf_ab_idx[k]]=2
# Температуры частей тела + состояние печки обсидчиков
T_surf_bees = np.tile(np.array([Thd_i, Tth_i, Tab_i]), (N_surf, 1)).astype(np.float64)
surf_heater_on = np.zeros(N_surf, np.int32)
surf_heater_timer = np.zeros(N_surf, np.float64)
# ==========================================================
# 6. ЯДРО: СОТ (Te из T_street) + УЗЛЫ УЛОЧКИ
# ==========================================================
@njit(fastmath=True, cache=True)
def do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid,
P_extra, P_brood_arr, C_brood_arr, heater_on, is_capped,
Cover_frac, Cover_bee, Cover_part, T_surf_bees, surf_heater_on,
T_street, street_idx, dt, T_amb):
Q_entrance_arr = np.zeros(len(T_cur))
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]
si=street_idx[i]
Te = T_street[si] if si>=0 else T_amb
# --- T воздуха ---
num=0.0; den=0.0
for f in range(9):
num+=T_cur[i,4+f]/R_ai_wall_arr[f]; den+=1.0/R_ai_wall_arr[f]
if state==CELL_BEE:
num+=Tth/R_th_ai+Thd/R_hd_ai+Tab/R_ab_ai+Te/R_ab_
den+=1.0/R_th_ai+1.0/R_hd_ai+1.0/R_ab_ai+1.0/R_ab_
else:
cover = Cover_frac[i]
of = 1.0 - cover
if of > 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
# --- поток через вход (в узел улочки или в T_amb) ---
if state==CELL_BEE:
Qtot=hA_*(Tai-Te)+(Tab-Te)/R_ab_
else:
Qtot = h_entrance * (1.0 - Cover_frac[i]) * A_entrance * (Tai - Te)
Q_entrance_arr[i]=Qtot
# --- средняя T стенок ---
Tw=0.0
for f in range(9): Tw+=T_cur[i,4+f]
Tw/=9.0
# --- излучение ---
p_rad=0.0
if state==CELL_BEE:
p_s_th=-se_th*((Tth-t0)**4-(Tw-t0)**4); p_s_hd=-se_hd*((Thd-t0)**4-(Tw-t0)**4)
p_s_ab=-se_ab_*((Tab-t0)**4-(Tw-t0)**4); p_rad=p_s_th+p_s_hd+p_s_ab
if state==CELL_BROOD:
p_rad=se_brood*((Tth-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=(Tsense 0.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_entrance_arr
@njit(fastmath=True, cache=True)
def do_street_step(T_street, T_street_new, Q_entrance_arr, Q_dorsal_node,
street_cell0, street_cell1, street_neighbors,
is_boundary, dt, T_amb, C_node_val, G_lat_val):
for si in range(len(T_street)):
if is_boundary[si]:
T_street_new[si]=T_amb; continue
Q=Q_entrance_arr[street_cell0[si]]+Q_entrance_arr[street_cell1[si]]+Q_dorsal_node[si]
Ts=T_street[si]
for j in range(6):
nb=street_neighbors[si,j]
if nb>=0: Q+=G_lat_val*(T_street[nb]-Ts)
T_street_new[si]=Ts+(dt/C_node_val)*Q
@njit(fastmath=True, cache=True)
def do_surface_bees_step(T_surf, T_cur, T_street, dt,
head_idx, th_idx, ab_idx, node_head, node_th, node_ab,
N_surf, States, is_capped, surf_heater_on, surf_heater_timer):
Q_dorsal_node = np.zeros(len(T_street))
Qst=0.0; Qce=0.0; Qev=0.0; dU=0.0; maxres=0.0
for k in range(N_surf):
hi=head_idx[k]; ti=th_idx[k]; ai=ab_idx[k]
nhi=node_head[k]; nti=node_th[k]; nai=node_ab[k]
Thd=T_surf[k,0]; Tth=T_surf[k,1]; Tab=T_surf[k,2]
Tsense=T_cur[hi,I_AI]
# Температура улочки под каждой частью тела
Tst_hd=T_street[nhi]; Tst_th=T_street[nti]; Tst_ab=T_street[nai]
# Средние температуры стенок под частями тела
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
return Q_dorsal_node, Qst, Qce, Qev, dU, maxres
# ==========================================================
# 7. ПРОГОН
# ==========================================================
dt_sim=0.1; time_total=7200.0
# ==========================================================
# БАЛАНС МОЩНОСТИ: термогенез vs потери во вне
# ==========================================================
def calc_power_balance(T_cur, T_prev, T_s, hon, sho, dt_window, Q_dorsal_sum=0.0):
"""Термогенез (расплод+грелки) и потери во вне от каждого сота + периметр улочки."""
# --- термогенез расплода ---
P_brood = float(P_brood_arr[brood_idx].sum())
# --- термогенез грелок (метаболизм + печка) ---
P_heat = 0.0
for idx in heater_idx:
Tai = T_cur[idx, I_AI]
Pmet = max((29.133 - 0.739*Tai)*1e-3, Pmin)
P_heat += Pmet + (P_heater if hon[idx]==1 else 0.0)
# --- потери во вне от каждого сота: внешние входы + граничные стенки ---
P_out = [0.0, 0.0]
for i in range(N_CELLS):
fr = cells[i][0]
si = street_idx[i]
Tai = T_cur[i, I_AI]
if si < 0: # внешняя сторона: вход -> T_ambient
if States[i]==CELL_BEE:
P_out[fr] += hA_*(Tai-T_ambient) + (T_cur[i,I_AB]-T_ambient)/R_ab_
else:
P_out[fr] += hA*(Tai-T_ambient)
for f in range(9): # граничные стенки -> T_ambient
if Topology[i,f]==-1:
P_out[fr] += (T_cur[i,4+f]-T_ambient)/R_wall_self[f]
# --- потери через периметр улочки (граничные узлы -> T_ambient) ---
P_perim = 0.0
for si in range(N_street):
if not is_boundary[si]:
Ts = T_s[si]
for j in range(6):
nb = street_neighbors[si,j]
if nb>=0 and is_boundary[nb]:
P_perim += G_lat*(Ts - T_ambient)
# --- накопление энергии в телах грелок (dU_bee/dt) ---
dU_bee = 0.0
for idx in heater_idx:
for j in range(3): # I_HD, I_TH, I_AB
dU_bee += C_grid[idx, j] * (T_cur[idx, j] - T_prev[idx, j]) / dt_window
P_heat_net = P_heat - dU_bee # в наружу идёт только это
# --- термогенез обсидчиков ---
P_surf = 0.0
for k in range(N_surf):
Tsense = T_cur[surf_head_idx[k], I_AI]
Pmet = max((29.133 - 0.739 * Tsense) * 1e-3, Pmin)
P_surf += Pmet + (P_heater if sho[k] == 1 else 0.0)
# --- накопление энергии в телах обсидчиков ---
dU_surf = 0.0
for k in range(N_surf):
for j in range(3): # голова, торакс, брюшко
dU_surf += C_grid[surf_head_idx[k], j] * (T_cur[surf_head_idx[k], j] - T_prev[surf_head_idx[k], j]) / dt_window
# --- накопление энергии в расплоде ---
dU_brood = 0.0
for idx in brood_idx:
dU_brood += C_brood_arr[idx] * (T_cur[idx, I_TH] - T_prev[idx, I_TH]) / dt_window
P_surf = Q_dorsal_sum
P_surf_net = P_surf - dU_surf
P_brood_net = P_brood - dU_brood
return P_brood, P_heat, P_heat_net, P_surf, P_surf_net, P_brood_net, P_out[0], P_out[1], P_perim
def total_energy(T_cur, Ts):
"""Полная тепловая энергия: все узлы ячеек + воздух улочки."""
return float(np.sum(C_grid * T_cur)) + C_node * float(np.sum(Ts))
def run_sim(time_total, log=False, use_surf=True):
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()
heater_on_time = np.zeros(N_CELLS, np.float64)
heater_on_count = np.zeros(N_CELLS, np.int32)
hon_prev = hon.copy()
Ts=np.full(N_street,T_ambient); Ts_new=np.full(N_street,T_ambient)
# --- обсидчики ---
Tsb = 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()
surf_on_time = np.zeros(max(N_surf, 1), np.float64)
surf_on_count = np.zeros(max(N_surf, 1), np.int32)
sho_prev = sho.copy() if N_surf else np.zeros(1, np.int32)
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)
print(f" DEBUG run_sim: use_surf={use_surf}, Cover_frac.sum()={Cf.sum():.2f}")
T_cur, T_next = TA, TB
log_t, log_Tbr0, log_Tstr_c, log_Tstr_m, log_hon = [], [], [], [], []
log_Tht, log_TthS, log_surf_on = [], [], []
log_Psrc=[]; log_Ploss=[]; log_U=[]
T_prev = T_cur.copy()
for step in range(steps):
Qe=do_cluster_step(T_cur,T_next,Topology,Neighbor_Faces,States,C_grid,P_extra,
P_brood_arr,C_brood_arr,hon,is_capped,
Cf,Cb,Cp,Tsb,sho,Ts,street_idx,dt_sim,T_ambient)
T_cur,T_next = T_next,T_cur
# --- учёт включений грелок (Python, без изменения ядра) ---
if heater_idx:
hon_cur = hon[heater_idx]
heater_on_time[heater_idx] += np.where(hon_cur == 1, dt_sim, 0.0)
heater_on_count[heater_idx] += ((hon_cur == 1) & (hon_prev[heater_idx] == 0)).astype(np.int32)
hon_prev[heater_idx] = hon_cur
# --- учёт включений обсидчиков (Python, без изменения ядра) ---
if N_surf > 0 and sho_prev.size > 0:
sho_cur = sho[:N_surf]
surf_on_time[:N_surf] += np.where(sho_cur == 1, dt_sim, 0.0)
surf_on_count[:N_surf] += ((sho_cur == 1) & (sho_prev[:N_surf] == 0)).astype(np.int32)
sho_prev[:N_surf] = sho_cur
if use_surf and N_surf>0:
Q_dorsal_node,*_ = do_surface_bees_step(Tsb,T_cur,Ts,dt_sim,
surf_head_idx,surf_thorax_idx,surf_ab_idx,
surf_node_head,surf_node_thorax,surf_node_ab,
N_surf,States,is_capped,sho,sht)
else:
Q_dorsal_node = np.zeros(N_street)
do_street_step(Ts,Ts_new,Qe,Q_dorsal_node,street_cell0,street_cell1,street_neighbors,
is_boundary,dt_sim,T_ambient,C_node,G_lat)
Ts,Ts_new = Ts_new,Ts
# --- сброс счётчиков в середине прогона: статистика только за 2-й час ---
if step == int(time_total / 2 / dt_sim):
heater_on_time[:] = 0.0
heater_on_count[:] = 0
surf_on_time[:] = 0.0
surf_on_count[:] = 0
# hon_prev и sho_prev НЕ сбрасываем, чтобы не создать ложный переход 0→1
if log and step%100==0:
log_t.append(step*dt_sim)
log_Tbr0.append(T_cur[brood_idx,I_TH].mean() if brood_idx else T_ambient)
ci_c,ri_c=n_cols//2,n_rows//2
log_Tstr_c.append(Ts[ci_c*n_rows+ri_c])
log_Tstr_m.append(Ts[~is_boundary].mean())
log_hon.append(hon[heater_idx].mean() if heater_idx else 0.0)
log_Tht.append(T_cur[heater_idx, I_TH].mean() if heater_idx else T_ambient)
log_TthS.append(Tsb[:,1].mean() if N_surf else T_ambient) # <-- 0.0="" 0:="" 6000="=" and="" b="" b_net="" cur="" dt_sim="" else="" h="" h_net="" if="" log="" log_ploss.append="" log_psrc.append="" log_surf_on.append="" log_u.append="" mean="" n_surf="" o0="" o1="" p="" parts="[f" pb="" s="" s_net="" sho="" step="" surf="" t="{t_sec:.0f}" t_prev="" t_sec="step" total_energy=""> 0:
ns_on = int(sho[:N_surf].sum())
Tsense_arr = T_cur[surf_head_idx, I_AI]
parts.append(f"обсидчиков ВКЛ={ns_on}/{N_surf}, "
f"Tsense min={Tsense_arr.min():.2f}, max={Tsense_arr.max():.2f}")
if len(heater_idx) > 0:
n_on = int(hon[heater_idx].sum())
parts.append(f"грелок ВКЛ={n_on}/{len(heater_idx)}, "
f"Tth={T_cur[heater_idx, I_TH].mean():.1f}, "
f"Tai={T_cur[heater_idx, I_AI].mean():.1f}, "
f"Thd={T_cur[heater_idx, I_HD].mean():.1f}")
print(", ".join(parts))
# извлечение полей
T_air_0h=np.full((n_cols,n_rows),T_ambient); T_br_0h=np.full((n_cols,n_rows),np.nan)
T_air_1Y=np.full((n_cols,n_rows),T_ambient)
for i,cell in enumerate(cells):
frame,side,col,row=cell; ci,ri=col-cmin,row-rmin
if frame==0 and side=='^':
T_air_0h[ci,ri]=T_cur[i,I_AI]
if States[i]==CELL_BROOD: T_br_0h[ci,ri]=T_cur[i,I_TH]
elif frame==1 and side=='Y':
T_air_1Y[ci,ri]=T_cur[i,I_AI]
T_br_1Y = np.full((cols, rows), np.nan)
for i, cell in enumerate(cells):
fr, side, col, row = cell
if fr == 1 and side == 'Y':
ci, ri = col - cmin, row - rmin
if States[i] == CELL_BROOD:
T_br_1Y[ci, ri] = T_cur[i, I_TH]
Tstr=np.full((n_cols,n_rows),T_ambient)
for ci in range(n_cols):
for ri in range(n_rows): Tstr[ci,ri]=Ts[ci*n_rows+ri]
# --- БАЛАНС и dU/dt за последние 600 с ---
if log and len(log_U) > 60:
N = min(len(log_U), len(log_Psrc), len(log_Ploss))
Psrc_avg = np.mean(log_Psrc[-N:])
Ploss_avg = np.mean(log_Ploss[-N:])
U_start = log_U[-N]
U_end = log_U[-1]
dt_win = (N - 1) * 100 * dt_sim
dU_dt = (U_end - U_start) / dt_win
print(f"\n БАЛАНС (последние {dt_win:.0f} с):")
print(f" Источники: {Psrc_avg*1000:>8.1f} мВт")
print(f" Потери: {Ploss_avg*1000:>8.1f} мВт")
# print(f" Невязка потоков: {(Psrc_avg-Ploss_avg)*1000:>8.2f} мВт")
print(f" dU/dt (по энергии):{dU_dt*1000:>8.2f} мВт")
print(f" Длина логов: log_U={len(log_U)}, log_Psrc={len(log_Psrc)}, log_Ploss={len(log_Ploss)}")
if N_surf > 0:
print(f" Tsense обсидчиков: min={T_cur[surf_head_idx, I_AI].min():.2f}, "
f"max={T_cur[surf_head_idx, I_AI].max():.2f}")
Qe_final = Qe.copy()
return dict(T_cur=T_cur, T_prev=T_prev, Ts=Ts, hon=hon, T_air_0h=T_air_0h, T_br_0h=T_br_0h,
T_air_1Y=T_air_1Y, Tstr=Tstr, T_br_1Y=T_br_1Y, Qe_final=Qe_final,
log_t=log_t, log_Tbr0=log_Tbr0, log_Tstr_c=log_Tstr_c, heater_on_time=heater_on_time,
heater_on_count=heater_on_count, surf_on_time=surf_on_time, surf_on_count=surf_on_count,
log_Tstr_m=log_Tstr_m, log_hon=log_hon, log_Tht=log_Tht, log_TthS=log_TthS,
T_surf_bees=T_surf_bees, log_surf_on=log_surf_on, surf_heater_timer=surf_heater_timer, sho=sho)
# компиляция
dt_window = 100 * dt_sim # 10 с между вызовами
TA0=np.full((N_CELLS,13),T_ambient); TB0=np.full((N_CELLS,13),T_ambient)
hon0=heater_on.copy(); Ts0=np.full(N_street,T_ambient)
Tsb0=T_surf_bees.copy(); sho0=surf_heater_on.copy(); sht0=surf_heater_timer.copy()
Qe0 = do_cluster_step(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,Tsb0,sho0,
Ts0,street_idx,dt_sim,T_ambient)
if N_surf > 0:
Q_dorsal_node,Qst0,Qce0,Qev0,dUk0,mr0 = do_surface_bees_step(
Tsb0,TA0,Ts0,dt_sim,
surf_head_idx,surf_thorax_idx,surf_ab_idx,
surf_node_head,surf_node_thorax,surf_node_ab,
N_surf,States,is_capped,sho0,sht0)
else:
Q_dorsal_node = np.zeros(N_street)
Ts0n=np.full(N_street,T_ambient)
do_street_step(Ts0,Ts0n,Qe0,Q_dorsal_node,street_cell0,street_cell1,street_neighbors,
is_boundary,dt_sim,T_ambient,C_node,G_lat)
print(f">>> ПРОВЕРКА: грелок={len(heater_idx)} (ENABLE_HEATERS={ENABLE_HEATERS}), "
f"обсидчиков={N_surf} (ENABLE_OBSIDCHIKI={ENABLE_OBSIDCHIKI})")
if ENABLE_OBSIDCHIKI and N_surf == 0:
print(" ВНИМАНИЕ: обсидчики разрешены, но не сгенерировались!")
if ENABLE_HEATERS and len(heater_idx) == 0:
print(" ВНИМАНИЕ: грелки разрешены, но не сгенерировались!")
# --- Шаг 1а: верификация (без расплода) ---
print(f"\n=== Шаг 1а: верификация (без расплода, 60 с) ===")
States_save=States.copy(); P_extra_save=P_extra.copy(); P_brood_save=P_brood_arr.copy()
C_brood_save=C_brood_arr.copy(); heater_on_save=heater_on.copy(); is_capped_save=is_capped.copy()
States[:]=CELL_EMPTY; P_extra[:]=0; P_brood_arr[:]=0; C_brood_arr[:]=0; heater_on[:]=0; is_capped[:]=0
t0t=time.time(); res_verif=run_sim(60.0, log=False, use_surf=False)
dev=np.max(np.abs(res_verif['Ts']-T_ambient))
print(f" за {time.time()-t0t:.2f} с. max|T_street - T_amb| = {dev:.2e} °C {'✓' if dev<1e-10 ---="" 1="" and="" brood_save="" c_brood_arr="" cover_frac.sum="" debug:="" elapsed:.2f="" elapsed="time.time()-t0t" else="" enable_obsidchiki="" extra_save="" f="" heater_on="" heater_on_save="" if="" is_capped="" is_capped_save="" log="True," n="==" n_surf="" over_frac.sum="" p_brood_arr="" p_extra="" print="" res="run_sim(time_total," res_base="run_sim(time_total," states="" t0t="time.time();" tates_save="" time_total:.0f="" use_surf="False)" verif=""> 0:
print("\n=== Прогон с обсидчиками ===")
res_surf = run_sim(time_total, log=True, use_surf=True)
has_surf = True
else:
print("\n=== Обсидчики выключены — второй прогон пропущен ===")
res_surf = res_base
has_surf = False
# ==========================================================
# ВИЗУАЛИЗАЦИЯ РЕЗУЛЬТАТОВ (для блога, в столбик)
# ==========================================================
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')
Tstr=res_base['Tstr']; T_air_0h=res_base['T_air_0h']; T_air_1Y=res_base['T_air_1Y']; T_br0=res_base['T_br_0h']; T_br_1Y=res_base['T_br_1Y']
# общая шкала для температурных полей
allv=np.concatenate([Tstr.ravel(),T_air_0h.ravel(),T_air_1Y.ravel()])
if np.any(~np.isnan(0)): allv=np.concatenate([allv,T_br0[~np.isnan(T_br0)]])
gmin,gmax=float(allv.min()),float(allv.max())
ext=extent(); fw=10; fh=fw*(ext[3]-ext[2])/(ext[1]-ext[0])
# --- РИС. Поля температур: улочка + обе стороны (в столбик) ---
fig,(ax1,ax2,ax3)=plt.subplots(3,1,figsize=(fw/2,fh*2))
draw_hex(ax2,fig,Tstr,'inferno',gmin,gmax,title=f'T улочки (d={d_street*1e3:.0f} мм)',cb='°C')
draw_hex(ax1,fig,T_air_0h,'inferno',gmin,gmax,title="T воздуха, сот 0 сторона '^'",cb='°C',brood_field=T_br0)
draw_hex(ax3,fig,T_air_1Y,'inferno',gmin,gmax,title="T воздуха, сот 1 сторона 'Y'",cb='°C',brood_field=res_base['T_br_1Y'])
plt.tight_layout(); plt.savefig('fig_street_fields.png',dpi=150,bbox_inches='tight'); plt.show()
# --- РИС. Расплод: T_brood + dT от T_amb ---
fig,(ax1,ax2)=plt.subplots(2,1,figsize=(fw/2,fh*1.8))
draw_hex(ax1,fig,T_br0,'magma',np.nanmin(T_br0),np.nanmax(T_br0),title='T расплода',cb='°C')
dT_str=Tstr-T_ambient; dmax=max(np.nanmax(dT_str),1e-3)
draw_hex(ax2,fig,dT_str,'YlOrRd',0,dmax,title='Нагрев улочки (dT от T_amb)',cb='°C')
plt.tight_layout(); plt.savefig('fig_street_brood_dT.png',dpi=150,bbox_inches='tight'); plt.show()
# --- РИС. Динамика ---
log=res_base
fig,(a1,a2,a3)=plt.subplots(3,1,figsize=(9,9),sharex=True)
a1.plot(log['log_t'],log['log_Tbr0'],'g-',lw=2,label='T расплода (средняя)')
a1.axhline(T_ambient,color='gray',ls=':',lw=1)
a1.axhline(34.5,color='blue',ls='--',lw=1.2,label='минимум 34.5°C')
a1.axhline(35.5,color='red',ls='--',lw=1.2,label='цель 35.5°C')
a1.set_ylabel('°C',fontsize=FS); a1.set_title('T расплода',fontsize=FS); a1.legend(fontsize=FS-1); a1.grid(ls=':')
a2.plot(log['log_t'],log['log_Tstr_c'],'r-',lw=2,label='T улочки (центр)')
a2.plot(log['log_t'],log['log_Tstr_m'],'b-',lw=2,label='T улочки (среднее)')
a2.axhline(T_ambient,color='gray',ls=':',lw=1)
a2.set_ylabel('°C',fontsize=FS); a2.set_title('T улочки',fontsize=FS); a2.legend(fontsize=FS-1); a2.grid(ls=':')
a3.plot(log['log_t'],np.array(log['log_hon'])*100,'r-',lw=2,label='доля грелок в режиме "печка"')
a3.set_ylabel('%',fontsize=FS); a3.set_xlabel('время, с',fontsize=FS)
a3.set_title('Грелки',fontsize=FS); a3.legend(fontsize=FS-1); a3.grid(ls=':')
plt.tight_layout(); plt.savefig('fig_street_dynamics.png',dpi=150,bbox_inches='tight'); plt.show()
# ==========================================================
# СВОДНАЯ ТАБЛИЦА (числа для текста раздела 3)
# ==========================================================
scenarios = [('Грелки (без обсидчиков)', res_base)]
if has_surf:
scenarios.append(('Грелки + обсидчики', res_surf))
def heater_stats(name, on_time, on_count, idx_arr, total_time):
"""Средняя длительность включения и число включений на пчелу."""
if idx_arr is None or len(idx_arr) == 0:
print(f" {name:<32 100="" avg_dur="t_sum" duty="t_sum" if="" n_bees="" n_sum="" return="" t_sum="float(on_time[idx_arr].sum())" total_time=""> 0 else 0.0
avg_cnt = n_sum / n_bees
print(f" {name:<32 duty:="">5.1f}% ср.длит. {avg_dur:>7.1f} с включ./пчелу {avg_cnt:>5.1f}")
print(f"\n{'='*78}")
print(f" СТАТИСТИКА РЕЖИМА ПЕЧКИ (t={time_total:.0f} с)")
print(f"{'='*78}")
heater_stats("Грелки (без обсидчиков)",
res_base['heater_on_time'], res_base['heater_on_count'],
heater_idx, time_total)
heater_stats("Грелки (с обсидчиками)",
res_surf['heater_on_time'], res_surf['heater_on_count'],
heater_idx, time_total)
if ENABLE_OBSIDCHIKI and N_surf > 0:
heater_stats("Обсидчики",
res_surf['surf_on_time'], res_surf['surf_on_count'],
np.arange(N_surf), time_total)
print(f"{'='*78}")
print(f"\n{'='*72}")
print(f" СВОДНАЯ ТАБЛИЦА: T расплода ({len(brood_idx)} ячеек), T_amb={T_ambient}°C")
print(f"{'='*72}")
print(f" {'Сценарий':<30>7} {'мин':>7} {'макс':>7} {'ΔT':>7}")
print(f" {'-'*62}")
for name, res in scenarios:
Tb = res['T_cur'][brood_idx, I_TH]
print(f" {name:<30 b.mean="">7.2f} {Tb.min():>7.2f} {Tb.max():>7.2f} {Tb.mean()-T_ambient:>+7.2f}")
print(f"{'='*72}")
Q_dorsal_sum = float(np.sum(Q_dorsal_node[~is_boundary]))
P_brood, P_heat, P_heat_net, P_surf, P_surf_net, P_brood_net, P_out0, P_out1, P_perim = calc_power_balance(res_surf['T_cur'],
res_surf['T_prev'], res_surf['Ts'], res_surf['hon'], res_surf['sho'], dt_window, Q_dorsal_sum)
P_src_net = P_brood_net + P_heat_net + P_surf_net
P_loss = P_out0 + P_out1 + P_perim
print(f" Термогенез расплода (чистый): {P_brood_net*1000:>8.1f} мВт")
print(f" Термогенез грелок (чистый): {P_heat_net*1000:>8.1f} мВт")
print(f" Термогенез обсидчиков (чистый): {P_surf_net*1000:>8.1f} мВт")
print(f" ВСЕГО источников (чистых): {P_src_net*1000:>8.1f} мВт")
print(f" ВСЕГО потерь: {P_loss*1000:>8.1f} мВт")
print(f" Невязка: {(P_src_net - P_loss)*1000:>8.2f} мВт")
Q_entr_sum = 0.0
Qe_final = res_surf['Qe_final']
for si in range(N_street):
if not is_boundary[si]:
Q_entr_sum += Qe_final[street_cell0[si]] + Qe_final[street_cell1[si]]
print(f" ΣQ_entrance={Q_entr_sum*1000:.1f} мВт, ΣQ_dorsal={Q_dorsal_sum*1000:.1f} мВт, P_perim={P_perim*1000:.1f} мВт")
if has_surf:
log_b = res_base; log_s = res_surf
fig,(a1,a2,a3)=plt.subplots(3,1,figsize=(9,9),sharex=True)
# --- T расплода ---
a1.plot(log_b['log_t'],log_b['log_Tbr0'],'g--',lw=1.5,alpha=0.7,label='без обсидчиков')
a1.plot(log_s['log_t'],log_s['log_Tbr0'],'g-', lw=2, label='с обсидчиками')
a1.axhline(T_ambient,color='gray',ls=':',lw=1)
a1.axhline(34.5,color='blue',ls='--',lw=1.2,label='минимум 34.5°C')
a1.axhline(35.5,color='red', ls='--',lw=1.2,label='цель 35.5°C')
a1.set_ylabel('°C',fontsize=FS); a1.set_title('T расплода (средняя)',fontsize=FS)
a1.legend(fontsize=FS-1); a1.grid(ls=':')
# --- T улочки ---
a2.plot(log_b['log_t'],log_b['log_Tstr_m'],'b--',lw=1.5,alpha=0.7,label='без обсидчиков')
a2.plot(log_s['log_t'],log_s['log_Tstr_m'],'b-', lw=2, label='с обсидчиками')
a2.axhline(T_ambient,color='gray',ls=':',lw=1)
a2.set_ylabel('°C',fontsize=FS); a2.set_title('T улочки (средняя)',fontsize=FS)
a2.legend(fontsize=FS-1); a2.grid(ls=':')
# --- доля в печке ---
a3.plot(log_b['log_t'],np.array(log_b['log_hon'])*100,'r--',lw=1.5,alpha=0.7,label='грелки (без обсидчиков)')
a3.plot(log_s['log_t'],np.array(log_s['log_hon'])*100,'r-', lw=2, label='грелки (с обсидчиками)')
a3.plot(log_s['log_t'],np.array(log_s['log_surf_on'])*100,'m-',lw=2,label='обсидчики')
a3.set_ylabel('%',fontsize=FS); a3.set_xlabel('время, с',fontsize=FS)
a3.set_title('Доля пчёл в режиме печки',fontsize=FS)
a3.legend(fontsize=FS-1); a3.grid(ls=':')
plt.tight_layout()
plt.savefig('fig_street_dynamics_cmp.png',dpi=150,bbox_inches='tight'); plt.show()
fig, ax = plt.subplots(figsize=(9, 5))
# Грелки (без обсидчиков)
ax.plot(log_b['log_t'], log_b['log_Tht'], 'r--', lw=1.5, alpha=0.7,
label='грелки (без обсидчиков)')
# Грелки (с обсидчиками)
ax.plot(log_s['log_t'], log_s['log_Tht'], 'r-', lw=2,
label='грелки (с обсидчиками)')
# Обсидчики
ax.plot(log_s['log_t'], log_s['log_TthS'], 'm-', lw=2,
label='обсидчики (средняя)')
ax.axhline(T_ambient, color='gray', ls=':', lw=1)
ax.set_ylabel('T торакса, °C', fontsize=FS)
ax.set_xlabel('время, с', fontsize=FS)
ax.set_title('Средняя температура торакса', fontsize=FS)
ax.legend(fontsize=FS-1)
ax.grid(ls=':')
plt.tight_layout()
plt.savefig('fig_thorax_temp_cmp.png', dpi=150, bbox_inches='tight')
plt.show()
Комментарии
Отправить комментарий