В самом начале работы над моделью я, конечно же, искал готовые решения, и видел, что есть статья J.A.C. Humphrey, E.S. Dykes Thermal energy conduction in a honey bee comb due to cell-heating bees [30], где решалась задача нестационарного теплопереноса для кластера 20х20 ячеек с разным содержимым ячеек и пчелами-грелками. Однако полного доступа к статье у меня не было, и пришлось делать модель самому, чем я и занимаюсь уже почти год. Как я уже писал в заметке Статус ячейки с расплодом 2, мне очень повезло, что блог посетил В.Г. Маршаков и предложил помощь с первоисточниками. В его огромном архиве нашлась копия этой статьи, и он мне ее прислал. Еще раз огромное спасибо, Владимир Георгиевич!
После того, как сделаны модели всех статусов ячейки и модель теплового взаимодействия ячеек в соте можно попытаться промоделировать участок сота с разными типами ячеек. Тем более, что теперь есть с чем сравнить. Поэтому возьмем такой же кластер 20х20 ячеек и заполним его аналогичным образом: в центре расплод и пять пустых ячеек (для помещения туда в дальнейшем пчел), далее - кольцо перги, а за ним пустые ячейки, если они ниже середины по высоте, или мед - если выше:
И рассмотрим аналогичные три случая.
Случай 1
Там описан так: В этом случае проводится расчёт для установления стационарного температурного распределения в соте для распределения ячеек, показанного на рис. 1 , при условии, что в зоне выводка нет пчёл, нагревающих клетки, но куколки генерируют тепло с относительно низкой скоростью — 103 W/m³ (1W/kg). Стационарное состояние — это то, при котором максимальное изменение температуры в любой точке гребня составляет меньше 0,011°C с одного шага времени на другой. Продолжительность шага не указана.
Заданный уровень термогенеза в 1 Вт/кг соответствует примерно 0.1 мВт на куколку. Зададим возраст расплода, для примера, 390-400 часов, и, согласно Статус ячейки с расплодом 2, на каждую куколку в этом возрасте приходится более 0.45 мВт. Тогда общая мощность термогенеза расплода (2 стороны по 40 ячеек) составит более 38 мВт. Почему две стороны, ведь в статье - только одна? Дело в том, что в статье предполагается, что вторая сторона - симметрична первой, средостение - адиабатическое, поэтому между сторонами практически нет теплообмена. В моей модели есть возможность задать вторую сторону симметричной.
Распределение температуры по ячейкам можно увидеть на рисунке:
Так же, как и авторы статьи, наблюдаем асимметрию распределения температур относительно центральной горизонтальной оси - вниз тепла идет больше и оно распространяется дальше, чем вверх, т.к. вверху мед, а внизу пустые ячейки. Как видим, несмотря на то, что заданная нами мощность термогенеза расплода в 4.5 раза больше, температура в центральной части поднялась всего до 34.34°С, тогда как в статье - до 34.5°С. Дело, видимо, в том, что авторы посчитали, что в улочки потерь нет, а они есть:
В каждую сторону сота уходит примерно 4.8 мВт, т.е. суммарно в улочки уходит 9.6 мВт из 39 - 25%.
Случай 2
Используя стационарное температурное распределение, полученное из Случая 1, в качестве начального граничного условия, выполняется расчёт с одной пчелой, нагревающей ячейку, расположенной в центре соты, генерирующей тепло со скоростью 2,9105 W/m³ (290W/kg). Пчела генерирует тепло в клетке в течение 10 минут, после чего рассматриваются два варианта: (i) пчела остаётся в клетке ещё 10 минут, но при этом вырабатывает тепло с меньшей скоростью 2,9104 W/m³ (29W/kg) благодаря нормальным метаболическим процессам; (ii) пчела покидает камеру. В обоих случаях 10-минутная фаза нагрева сменяется 10-минутной фазой охлаждения общей продолжительностью 20 минут.
В своем эксперименте сохраним ту же логику, но мощность генерируемая пчелой в режиме печки составляет 45 мВт плюс мощность метаболизма в зависимости от температуры в ячейке: $$P_{met} = 29.133 - 0.739 * T$$ (мВт), и эта же мощность термогенеза (но уже без 45 мВт печки) остается на вторые 10 минут (i), а во втором варианте (ii) пчела также просто удаляется и мощность термогенеза пустой ячейки становится равной нулю.
Сначала посмотрим режим разогрева - первые 10 минут. Поскольку установившийся режим был достигнут через 1 час, т.е. к 3600 секундам, то +10 минут заканчиваются к 4200 секундам:
Видим, что к концу десятой минуты нагрева качественно тепловая карта похожа на цитируемую - пчела греет не только соседние ячейки, но и ячейки второго, третьего и даже четвертого от нее ряда. Температура воздуха в центральной ячейке, в которую вошла пчела-печка пришла к значению 37.21°С, что несколько ниже, чем в статье - 37.4°С, но довольно близко. Но вот по динамике есть качественное отличие. Дело в том, что авторы говорят о температуре в ячейке, как о температуре в центре входа в ячейку, а там воздух, поэтому и я вывожу на график температуру воздуха в ячейке. Ввиду малой теплоемкости воздуха в ячейке (тем более с пчелой) я полагаю ее одинаковой по всему объему ячейки. При входе горячей пчелы в ячейку воздух быстро нагревается, а потом отдает тепло стенкам ячейки и достигает равновесной температуры. В статье температура в ячейке растет плавно, без первоначального скачка. Это, считаю, не важно для итогового результата.
Второе существенное отличие - в нашей модели видно, сколько тепла уходит в улочки. Как было показано ранее, из генерируемых расплодом 39 мВт примерно 25% уходило в улочки. При включении одной пчелы-печки из 84 мВт (39+45) в улочки суммарно уходит 18 мВт - около 21%, но видно, что за 10 минут установившееся состояние не было достигнуто - при рассмотрении "случая 1" было видно, что для его достижения нужно около часа.
Третье отличие заключается в температуре соседних ячеек. В статье в соседних с пчелой ячейках температура достигла 37°С, а здесь - 35.27; в следующем ряду - 35.5, а здесь - 34.78. Дело, видимо, в том, что здесь в модели ячейки использованы коричневые соты с большим тепловым сопротивлением, чем в статье.
При "выключении" грелки на 10 минут (i):
К окончанию десятой минуты (двадцатой от входа пчелы в ячейку) температура в центральной ячейке упала до 34.54°С, что выше первоначально установившейся температуры - 34.34°С, но не так значительно, как в статье - там более 36.
При удалении пчелы после 10-минутного прогрева (ii):
Здесь остывание центральной ячейки прошло до температуры 34.46°С (в статье - до 35.4), а в остальном не сильно отличается от охлаждения не с удалением, а "выключением" пчелы - часть тепловой энергии осталась (накопилась) в соте.
Случай 3
Этот случай такой же, как и случай 2, за исключением того, что пять пчёл, нагревающих клетки, занимают пять клеток в области соты в течение первых 10 минут, после чего альтернативы (i) и (ii) выше рассматриваются ещё 10 минут с общим временем 20 минут. Обратите внимание, что скорости выработки энергии, предполагаемые для пассивно и активно нагревающих взрослых рабочих пчёл, находятся в пределах диапазонов, приведённых Сили (1985) и Фернхольца и др. (1989); См. обсуждение выше. При отсутствии сопоставимых данных по куколкам мы предположили, что скорость генерации энергии равна 3,5% от скорости пассивно нагревающейся взрослой пчелы. Главный момент в том, что используемые скорости генерации энергии представляют разумные приближения к природной ситуации и, как показано ниже, приводят к интересным, физически реалистичным результатам.
После 10-минутного разогрева пятью пчелами:
Максимальная температура в центре достигла 38.94°С, в статье - 41°С.
Для обеспечения условия адиабатичности средостения нужно на стороне Х также разместить 5 пчел, тогда условия будут ближе, к рассмотренным в статье:
В этом случае, как видим, максимальная температура в центре достигла 40.41°С.
Для сокращения объема заметки оставим только вариант (ii), т.е. 10 минут нагрева десятью пчелами, с последующим их "выключением" без удаления:
Охлаждение центральной ячейки дошло до 35.39°С, в статье - до 38.7.
Выводы
1. Разработанная модель позволяет проводить эксперименты с различной комбинацией состава сот, включая мед, пергу, разновозрастный расплод, пчел.
2. Сравнение с ранее разработанной моделью, авторы которой указывают на ее непротиворечивость имеющимся экспериментальным данным, подтверждает адекватность предлагаемой здесь модели. Расхождение в некоторых абсолютных значениях объясняются более точным учетом здесь мощности термогенеза расплода, тепловых свойств ячеек, потерь тепла в улочки, мощности пчелы-печки.
🐝 Python-скрипт
""" Программа к заметке Моделирование кластера с разными типами ячеек """
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math
import time
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
import pickle
import os
# ==========================
# 1. КОНСТАНТЫ ТИПОВ ЯЧЕЕК
# ==========================
CELL_EMPTY = 0 # Пустая ячейка
CELL_HONEY_SEALED = 1 # Запечатанный мед
CELL_HONEY_OPEN = 2 # Открытый мед
CELL_POLLEN = 3 # Перга (пыльца)
CELL_BROOD = 4 # Расплод
CELL_BEE = 5 # Пчела
# ==========================
# 2. ФИЗИЧЕСКИЕ ПАРАМЕТРЫ
# ==========================
# Параметры пчелы
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 = 3500
C_hd, C_th, C_ab = m_head * c, m_thorax * c, m_ab * c
C_bee = C_hd + C_th + C_ab
# Параметры воздуха
ρ_air, k_air, c_air = 1.205, 0.025, 1005
T_ambient = 34.0
# Параметры ячейки
variant = 3 # 1 ... 3; вариант 1 соответствует первоначальной гипотезе eps_wall = eps = 0.22
L_w = [0.0001, 0.00014, 0.0002]
L_b = [0.00013, 0.00055, 0.0012]
C_w = [53.324, 103.878, 140.235]
c_w = [3300., 3000., 2700.] # удельная теплоемкость, Дж/(кг·°С)
λ_w = [0.235, 0.1608, 0.055] # теплопроводность Вт/(м·К)
L_wall = L_w[variant - 1]
L_bottom = L_b[variant - 1]
thread = 5.4e-3 # шаг решетки
d_cell = thread - L_w[variant - 1] # расстояние между гранями
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. # площадь ромба
# Параметры воска
ρ_wax, λ_wax, c_wax = 960, 0.25, c_w[variant - 1]
M_lateral = ρ_wax * A_lateral * L_w[variant - 1] / 2
M_bottom = ρ_wax * A_bottom * L_b[variant - 1] / 2
# Параметры излучения
sigma = 5.67e-8
t0 = -273.15
eps = 0.22
eps_w = [0.22, 0.5, 0.9]
eps_wall = eps_w[variant - 1]
# Параметры меда
ρ_honey = 1420 # кг/м³, плотность меда
c_honey = 1380 # Дж/(кг·К), удельная теплоемкость меда
λ_honey = 0.5 # Вт/(м·К), теплопроводность меда
# Параметры перги (пыльцы)
ρ_pollen = 600 # кг/м³, плотность перги
c_pollen = 2000 # Дж/(кг·К), удельная теплоемкость перги
λ_pollen = 0.2 # Вт/(м·К), теплопроводность перги
# параметры нагрева
P_warm_up = 45e-3 # в режиме печки
# P_warm_up = 0.
def metabolic_power(T):
power = (29.133 - 0.739 * T) * 1e-3
return max(power, Pmin)
# ==========================
# 3. ГЕОМЕТРИЯ И ТЕПЛОЁМКОСТИ
# ==========================
a = d_cell / (2 * math.sin(math.pi / 3))
area_hex = (3 * math.sqrt(3) / 2) * a ** 2
A_entrance = area_hex
V_cell = area_hex * h_total
# Базовые теплоемкости
C_air_cell = ρ_air * V_cell * c_air
C_air_bee = ρ_air * (V_cell - V_bee) * c_air
C_wax_lateral = ρ_wax * A_lateral * L_w[variant - 1] * c_wax
C_wax_bottom = ρ_wax * A_bottom * L_b[variant - 1] * c_wax
C_wall = [C_wax_lateral] * 6 + [C_wax_bottom] * 3
print(f'Базовые теплоемкости (Дж/°С):')
print(f' Воздух в пустой ячейке: {C_air_cell:.6f}')
print(f' Пчела: {C_bee:.6f}')
print(f' Воздух в ячейке с пчелой: {C_air_bee:.6f}')
print(f' Воск ячейки: {sum(C_wall):.6f}')
# Теплоемкости содержимого ячеек
C_honey_full = ρ_honey * V_cell * c_honey # Полностью заполненная медом
C_pollen_full = ρ_pollen * V_cell * c_pollen # Полностью заполненная пергой
print(f' Мед в ячейке: {C_honey_full:.6f}')
print(f' Перга в ячейке: {C_pollen_full:.6f}')
# Коэффициенты теплообмена
h_in = k_air / ((d_cell - d_bee) / 2)
h_entrance = 11.0
h_walls = 5.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)
# Тепловые сопротивления от воздуха к стенкам
R_ai_wall = [L_wall / (k_air * A) for A in [A_lateral] * 6 + [A_bottom] * 3]
# Тепловые сопротивления восковых стенок
R_lateral = L_wall / (λ_wax * A_lateral)
R_bottom = L_bottom / (λ_wax * A_bottom)
# ==========================
# 4. МОДЕЛЬ РАСПЛОДА
# ==========================
# Параметры развития расплода (в часах)
t_egg = 0.0 # начало - откладка яйца
t_sealed = 203.0 # запечатывание
t_transition = 227.0 # переход к полиномиальной зависимости
t_emergence = 503.0 # выход пчелы
# Массовые параметры (в мг)
mass_initial = 0.11 # начальная масса яйца
mass_sealed = 140.0 # масса при запечатывании
mass_final = 100.0 # масса при выходе
def brood_mass(t_hours):
"""Масса расплода в мг в зависимости от времени t (часы)"""
if t_hours < t_sealed:
# Экспоненциальный рост от 0.11 мг до 140 мг
k = np.log(mass_sealed / mass_initial) / t_sealed
return mass_initial * np.exp(k * t_hours)
else:
# Линейный спад от 140 мг до 100 мг
slope = (mass_sealed - mass_final) / (t_emergence - t_sealed)
return mass_sealed - slope * (t_hours - t_sealed)
def brood_metabolic_power_microW(t_hours):
"""Тепловыделение расплода в мкВт"""
# Рассчитаем массу
m = brood_mass(t_hours)
if t_hours < t_sealed:
# Режим 1: до запечатывания - пропорционально массе
return m * 8.0 # мкВт
elif t_hours < t_transition:
# Режим 2: линейный спад в течение 24 часов после запечатывания
# Максимальное тепловыделение в момент запечатывания
P_max = brood_mass(t_sealed) * 8.0 # мкВт
# Значение по полиному в точке 227 часов
x = t_transition
P_poly_at_227 = (0.000010 * x ** 2 - 0.005019 * x + 0.877249) * 1000 # мВт → мкВт
# Линейная интерполяция
fraction = (t_hours - t_sealed) / (t_transition - t_sealed)
return P_max + fraction * (P_poly_at_227 - P_max)
else:
# Режим 3: полиномиальная зависимость
x = t_hours
P_mW = 0.000010 * x ** 2 - 0.005019 * x + 0.877249
P_microW = P_mW * 1000
return max(P_microW, 0)
def brood_metabolic_power(t_hours):
"""Тепловыделение расплода в Ваттах"""
return brood_metabolic_power_microW(t_hours) * 1e-6
def brood_heat_capacity(t_hours):
"""Теплоемкость расплода в Дж/°C"""
m_kg = brood_mass(t_hours) * 1e-6 # мг → кг
c_brood = 3500 # Дж/(кг·К)
return m_kg * c_brood
# ==========================
# 5. ФУНКЦИИ ДЛЯ РАЗНЫХ ТИПОВ ЯЧЕЕК
# ==========================
def get_cell_heat_capacity(cell_type, fill_factor=1.0, brood_time=0.0):
"""Возвращает полную теплоемкость ячейки в зависимости от типа"""
# Теплоемкость стенок ячейки (всегда одинаковая)
C_walls_total = sum(C_wall)
if cell_type == CELL_EMPTY: # or cell_type == CELL_BROOD_EMPTY:
# Пустая ячейка: только воздух и стенки
return C_air_cell + C_walls_total
elif cell_type == CELL_HONEY_SEALED or cell_type == CELL_HONEY_OPEN:
# Ячейка с медом: мед + стенки
return C_honey_full + C_walls_total
elif cell_type == CELL_POLLEN:
# Ячейка с пергой: перга + стенки
return C_pollen_full * fill_factor + C_walls_total
elif cell_type == CELL_BROOD:
# Ячейка с расплодом: расплод + воздух + стенки
return brood_heat_capacity(brood_time) + C_air_cell + C_walls_total
elif cell_type == CELL_BEE:
# Ячейка с пчелой: пчела + воздух + стенки
return C_bee + C_air_bee + C_walls_total
return C_air_cell + C_walls_total # по умолчанию
def get_cell_thermal_conductivity(cell_type, fill_factor=1.0):
"""Возвращает эффективную теплопроводность содержимого ячейки"""
if cell_type == CELL_EMPTY or cell_type == CELL_BEE or cell_type == CELL_BROOD:
# Воздух
return k_air
elif cell_type == CELL_HONEY_SEALED or cell_type == CELL_HONEY_OPEN:
# Мед
return λ_honey * fill_factor + k_air * (1 - fill_factor)
elif cell_type == CELL_POLLEN:
# Перга
return λ_pollen * fill_factor + k_air * (1 - fill_factor)
return k_air
def has_convection(cell_type):
"""Определяет, есть ли конвекция через вход ячейки"""
# Конвекция есть только у открытых ячеек
return cell_type in [CELL_EMPTY, CELL_HONEY_OPEN, CELL_BROOD, CELL_BEE]
# ==========================
# 6. ПРАВИЛА СОСЕДСТВА
# ==========================
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):
"""Возвращает (neighbor_addr, neighbor_face_id) для данной грани."""
scomb_id = 0
if face_id < 6:
if side == 'Y':
deltas = ODD_ROW_DELTAS if (row % 2 == 1) else EVEN_ROW_DELTAS
else:
deltas = EVEN_ROW_DELTAS if (row % 2 == 1) else ODD_ROW_DELTAS
dc, dr = deltas[face_id]
n_side = side
n_col = int(col + dc)
n_row = int(row + dr)
else:
k = face_id - 6
if side == 'Y':
if row % 2 == 0:
bottom_list = [(col, row - 1), (col, row), (col - 1, row)]
else:
bottom_list = [(col, row - 1), (col, row), (col + 1, row)]
n_col, n_row = bottom_list[k]
n_side = '^'
else:
if row % 2 == 0:
bottom_list = [(col, row), (col, row + 1), (col + 1, row)]
else:
bottom_list = [(col, row), (col, row + 1), (col - 1, row)]
n_col, n_row = bottom_list[k]
n_side = 'Y'
n_addr = (scomb_id, n_side, n_col, n_row)
if face_id < 6:
if n_side == 'Y':
n_deltas = ODD_ROW_DELTAS if (n_row % 2 == 1) else EVEN_ROW_DELTAS
else:
n_deltas = EVEN_ROW_DELTAS if (n_row % 2 == 1) else ODD_ROW_DELTAS
for rev_id, (ndc, ndr) in enumerate(n_deltas):
if (n_col + ndc, n_row + ndr) == (col, row):
return n_addr, rev_id
else:
if n_side == 'Y':
if n_row % 2 == 0:
n_bottom = [(n_col, n_row - 1), (n_col, n_row), (n_col - 1, n_row)]
else:
n_bottom = [(n_col, n_row - 1), (n_col, n_row), (n_col + 1, n_row)]
else:
if n_row % 2 == 0:
n_bottom = [(n_col, n_row), (n_col, n_row + 1), (n_col + 1, n_row)]
else:
n_bottom = [(n_col, n_row), (n_col, n_row + 1), (n_col - 1, n_row)]
for rev_id, (nc, nr) in enumerate(n_bottom):
if (nc, nr) == (col, row):
return n_addr, 6 + rev_id
return n_addr, 0
# ==========================
# 7. ГЕНЕРАЦИЯ КЛАСТЕРА 20x20
# ==========================
cells = []
cols, rows = 20, 20
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)}
index_to_cell = {i: cell for i, cell in enumerate(cells)}
all_cells_set = set(cells)
# Определяем центральную ячейку и ее индекс
center_cell = (0, 'Y', 0, 0)
center_idx = cell_to_index[center_cell]
# ==========================
# 8. КОНФИГУРАЦИЯ КЛАСТЕРА
# ==========================
cell_states = {} # cell -> (type, fill_factor или brood_time)
brood_times = {} # cell -> время развития расплода (часы)
# Создаем конфигурацию согласно заданию
for cell in cells:
scomb, side, col, row = cell
'''
# По умолчанию все ячейки на стороне '^' пустые
if side == '^':
cell_states[cell] = (CELL_HONEY_SEALED, 0.0)
else:
# На стороне 'Y'
'''
if row < 0:
# Верхняя половина (row < 0) - запечатанный мед
cell_states[cell] = (CELL_HONEY_SEALED, 1.0)
else:
# Нижняя половина - пустые ячейки
cell_states[cell] = (CELL_EMPTY, 0.0)
# Кольцо перги (расстояние от центра меньше 6)
if np.sqrt(col ** 2 + row ** 2) < 6:
cell_states[cell] = (CELL_POLLEN, 1.0)
# Центральная часть с расплодом (расстояние от центра меньше 4)
if np.sqrt(col ** 2 + row ** 2) < 4:
cell_states[cell] = (CELL_BROOD, np.random.uniform(390, 400))
brood_times[cell] = cell_states[cell][1]
# Пять пустых ячеек в центре
if cell in [(0, 'Y', 0, 0), (0, 'Y', -2, -1), (0, 'Y', -2, 1), (0, 'Y', 1, -1), (0, 'Y', 1, 1),
(0, '^', 0, 0), (0, '^', -2, -1), (0, '^', -2, 1), (0, '^', 1, -1), (0, '^', 1, 1)]:
cell_states[cell] = (CELL_EMPTY, 0.0)
if cell in brood_times:
del brood_times[cell]
if cell in [(0, 'Y', 0, 0), (0, '^', 0, 0)]:
cell_states[cell] = (CELL_BEE, 1.0)
# Добавляем пчел в некоторые пустые ячейки в центральной области
bee_cells = [(0, 'Y', 0, 0), (0, '^', 0, 0)]
# ==========================
# 9. КЭШИРОВАНИЕ СВЯЗЕЙ ГРАНЕЙ
# ==========================
face_connections = []
for i, cell in enumerate(cells):
connections_i = []
for face_id in range(9):
n_addr, n_face = get_neighbor_for_face(cell[1], cell[2], cell[3], face_id)
if n_addr in all_cells_set:
j = cell_to_index[n_addr]
connections_i.append((j, n_face))
else:
connections_i.append(None)
face_connections.append(connections_i)
# ==========================
# 10. МОДЕЛЬ ПЧЕЛЫ С 9 ГРАНЯМИ
# ==========================
se = sigma * eps
se_th, se_hd, se_ab_ = se * A_th, se * A_hd, se * (A_ab - s_r)
A_wall_total = 6 * A_lateral + 3 * A_bottom
hA_ = h_entrance * (A_entrance - s_r)
hA = h_entrance * A_entrance
def bee_model_9walls(T_prev, neighbor_T_list, T_sensed, dt):
Thd_prev, Tth_prev, Tab_prev, Tai_prev = T_prev[0:4]
T_wall_prev = T_prev[4:13]
l = (1.3493 * np.exp(-0.058 * T_sensed)) * 1e-3
R_th_hd = l / (k_hd * s_r)
R_th_ab = l / (k_ab * s_r)
def equations(vars):
Thd, Tth, Tab, Tai = vars[0:4]
T_wall = vars[4:13]
q_evap_hd = ke0 * Tth
q_evap_ab = ke0 * Tth
P_rest = metabolic_power(T_sensed)
T_wall_avg = sum(T_wall) / 9
p_s_th = -se_th * ((Tth - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_hd = -se_hd * ((Thd - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_ab = -se_ab_ * ((Tab - t0) ** 4 - (T_wall_avg - t0) ** 4)
p_s_total = p_s_th + p_s_hd + p_s_ab
p_s_to_wall = [(A_lateral if k < 6 else A_bottom) / A_wall_total * p_s_total for k in range(9)]
eq1 = C_th * (Tth - Tth_prev) - dt * (
P_rest + P_warm_up + p_s_th - (Tth - Thd) / R_th_hd - (Tth - Tab) / R_th_ab - (Tth - Tai) / R_th_ai)
eq2 = C_hd * (Thd - Thd_prev) - dt * (
p_s_hd + (Tth - Thd) / R_th_hd - (Thd - Tai) / R_hd_ai - q_evap_hd)
eq3 = C_ab * (Tab - Tab_prev) - dt * (
p_s_ab + (Tth - Tab) / R_th_ab - (Tab - Tai) / R_ab_ai - (Tab - T_ambient) / R_ab_ - q_evap_ab)
Q_ai_to_wall = sum((Tai - T_wall[k]) / R_ai_wall[k] for k in range(9))
eq4 = C_air_bee * (Tai - Tai_prev) - dt * (
(Tth - Tai) / R_th_ai + (Thd - Tai) / R_hd_ai + (Tab - Tai) / R_ab_ai
- Q_ai_to_wall - hA_ * (Tai - T_ambient))
eqs_wall = []
for k in range(9):
R = R_lateral if k < 6 else R_bottom
Q_from_neighbor = (neighbor_T_list[k] - T_wall[k]) / R
Q_from_air = (Tai - T_wall[k]) / R_ai_wall[k]
eq_wall = C_wall[k] * (T_wall[k] - T_wall_prev[k]) - dt * (
Q_from_air + Q_from_neighbor + p_s_to_wall[k])
eqs_wall.append(eq_wall)
return [eq1, eq2, eq3, eq4] + eqs_wall
x0 = T_prev
sol = fsolve(equations, x0, xtol=1e-6, maxfev=300)
if not np.all(np.isfinite(sol)):
return x0
return sol
# ==========================
# 11. МОДЕЛЬ ОБЩЕЙ ЯЧЕЙКИ (ДЛЯ ВСЕХ ТИПОВ КРОМЕ ПЧЕЛЫ)
# ==========================
def general_cell_update(cell_type, T_air_prev, T_wall_prev, neighbor_T_list, fill_factor, brood_time, dt):
"""
Обновление состояния ячейки любого типа (кроме пчелы)
"""
# Определяем теплоемкость ячейки
C_cell = get_cell_heat_capacity(cell_type, fill_factor, brood_time)
# Определяем, есть ли тепловыделение
P_heat = 0.0
if cell_type == CELL_BROOD:
P_heat = brood_metabolic_power(brood_time)
# Поток тепла от содержимого к стенкам
# Для ячеек с содержимым (мед, перга) используем эффективную теплопроводность
λ_eff = get_cell_thermal_conductivity(cell_type, fill_factor)
R_eff = [L_wall / (λ_eff * A) if λ_eff > 0 else 1e10
for A in [A_lateral] * 6 + [A_bottom] * 3]
Q_content_to_wall = 0.0
for k in range(9):
Q_content_to_wall += (T_air_prev - T_wall_prev[k]) / R_eff[k]
# Конвекция через вход (только для открытых ячеек)
Q_conv_entrance = 0.0
if has_convection(cell_type):
Q_conv_entrance = h_entrance * A_entrance * (T_air_prev - T_ambient)
# Обновление температуры воздуха/содержимого
T_air_new = T_air_prev + dt / C_cell * (P_heat - Q_content_to_wall - Q_conv_entrance)
# Обновление температур стенок
T_wall_new = []
for k in range(9):
R = R_lateral if k < 6 else R_bottom
Q_from_neighbor = (neighbor_T_list[k] - T_wall_prev[k]) / R
Q_from_content = (T_air_prev - T_wall_prev[k]) / R_eff[k]
dT = dt / C_wall[k] * (Q_from_content + Q_from_neighbor)
T_wall_new.append(T_wall_prev[k] + dT)
return T_air_new, T_wall_new, P_heat
# ==========================
# 12. НАЧАЛЬНЫЕ УСЛОВИЯ
# ==========================
h_open = 11
R_head_air_open = 1 / (h_open * A_hd)
R_thorax_air_open = 1 / (h_open * A_th)
R_ab_air_open = 1 / (h_open * A_ab)
def calculate_initial_temperatures(T_air):
l = (1.3493 * np.exp(-0.058 * T_air)) * 1e-3
R_th_hd = l / (k_hd * s_r)
R_th_ab = l / (k_ab * s_r)
def equations(vars):
Thd, Tth, Tab = vars
q_evap_hd = ke0 * Tth
q_evap_ab = ke0 * Tth
P_rest = metabolic_power(T_air)
p_s_hd = -eps * sigma * ((Thd - t0) ** 4 - (T_air - t0) ** 4) * A_hd
p_s_th = -eps * sigma * ((Tth - t0) ** 4 - (T_air - t0) ** 4) * A_th
p_s_ab = -eps * sigma * ((Tab - t0) ** 4 - (T_air - t0) ** 4) * A_ab
eq1 = p_s_hd + (Tth - Thd) / R_th_hd - (Thd - T_air) / R_head_air_open - q_evap_hd
eq2 = (P_rest + P_warm_up + p_s_th
- (Tth - Thd) / R_th_hd - (Tth - Tab) / R_th_ab - (Tth - T_air) / R_thorax_air_open)
eq3 = p_s_ab + (Tth - Tab) / R_th_ab - (Tab - T_air) / R_ab_air_open - q_evap_ab
return [eq1, eq2, eq3]
sol = fsolve(equations, [T_air + 10, T_air + 15, T_air + 10], xtol=1e-8)
return sol
# ==========================
# 13. ФУНКЦИИ ДЛЯ СОХРАНЕНИЯ И ЗАГРУЗКИ СОСТОЯНИЯ
# ==========================
def save_simulation_state(filename):
"""
Сохраняет текущее состояние моделирования в файл
"""
global T_air, T_wall, T_bee_states, cell_states, brood_times, bee_cells
global log_time, log_T_center, log_bee_power, log_Q_Y, log_Q_X, log_brood_heat, current_time
state = {
'current_time': current_time,
'T_air': T_air.copy() if T_air is not None else None,
'T_wall': [wall.copy() for wall in T_wall] if T_wall is not None else None,
'T_bee_states': {k: v.copy() for k, v in T_bee_states.items()} if T_bee_states is not None else None,
'cell_states': cell_states.copy() if cell_states is not None else None,
'brood_times': brood_times.copy() if brood_times is not None else None,
'bee_cells': bee_cells.copy() if bee_cells is not None else None,
'log_time': log_time.copy() if log_time is not None else None,
'log_T_center': log_T_center.copy() if log_T_center is not None else None,
'log_bee_power': log_bee_power.copy() if log_bee_power is not None else None,
'log_Q_Y': log_Q_Y.copy() if log_Q_Y is not None else None,
'log_Q_X': log_Q_X.copy() if log_Q_X is not None else None,
'log_brood_heat': log_brood_heat.copy() if log_brood_heat is not None else None,
'variant': variant,
'T_ambient': T_ambient,
'dt': dt
}
with open(filename, 'wb') as f:
pickle.dump(state, f)
print(f"\nСостояние сохранено в файл: {filename}")
print(f"Текущее время моделирования: {current_time:.1f} с")
return True
def load_simulation_state(filename):
"""
Загружает состояние моделирования из файла
"""
global T_air, T_wall, T_bee_states, cell_states, brood_times, bee_cells
global log_time, log_T_center, log_bee_power, log_Q_Y, log_Q_X, log_brood_heat, current_time
if not os.path.exists(filename):
print(f"Файл {filename} не найден. Начинаем с нуля.")
return None
with open(filename, 'rb') as f:
state = pickle.load(f)
print(f"\nСостояние загружено из файла: {filename}")
print(f"Текущее время моделирования: {state['current_time']:.1f} с")
# Обновляем все глобальные переменные
current_time = state['current_time']
T_air = state['T_air']
T_wall = state['T_wall']
T_bee_states = state['T_bee_states']
cell_states = state['cell_states']
brood_times = state['brood_times']
bee_cells = state['bee_cells']
log_time = state['log_time']
log_T_center = state['log_T_center']
log_bee_power = state['log_bee_power']
log_Q_Y = state['log_Q_Y']
log_Q_X = state['log_Q_X']
log_brood_heat = state['log_brood_heat']
return current_time
def add_bees_to_cells(bee_positions):
"""
Добавляет пчел в указанные ячейки
"""
global bee_cells, cell_states, T_bee_states, T_ambient
for cell in bee_positions:
if cell in cell_states and cell not in bee_cells:
# Меняем тип ячейки на пчелу
cell_states[cell] = (CELL_BEE, 1.0)
bee_cells.append(cell)
# Инициализируем состояние пчелы
Thd_init, Tth_init, Tab_init = calculate_initial_temperatures(T_ambient)
T_bee_states[cell] = [Thd_init, Tth_init, Tab_init, T_ambient] + [T_ambient] * 9
print(f"Пчела добавлена в ячейку {cell}")
print(f"Всего пчел: {len(bee_cells)}")
def remove_bees_from_cells(bee_positions):
"""
Удаляет пчел из указанных ячеек
"""
global bee_cells, cell_states, T_bee_states
for cell in bee_positions:
if cell in bee_cells:
# Возвращаем ячейку в исходное состояние (пустую)
cell_states[cell] = (CELL_EMPTY, 0.0)
bee_cells.remove(cell)
# Удаляем состояние пчелы
if cell in T_bee_states:
del T_bee_states[cell]
print(f"Пчела удалена из ячейки {cell}")
print(f"Всего пчел: {len(bee_cells)}")
# ==========================
# 14. ИНИЦИАЛИЗАЦИЯ И МОДЕЛИРОВАНИЕ
# ==========================
# Глобальные переменные, которые будут использоваться в функциях
T_air = None
T_wall = None
T_bee_states = {}
log_time = []
log_T_center = []
log_bee_power = []
log_Q_Y = []
log_Q_X = []
log_brood_heat = []
current_time = 0.0
# Параметры моделирования
dt = 0.01
time_total = 600.0 # 5 минут для тестирования
# Проверяем, есть ли сохраненное состояние
state_file = 'simulation_state.pkl'
if os.path.exists(state_file):
choice = input(f"Найден файл состояния '{state_file}'. Загрузить? (y/n): ")
if choice.lower() == 'y':
current_time = load_simulation_state(state_file)
bee_cells = [(0, 'Y', 0, 0), (0, '^', 0, 0)] # добавляем две пчелы в центре
add_bees_to_cells(bee_cells)
if current_time is not None:
# Корректируем общее время
time_total = current_time + 600.0 # Продолжаем еще 10 минут
print(f"Продолжаем моделирование до {time_total:.1f} с")
else:
current_time = 0.0
print("Начинаем новое моделирование")
else:
current_time = 0.0
# Удаляем старый файл, если начинаем заново
if os.path.exists(state_file):
os.remove(state_file)
print("Начинаем новое моделирование")
else:
current_time = 0.0
print("Начинаем новое моделирование")
# Инициализация температур (только если начинаем с нуля)
if current_time == 0:
T_air = [T_ambient] * len(cells)
T_wall = [[T_ambient] * 9 for _ in range(len(cells))]
# Инициализация состояний пчел
T_bee_states = {}
for cell in cells:
cell_type, param = cell_states[cell]
if cell_type == CELL_BEE:
Thd_init, Tth_init, Tab_init = calculate_initial_temperatures(T_ambient)
T_bee_init = [Thd_init, Tth_init, Tab_init, T_ambient] + [T_ambient] * 9
T_bee_states[cell] = T_bee_init.copy()
# Инициализация логов
log_time = []
log_T_center = []
log_bee_power = []
log_Q_Y = []
log_Q_X = []
log_brood_heat = []
# Начальные значения
log_time.append(0)
# center_idx уже определен в разделе 7
log_T_center.append(T_air[center_idx])
log_Q_Y.append(0)
log_Q_X.append(0)
# Рассчитываем начальное тепловыделение
initial_heat_brood = 0.0
initial_heat_bee = 0.0
for cell in cells:
cell_type, param = cell_states[cell]
if cell_type == CELL_BROOD:
initial_heat_brood += brood_metabolic_power(param)
elif cell_type == CELL_BEE:
initial_heat_bee += metabolic_power(T_ambient) + P_warm_up
log_brood_heat.append(initial_heat_brood * 1000) # мВт
log_bee_power.append(initial_heat_bee * 1000) # мВт
# Основной цикл моделирования
steps_total = int((time_total - current_time) / dt)
start_time_global = time.time()
print(f"\nНачинаем моделирование:")
print(f" Текущее время: {current_time:.1f} с")
print(f" Целевое время: {time_total:.1f} с")
print(f" Осталось шагов: {steps_total}")
print(f" Шаг по времени: {dt} с")
for step in range(steps_total):
t = current_time + step * dt
Q_entrance_Y = 0.0
Q_entrance_X = 0.0
total_heat_generation = 0.0
# Обновление всех ячеек
new_T_air = T_air.copy()
new_T_wall = [T_wall[i].copy() for i in range(len(cells))]
for i in range(len(cells)):
cell = index_to_cell[i]
side = cell[1]
cell_type, param = cell_states[cell]
# Сбор температур за гранями
neighbor_T_list = []
for k in range(9):
conn = face_connections[i][k]
if conn is None:
neighbor_T_list.append(T_ambient)
else:
j, n_face = conn
neighbor_cell = index_to_cell[j]
n_type, n_param = cell_states[neighbor_cell]
# Получаем температуру соседней грани
if n_type == CELL_BEE:
neighbor_T_list.append(T_bee_states[neighbor_cell][4 + n_face])
else:
neighbor_T_list.append(T_wall[j][n_face])
if cell_type == CELL_BEE:
# Ячейка с пчелой
Tai = T_bee_states[cell][3]
Tab = T_bee_states[cell][2]
T_bee_new = bee_model_9walls(T_bee_states[cell], neighbor_T_list, Tai, dt)
T_bee_states[cell] = T_bee_new
new_T_air[i] = T_bee_states[cell][3]
new_T_wall[i] = T_bee_states[cell][4:13]
Q_cell = hA_ * (Tai - T_ambient)
Q_ab = (Tab - T_ambient) / R_ab_
heat_gen = metabolic_power(Tai) + P_warm_up
if side == 'Y':
Q_entrance_Y += Q_cell + Q_ab
else:
Q_entrance_X += Q_cell + Q_ab
total_heat_generation += heat_gen
else:
# Все другие типы ячеек
fill_factor = param if cell_type in [CELL_HONEY_SEALED, CELL_HONEY_OPEN, CELL_POLLEN] else 0.0
brood_time = param if cell_type == CELL_BROOD else 0.0
T_air_new, T_wall_new, P_heat = general_cell_update(
cell_type, T_air[i], T_wall[i], neighbor_T_list,
fill_factor, brood_time, dt
)
new_T_air[i] = T_air_new
new_T_wall[i] = T_wall_new
# Конвекция для открытых ячеек
if has_convection(cell_type):
Q_cell = hA * (T_air[i] - T_ambient)
if side == 'Y':
Q_entrance_Y += Q_cell
else:
Q_entrance_X += Q_cell
total_heat_generation += P_heat
# Обновление времени развития расплода
if cell_type == CELL_BROOD:
brood_times[cell] += dt / 3600.0 # секунды → часы
cell_states[cell] = (CELL_BROOD, brood_times[cell])
T_air = new_T_air
T_wall = new_T_wall
# Логирование каждые 100 шагов
if step % 100 == 0 and step > 0:
log_time.append(t)
log_T_center.append(T_air[center_idx])
# Суммарная мощность пчел
bee_power = 0.0
for cell in bee_cells:
if cell in T_bee_states:
Tai = T_bee_states[cell][3]
bee_power += metabolic_power(Tai) + P_warm_up
log_bee_power.append(bee_power * 1000) # мВт
log_Q_Y.append(Q_entrance_Y * 1000) # мВт
log_Q_X.append(Q_entrance_X * 1000) # мВт
log_brood_heat.append(total_heat_generation * 1000) # мВт
# Автосохранение каждые 5000 шагов (50 секунд)
if step % 5000 == 0 and step > 0:
save_simulation_state('autosave.pkl')
# Прогресс каждые 10000 шагов
if step % 10000 == 0 and step > 0:
elapsed = time.time() - start_time_global
steps_per_sec = step / elapsed if elapsed > 0 else 0
remaining_steps = steps_total - step
remaining_time = remaining_steps / steps_per_sec if steps_per_sec > 0 else 0
print(f"Прогресс: {t:.1f}/{time_total:.1f} с ({100 * t / time_total:.1f}%)")
print(f" Скорость: {steps_per_sec:.0f} шаг/с, осталось: {remaining_time / 60:.1f} мин")
# Проверка на установление
if len(log_T_center) > 2:
recent_changes = np.abs(np.diff(log_T_center[-2:]))
avg_change = np.mean(recent_changes)
if recent_changes < 0.001: # Медленное изменение
print(f" Температура стабилизируется (текущее изменение: {avg_change:.7f} °C/шаг)")
# Обновляем текущее время
current_time = time_total
# Финальное сохранение
save_simulation_state(state_file)
print(f"\nМоделирование завершено за {time.time() - start_time_global:.1f} секунд")
# ==========================
# 15. ВИЗУАЛИЗАЦИЯ РЕЗУЛЬТАТОВ
# ==========================
# 15.1. Карта распределения температур
temp_Y = {}
temp_X = {}
cell_types_Y = {}
cell_types_X = {}
for cell in cells:
scomb, side, col, row = cell
if side == 'Y':
temp_Y[(col, row)] = T_air[cell_to_index[cell]]
cell_types_Y[(col, row)] = cell_states[cell][0]
else:
temp_X[(col, row)] = T_air[cell_to_index[cell]]
cell_types_X[(col, row)] = cell_states[cell][0]
size = 0.8 # Уменьшаем размер для 20x20
width = np.sqrt(3) * size
dy = 1.5 * size
def hex_corner_pointy(center, i):
angle_deg = 60 * i - 30
angle_rad = np.pi / 180 * angle_deg
return (center[0] + size * np.cos(angle_rad), center[1] + size * np.sin(angle_rad))
def primal_to_pixel(col, row):
x = width * (col + 0.5 * (row & 1))
y = -dy * row
return (x, y)
def dual_to_pixel(col, row):
x = width * (col + 0.5 * (1 - (row & 1)))
y = -dy * row - size / 2
return (x, y)
fig, ax = plt.subplots(figsize=(12, 10))
ax.set_aspect('equal')
all_temps = list(temp_Y.values()) + list(temp_X.values())
vmin, vmax = min(all_temps), max(all_temps)
# Цвета для разных типов ячеек
type_colors = {
CELL_EMPTY: 'lightgray',
CELL_HONEY_SEALED: 'goldenrod',
CELL_HONEY_OPEN: 'gold',
CELL_POLLEN: 'peru',
CELL_BROOD: 'lightcoral',
CELL_BEE: 'red',
}
type_names = {
CELL_EMPTY: 'Пустая',
CELL_HONEY_SEALED: 'Мед (зап.)',
CELL_HONEY_OPEN: 'Мед (откр.)',
CELL_POLLEN: 'Перга',
CELL_BROOD: 'Расплод',
CELL_BEE: 'Пчела',
}
# Сетка '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:
x, y = primal_to_pixel(col, row)
corners = [hex_corner_pointy((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
cell_type = cell_types_Y[(col, row)]
color = type_colors.get(cell_type, 'lightgray')
ax.fill(hex_x, hex_y, color=color, edgecolor='black', linewidth=0.5, alpha=0.8)
# Температура
# temp = temp_Y[(col, row)]
# ax.text(x, y, f"{temp:.2f}", ha='center', va='center', fontsize=6, color='black', fontweight='bold')
# Сетка '^' (только контуры)
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:
x, y = dual_to_pixel(col, row)
corners = [hex_corner_pointy((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
cell_type = cell_types_X[(col, row)]
color = 'gray' if cell_type == CELL_EMPTY else type_colors.get(cell_type, 'gray')
ax.plot(hex_x + (hex_x[0],), hex_y + (hex_y[0],), color=color, linewidth=0.5, alpha=0.5)
ax.axis('off')
# Легенда
legend_elements = []
for type_id, color in type_colors.items():
if type_id in type_names:
legend_elements.append(Patch(facecolor=color, edgecolor='black',
label=type_names[type_id]))
# ax.legend(handles=legend_elements, loc='upper left', bbox_to_anchor=(1.05, 1), borderaxespad=0., fontsize=10)
ax.legend(handles=legend_elements, loc='upper left', bbox_to_anchor=(0.9, 1), borderaxespad=0., fontsize=10)
plt.title(f'Распределение типов ячеек', fontsize=14)
plt.tight_layout()
plt.show()
# 15.2. График температур в центре
plt.figure(figsize=(10, 6))
plt.plot(log_time, log_T_center, 'b-', linewidth=2, label='Температура в центре (0,0)')
plt.axhline(T_ambient, color='r', linestyle='--', label=f'Температура наружного воздуха ({T_ambient}°C)')
plt.xlabel('Время (с)')
plt.ylabel('Температура (°C)')
plt.title(f'Изменение температуры в центральной ячейке кластера', fontsize=14)
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 15.3. График теплового баланса
fig, ax1 = plt.subplots(figsize=(10, 6))
ax1.plot(log_time, log_bee_power, 'orange', linewidth=2, label='Мощность пчел')
ax1.plot(log_time, log_brood_heat, 'green', linewidth=2, label='Тепловыделение расплода')
ax1.set_xlabel('Время (с)')
ax1.set_ylabel('Мощность (мВт)', color='black')
ax1.tick_params(axis='y', labelcolor='black')
ax1.set_ylim(bottom=0)
ax2 = ax1.twinx()
ax2.plot(log_time, log_Q_Y, 'steelblue', linewidth=1.5, label='Потери в улочку Y')
ax2.plot(log_time, log_Q_X, 'crimson', linewidth=1.5, label='Потери в улочку X')
ax2.set_ylabel('Теплопотери (мВт)', color='black')
ax2.tick_params(axis='y', labelcolor='black')
ax2.set_ylim(bottom=0)
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_title(f'Тепловой баланс кластера. Вариант сот {variant}', fontsize=14)
plt.tight_layout()
plt.show()
# 14.6. Температурная карта (тепловая карта) с гексагональной сеткой
fig, ax = plt.subplots(figsize=(12, 10))
ax.set_aspect('equal')
# Найдем диапазон температур для цветовой карты
all_temps_Y = list(temp_Y.values())
vmin_temp, vmax_temp = min(all_temps_Y), max(all_temps_Y)
# Создаем цветовую карту (используем тот же colormap 'hot')
norm = plt.Normalize(vmin=vmin_temp, vmax=vmax_temp)
cmap = plt.cm.hot
# Рисуем шестиугольники для каждой ячейки на стороне Y
for (col, row), temp in temp_Y.items():
x, y = primal_to_pixel(col, row)
corners = [hex_corner_pointy((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
# Цвет в зависимости от температуры
color = cmap(norm(temp))
# Закрашиваем шестиугольник
ax.fill(hex_x, hex_y, color=color, edgecolor='black', linewidth=0.5)
# Подписываем температуру (опционально)
# Можно уменьшить размер шрифта или убрать надписи, если их слишком много
ax.text(x, y, f"{temp:.2f}", ha='center', va='center',
fontsize=6, color='white' if temp < (vmin_temp + vmax_temp) / 2 else 'black',
fontweight='bold')
# Рисуем контуры ячеек на стороне ^ (опционально, для контекста)
for (col, row), temp in temp_X.items():
x, y = dual_to_pixel(col, row)
corners = [hex_corner_pointy((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='gray', linewidth=0.3, alpha=0.5)
ax.axis('off')
# Добавляем цветовую шкалу
sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([])
cbar = plt.colorbar(sm, ax=ax, orientation='vertical', pad=0.02)
cbar.set_label('Температура (°C)')
plt.title(f'Температурная карта кластера (сторона Y). Вариант сот {variant}', fontsize=14)
plt.tight_layout()
plt.show()
# 15.4. Статистика по температурам
print("\n" + "=" * 60)
print("СТАТИСТИКА ПО КЛАСТЕРУ 20x20")
print("=" * 60)
# Температуры по типам ячеек
temp_by_type = {}
count_by_type = {}
for cell in cells:
cell_type = cell_states[cell][0]
temp = T_air[cell_to_index[cell]]
if cell_type not in temp_by_type:
temp_by_type[cell_type] = []
count_by_type[cell_type] = 0
temp_by_type[cell_type].append(temp)
count_by_type[cell_type] += 1
print("\nСредние температуры по типам ячеек:")
for type_id in sorted(temp_by_type.keys()):
if type_id in type_names:
avg_temp = np.mean(temp_by_type[type_id])
min_temp = np.min(temp_by_type[type_id])
max_temp = np.max(temp_by_type[type_id])
count = count_by_type[type_id]
print(f" {type_names[type_id]:20s}: {avg_temp:.2f}°C "
f"(min: {min_temp:.2f}°C, max: {max_temp:.2f}°C, кол-во: {count})")
# Общая статистика
all_temps_flat = [T_air[i] for i in range(len(cells))]
print(f"\nОбщая статистика по всем {len(cells)} ячейкам:")
print(f" Средняя температура: {np.mean(all_temps_flat):.2f}°C")
print(f" Минимальная температура: {np.min(all_temps_flat):.2f}°C")
print(f" Максимальная температура: {np.max(all_temps_flat):.2f}°C")
print(f" Стандартное отклонение: {np.std(all_temps_flat):.2f}°C")
# Тепловыделение
final_bee_power = log_bee_power[-1] if log_bee_power else 0
final_brood_heat = log_brood_heat[-1] if log_brood_heat else 0
print(f"\nТепловыделение в конце моделирования:")
print(f" От пчел: {final_bee_power:.2f} мВт")
print(f" От расплода: {final_brood_heat:.2f} мВт")
print(f" Общее: {final_bee_power + final_brood_heat:.2f} мВт")
Раньше пустых ячеек в полях расплода не было. Пчёлы не грели расплод. Они сидели по периферии полей расплода.
ОтветитьУдалитьНе очень понятно, что Вы хотели сказать. Когда раньше и почему пчелы не грели расплод? Может быть есть ссылка на первоисточник, было бы интересно посмотреть.
ОтветитьУдалить