В заметке Пчела в ячейке была рассмотрена ячейка с пчелой. Ее взаимодействие с окружающей средой было представлено, как эквивалентное конвективное охлаждение воздухом общей (суммарной) внешней поверхности ячейки - это существенное упрощение. Оно не учитывает теплоемкость и теплопроводность восковых стенок соседних ячеек. А они, как будет видно далее, существенно изменяют текущую термодинамику, как в самой ячейке с пчелой, так и в соседних, и не только самых близких. Дело в том, что теплоемкость и теплопроводность воска значительно больше, чем у воздуха: теплоемкость воздуха в пустой ячейке равна 0.367 мДж/°С, а теплоемкость стенки (при учете половины толщины, т.к. вторая половина относится к соседу) - 68.178 мДж/°С, т.е. в 186 раз больше. Кроме того, соседние ячейки могут быть в разных состояниях (мед, расплод, пчела, пустая, перга, рамка), поэтому и взаимодействие с ними будет разным. Таким образом, для создания полноценной функциональной модели необходимо, во-первых, задать модель каждой соседней ячейки, а не только ячейки с пчелой, и, во-вторых, рассматривать теплообмен для каждой грани отдельно, определяя с какой именно гранью из соседней ячейки контакт.
Итак, перейдем от общей площади поверхности ячейки с пчелой к ее отдельным граням: шести боковым и трем донным. Это означает, что вместо уравнения $$ C_{wall}\frac{dT_{wall}}{dt}=p_{r\_wall}+\frac{T_{ai}-T_{wall}}{r_{ai\_wall}}-\frac{T_{wall}-T_{air}}{r_{wall\_air}} $$ нужно вводить девять частных уравнений для граней, каждое из которых имеет вид:$$C_{wall}[k]\frac{dT_{wall}[k]}{dt}=p_{r\_wall}[k]+\frac{T_{ai}-T_{wall}[k]}{r_{ai\_wall}[k]}-\frac{T_{wall}-T_{neighbor}[k]}{r},$$где k - индекс грани (от 0 до 8), r - тепловое сопротивление грани (принимает одно из двух значений в зависимости от того, боковая грань или донная), \(T_{neighbor}[k]\) - температура грани со стороны соседа.
Итого, получится система из 13 уравнений для ячейки с пчелой. Для модели пустых ячеек - это система из 10 уравнений, т.к. там нет уравнений для головы, торакса и брюшка пчелы. Таким образом, проходя по всем ячейкам с шагом, например, 0.01 секунды, можно увидеть динамику изменения температуры в любой ячейке, а в конце - итоговое распределение температуры по ячейкам.
В упомянутой упрощенной модели было показано, что в светлом соте пчеле было немного теплее, и войдя в такую ячейку при температуре окружающего воздуха 12°С она нагревала воздух в ячейке до 18°С. Поэтому здесь тоже начнем моделирование для светлого сота и исходной температуры 12°С. Пчелу разместим в центре кластера 7х7 ячеек на стороне Y.
Видим, что температура воздуха и стенки достигает установившегося состояния не за 4 минуты, а за 2; температура частей пчелы - не за полчаса, а за 5 минут. Воздух в ячейке нагрелся всего до 12.87°С (вместо 18°С).
Итоговое распределение температуры по ячейкам:
А теперь, также, как и в упрощенной модели, рассмотрим вариант темных сот:
Из итогового распределения видно, что пчела больше всего, естественно, согревает свою ячейку. Температура в соседних боковых ячейках поднялась больше, чем в соседних донных (на 0.35 и 0.16°С, соответственно), что тоже вполне естественно, т.к. площадь донных граней (ромбов) примерно в 4 раза меньше боковых, а их толщина больше. Температура немного поднялась во втором и третьем от пчелы ряду соседей: в боковом втором на 0.12-0.17°С, третьем - на 0.03-0.07°С, а в донном втором - на 0.09-0.11°С, третьем - на 0.03-0.06°С. В последующих рядах повышение температуры пренебрежимо мало.
Видим, что сот выполняет роль аккумулятора тепла, но в то же время и радиатора охлаждения. Т.е. часть тепла пчелы распределяется по соту, а часть уходит в улочку. Посмотрим как соотносятся эти доли.
Видно, что компактное размещение нескольких пчел позволяет заметно повысить температуру в ячейках с пчелами и в близлежащей зоне их размещения. 
Во-первых, видим, что в первые несколько секунд после входа пчелы в холодную ячейку мощность термогенеза падает с исходных 20.27 мВт до 18.25, а потом в течение нескольких минут устанавливается около 18.88 мВт. Т.е. в ячейке пчеле теплее, и она тратит примерно на 7% меньше энергии, чем вне ее.
Во-вторых, около 2.6 мВт из этих 18.88 уходит в улочки (пока не учитываем, что при повышении температуры в улочке это значение будет уменьшаться, но в дальнейшем при моделировании улочки это будет учтено). При этом более 3/4 уходит в улочку на стороне пчелы (в данном случае на стороне Y), а чуть менее четверти - на сторону ^.
Однако, продолжим моделирование и посмотрим, что показывает модель при нескольких пчелах. Добавим три пчелы на стороне Y и увеличим кластер до 11х11.
Таким образом, разработанная модель теплового взаимодействия ячеек показывает достаточно правдоподобные результаты.
Несмотря на то, что полученные данные существенно отличаются от результатов моделирования на упрощенной модели, не учитывавшей тепловое взаимодействие ячеек, я бы не стал сейчас отказываться от них и убирать заметку "Пчела в ячейке" из блога. Это был важный и, возможно, необходимый этап в становлении и развитии модели.
Когда я писал (и читал у некоторых авторов), что светлый сот теплее - это противоречило моим убеждениям: то, что темный сот теплее и его предпочитают пчелы - это же притча во языцех. Результаты полученные на модели подтверждают - старый темный сот теплее нового светлого.
🐝 Python-скрипт
""" Программа к заметке 'Тепловое взаимодействие ячеек в соте' """
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
import math
import time
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# ==========================
# 1. Физические параметры
# ==========================
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
k_hd, k_ab = 0.065, 0.042
ρ_air, k_air, c_air = 1.205, 0.025, 1005
T_ambient = 12.0
variant = 2 # 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 # шаг решетки
pitch = thread - L_w[variant-1] # расстояние между гранями
h_cell = 12e-3 # высота (глубина) ячейки
dh = 1e-3 # углубление
a = pitch / 2. / np.cos(np.radians(30)) # сторона внутреннего шестиугольника
A_lateral = a * (h_cell - dh/2.) # площадь боковой грани
d1, d2 = pitch, np.sqrt(a ** 2 + dh ** 2) # диагонали ромба
A_bottom = d1 * d2 / 2. # площадь ромба
ρ_wax, λ_wax, c_wax = 960, λ_w[variant-1], c_w[variant-1]
M_lateral = ρ_wax * A_lateral * L_w[variant-1]/2
M_bottom = ρ_wax * A_bottom * L_b[variant-1]/2
d_cell = pitch
sigma = 5.67e-8
t0 = -273.15
eps = 0.22
eps_w = [0.22, 0.5, 0.9]
eps_wall = eps_w[variant-1]
ke0 = 3e-5
Pmin = 0.2e-3
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
def metabolic_power(T):
power = (29.133 - 0.739 * T) * 1e-3
return max(power, Pmin)
R_lateral = L_wall / (λ_wax * A_lateral)
R_bottom = L_bottom / (λ_wax * A_bottom)
# ==========================
# 2. Геометрия и теплоёмкости
# ==========================
a = pitch / (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'Теплоемкости (Дж/°С)\n\t воздуха пустой ячейки\t\t {C_air_cell:.6f}\n\t пчелы\t\t\t\t\t {C_bee:.6f}'
f'\n\t воздуха в ячейке с пчелой\t {C_air_bee:.6f}\n\t воска ячейки \t {sum(C_wall):.6f}'
f'\nТепловые сопротивления (м²·°С/Вт)\n\t боковой стенки {R_lateral:.3f}\n\t донной грани {R_bottom:.3f}')
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]
# ==========================
# 3. Ваши правила соседства
# ==========================
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
# ==========================
# 4. Генерация кластера 11x11
# ==========================
cells = []
cols, rows = 11, 11
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))
bee_cells = [(0, 'Y', 0, 0), (0, 'Y', 0, 1), (0, 'Y', 0, 2), (0, 'Y', 2, -1)]
# bee_cells = [(0, 'Y', 0, 0)]
# Индексы для быстрого доступа
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)
# Типы ячеек
is_bee = [cell in bee_cells for cell in cells]
# ==========================
# 5. Кэширование связей граней
# ==========================
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)
# ==========================
# 6. Модель пчелы с 9 гранями (оставляем fsolve)
# ==========================
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_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
# ==========================
# 7. Модель пустой ячейки с 9 гранями
# ==========================
def empty_cell_update(T_air_prev, T_wall_prev, neighbor_T_list, dt):
Q_air_to_wall = sum((T_air_prev - T_wall_prev[k]) / R_ai_wall[k] for k in range(9))
Q_conv_entrance = h_entrance * A_entrance * (T_air_prev - T_ambient)
T_air_new = T_air_prev + dt / C_air_cell * (-Q_air_to_wall - Q_conv_entrance)
T_wall_new = []
for k in range(9):
R = R_lateral if k < 6 else R_bottom
A = A_lateral if k < 6 else A_bottom
# Тепло от соседей
Q_from_neighbor = (neighbor_T_list[k] - T_wall_prev[k]) / R
# Тепло от воздуха в ячейке
Q_from_air = (T_air_prev - T_wall_prev[k]) / R_ai_wall[k]
# Излучение наружу (только для граничных граней)
# Определяем, является ли грань граничной
if k < 6: # Боковые грани
# Проверяем, есть ли сосед через эту грань
# Для упрощения можно считать все боковые грани наружными
Q_rad = eps_wall * sigma * A * ((T_wall_prev[k] - t0) ** 4 - (T_ambient - t0) ** 4)
else: # Дно
Q_rad = eps_wall * sigma * A * ((T_wall_prev[k] - t0) ** 4 - (T_ambient - t0) ** 4)
dT = dt / C_wall[k] * (Q_from_air + Q_from_neighbor - Q_rad)
T_wall_new.append(T_wall_prev[k] + dT)
return T_air_new, T_wall_new
# ==========================
# 8. Начальные условия
# ==========================
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_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
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_avg_init = (m_thorax * Tth_init + m_head * Thd_init + m_ab * Tab_init) / (m_thorax + m_head + m_ab)
Power_init = metabolic_power(T_ambient)
# Инициализация (через индексы)
T_air = [T_ambient] * len(cells)
T_wall = [[T_ambient] * 9 for _ in range(len(cells))]
T_bee_states = {} # cell_addr -> state
for cell in bee_cells:
T_bee_states[cell] = T_bee_init.copy()
# ==========================
# 9. Моделирование
# ==========================
dt = 0.01
time_total = 600.0
steps = int(time_total / dt)
log_time = [0]
log_T_center = [Tth_init]
log_T_hd = [Thd_init]
log_T_ab = [Tab_init]
log_T_ai_bee = [T_ambient]
log_T_wall_bee = [T_ambient]
log_T_bee_avg = [T_bee_avg_init]
log_Power = [Power_init * 1e3]
log_Q_Y = [1e-9]
log_Q_X = [0.]
cumulative_heat_generated = 0.0
cumulative_heat_loss = 0.0
start_time = time.time()
for step in range(steps):
t = step * dt
Q_entrance_Y = 0.0
Q_entrance_X = 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]
# Сбор температур за гранями
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]
if is_bee[j]:
neighbor_T_list.append(T_bee_states[neighbor_cell][4 + n_face])
else:
neighbor_T_list.append(T_wall[j][n_face])
if is_bee[i]:
Tai = T_bee_states[cell][3]
P_met = metabolic_power(Tai)
cumulative_heat_generated += P_met * dt # Накопленное тепло
# Также нужно отслеживать изменение внутренней энергии пчелы
if cell == bee_cells[0]: # Для первой пчелы
# Сохраняем предыдущие температуры
global T_bee_prev
if 'T_bee_prev' not in globals():
T_bee_prev = T_bee_states[cell].copy()
else:
# Расчет изменения энергии
dE_bee = (C_hd * (T_bee_states[cell][0] - T_bee_prev[0]) +
C_th * (T_bee_states[cell][1] - T_bee_prev[1]) +
C_ab * (T_bee_states[cell][2] - T_bee_prev[2]) +
C_air_bee * (T_bee_states[cell][3] - T_bee_prev[3]))
T_bee_prev = T_bee_states[cell].copy()
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_
if side == 'Y':
Q_entrance_Y += Q_cell + Q_ab
else:
Q_entrance_X += Q_cell + Q_ab
else:
T_air_new, T_wall_new = empty_cell_update(T_air[i], T_wall[i], neighbor_T_list, dt)
new_T_air[i] = T_air_new
new_T_wall[i] = T_wall_new
Q_cell = hA * (T_air[i] - T_ambient)
if side == 'Y':
Q_entrance_Y += Q_cell
else:
Q_entrance_X += Q_cell
T_air = new_T_air
T_wall = new_T_wall
# Логирование (первая пчела)
first_bee = bee_cells[0]
if step % 100 == 0:
log_time.append(t)
log_T_center.append(T_bee_states[first_bee][1])
log_T_hd.append(T_bee_states[first_bee][0])
log_T_ab.append(T_bee_states[first_bee][2])
log_T_ai_bee.append(T_bee_states[first_bee][3])
log_T_wall_bee.append(T_bee_states[first_bee][4])
T_bee_avg = (m_head * T_bee_states[first_bee][0] + m_thorax * T_bee_states[first_bee][1] +
m_ab * T_bee_states[first_bee][2]) / (m_head + m_thorax + m_ab)
log_T_bee_avg.append(T_bee_avg)
log_Power.append(metabolic_power(T_bee_states[first_bee][3]) * 1000)
log_Q_Y.append(Q_entrance_Y * 1000) # в мВт
log_Q_X.append(Q_entrance_X * 1000) # в мВт
if step % 1000 == 0:
print(f"Шаг {step}, время на шаг: {(time.time() - start_time) / 100:.4f} с")
start_time = time.time()
print(f"\nИтоговый баланс за {time_total} с:")
print(f"Суммарное тепловыделение: {cumulative_heat_generated:.6f} Дж")
print(f"Температура пчелы изменилась с {T_bee_avg_init:.2f}°C до {log_T_bee_avg[-1]:.2f}°C")
# ==========================
# 10. Сбор финальных температур (воздух в ячейках)
# ==========================
final_temperatures = {}
for i, cell in enumerate(cells):
final_temperatures[cell] = T_air[i] # ← температура ВОЗДУХА
def check_energy_balance(step_num):
"""Проверка полного энергетического баланса системы"""
# 1. Общее тепловыделение от всех пчел
total_heat_generation = 0.0
for cell in bee_cells:
if cell in T_bee_states:
Tai = T_bee_states[cell][3]
P_met = metabolic_power(Tai)
total_heat_generation += P_met * dt # Дж за шаг
# 2. Тепло, аккумулированное в компонентах системы
heat_stored_total = 0.0
# а) В пчелах
for cell in bee_cells:
if cell in T_bee_states:
# Кинетическая энергия изменения температуры пчелы
T_curr = T_bee_states[cell]
# Нужно сохранять предыдущие температуры для расчета dT
# Пока пропустим, добавим позже
# 3. Теплопотери через ВСЕ границы кластера
total_heat_loss = 0.0
# а) Через входы ячеек
for i, cell in enumerate(cells):
if is_bee[i]:
if cell in T_bee_states:
Tai = T_bee_states[cell][3]
Tab = T_bee_states[cell][2]
Q_cell = hA_ * (Tai - T_ambient) * dt
Q_ab = (Tab - T_ambient) / R_ab_ * dt
total_heat_loss += Q_cell + Q_ab
else:
Q_cell = hA * (T_air[i] - T_ambient) * dt
total_heat_loss += Q_cell
# б) Через излучение ВСЕХ наружных стенок
for i, cell in enumerate(cells):
for face_id in range(9):
conn = face_connections[i][face_id]
if conn is None: # Граничная грань
A = A_lateral if face_id < 6 else A_bottom
T_wall_val = T_wall[i][face_id]
Q_rad = eps_wall * sigma * A * (
(T_wall_val - t0) ** 4 - (T_ambient - t0) ** 4) * dt
total_heat_loss += Q_rad
print(f"\n=== Энергетический баланс (шаг {step_num}) ===")
print(f"Тепловыделение за шаг: {total_heat_generation * 1000:.6f} мДж")
print(f"Теплопотери за шаг: {total_heat_loss * 1000:.6f} мДж")
print(f"Дисбаланс: {(total_heat_generation - total_heat_loss) * 1000:.6f} мДж")
# Вызов проверки баланса
check_energy_balance(steps)
def check_balance_simple():
"""Упрощенная проверка баланса за всё время моделирования"""
# Интегрируем мощность термогенеза
total_energy_generated = 0.0
# Для простоты - только для первой пчелы
if bee_cells:
first_bee = bee_cells[0]
if first_bee in T_bee_states:
# Используем логированные значения (каждые 100 шагов)
for i in range(1, len(log_time)):
dt_log = log_time[i] - log_time[i - 1]
# Средняя мощность за интервал
P_avg = (log_Power[i] + log_Power[i - 1]) / 2 / 1000 # Вт
total_energy_generated += P_avg * dt_log
# Расчет изменения внутренней энергии системы
# 1. Изменение температуры пчелы
if bee_cells:
first_bee = bee_cells[0]
if first_bee in T_bee_states:
T_final = T_bee_states[first_bee]
T_avg_final = (m_head * T_final[0] + m_thorax * T_final[1] + m_ab * T_final[2]) / (m_head + m_thorax + m_ab)
delta_E_bee = C_bee * (T_avg_final - T_bee_avg_init)
print(f"\n=== Упрощенный баланс за 600 с ===")
print(f"Выработано тепла: {total_energy_generated:.6f} Дж")
print(f"Изменение энергии пчелы: {delta_E_bee:.6f} Дж")
print(f"Разница (ушло в окружение): {total_energy_generated - delta_E_bee:.6f} Дж")
check_balance_simple()
# После основного цикла
def analyze_cluster_temperatures():
# Средняя температура воздуха в ячейках
T_air_avg = np.mean(T_air)
# Средняя температура стенок
T_wall_flat = [T_wall[i][k] for i in range(len(cells)) for k in range(9)]
T_wall_avg = np.mean(T_wall_flat)
# Температура в центральной области (3x3x2)
center_cells = []
for cell in cells:
scomb, side, col, row = cell
if abs(col) <= 1 and abs(row) <= 1:
center_cells.append(cell)
T_air_center = np.mean([T_air[cell_to_index[cell]] for cell in center_cells])
print(f"\n=== Температурный анализ кластера ===")
print(f"Средняя температура воздуха в кластере: {T_air_avg:.3f} °C")
print(f"Средняя температура стенок кластера: {T_wall_avg:.3f} °C")
print(f"Средняя температура в центральной области (3x3x2): {T_air_center:.3f} °C")
print(f"Повышение относительно начальной: {T_air_avg - T_ambient:.3f} °C")
# Теплосодержание кластера относительно начального состояния
Q_stored_wax = 0.0
for i in range(len(cells)):
for k in range(9):
Q_stored_wax += C_wall[k] * (T_wall[i][k] - T_ambient)
Q_stored_air = 0.0
for i in range(len(cells)):
if is_bee[i]:
C_air_cell_val = C_air_bee
else:
C_air_cell_val = C_air_cell
Q_stored_air += C_air_cell_val * (T_air[i] - T_ambient)
print(f"\n=== Накопленная энергия ===")
print(f"В воске кластера: {Q_stored_wax:.6f} Дж")
print(f"В воздухе кластера: {Q_stored_air:.6f} Дж")
print(f"Всего: {Q_stored_wax + Q_stored_air:.6f} Дж")
# Сколько времени нужно для нагрева на 1°C
time_for_1C = (Q_stored_wax + Q_stored_air) / (T_air_avg - T_ambient) / 0.0167 # при 16.7 мВт
print(f"\nОценочное время нагрева на 1°C при текущей мощности: {time_for_1C:.1f} с")
analyze_cluster_temperatures()
# ==========================
# 11. Визуализация
# ==========================
plt.figure(figsize=(7, 5))
plt.plot(log_time, log_T_center, 'r-', label='Торакс', linewidth=1.5)
plt.plot(log_time, log_T_hd, 'b-', label='Голова', linewidth=1.5)
plt.plot(log_time, log_T_ab, 'g-', label='Брюшко', linewidth=1.5)
plt.plot(log_time, log_T_ai_bee, 'c--', label='Воздух в ячейке', linewidth=1.5)
plt.plot(log_time, log_T_wall_bee, 'm--', label='Стенка ячейки', linewidth=1.5)
plt.plot(log_time, log_T_bee_avg, 'k:', label='Средняя температура пчелы', linewidth=2)
plt.axhline(T_ambient, color='k', linestyle=':', label=f'Наружный воздух ({T_ambient}°C)')
plt.xlabel('Время (с)')
plt.ylabel('Температура (°C)')
plt.title('Пчела в кластере. Вариант ' + str(variant))
plt.legend()
plt.grid(True, which='both', linestyle=':')
plt.tight_layout()
plt.show()
fig, ax1 = plt.subplots(figsize=(7, 5))
# Мощность термогенеза
ax1.plot(log_time, log_Power, 'orange', linewidth=2, label='Мощность термогенеза')
ax1.set_xlabel('Время (с)')
ax1.set_ylabel('Мощность (мВт)', color='orange')
ax1.tick_params(axis='y', labelcolor='orange')
# Теплопотери в улочки
ax2 = ax1.twinx()
ax2.plot(log_time, log_Q_Y, 'steelblue', linewidth=2, label='Потери в улочку Y')
ax2.plot(log_time, log_Q_X, 'crimson', linewidth=2, label='Потери в улочку ^')
ax2.set_ylabel('Теплопотери (мВт)', color='k')
log_Q_ratio = [y / x if x > 1e-6 else 0 for x, y in zip(log_Q_Y, log_Q_X)]
ax_inset = inset_axes(ax2, width="33%", height="33%", loc='lower center', borderpad=2)
ax_inset.plot(log_time, log_Q_ratio, 'purple', linewidth=1.5)
ax_inset.set_title('Q_^ / Q_Y', fontsize=10)
ax_inset.grid(True, linestyle=':')
ax_inset.set_ylim(0., 1.25)
# Общая легенда
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, which='both', linestyle=':')
ax1.set_title('Баланс мощности: генерация vs потери в улочки. Вариант ' + str(variant))
plt.grid(True, which='both', linestyle=':')
# plt.tight_layout()
plt.show()
# ==========================
# 12. 2D-визуализация
# ==========================
temp_Y = {}
temp_X = {}
for cell in cells:
scomb, side, col, row = cell
if side == 'Y':
temp_Y[(col, row)] = final_temperatures[cell]
else:
temp_X[(col, row)] = final_temperatures[cell]
size = 1.0
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=(7, 5))
ax.set_aspect('equal')
font_size = 9
all_temps = list(temp_Y.values()) + list(temp_X.values())
vmin, vmax = min(all_temps), max(all_temps)
# Сетка 'Y'
for col in range(-int(cols/2), int(cols/2) + 1):
for row in range(-int(rows/2), int(rows/2) + 1):
if (col, row) in temp_Y:
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 = plt.cm.viridis((temp_Y[(col, row)] - vmin) / (vmax - vmin))
ax.fill(hex_x, hex_y, color=color, edgecolor='steelblue', linewidth=1.2, alpha=0.8)
if (col, row) == (0, 0):
ax.fill(hex_x, hex_y, facecolor='none', edgecolor='black', linewidth=3)
ax.text(x, y, f"{temp_Y[(col, row)]:.2f}", ha='center', va='center',
fontsize=font_size-2, color='b')
# Сетка '^'
for col in range(-int(cols/2), int(cols/2) + 1):
for row in range(-int(rows/2), int(rows/2) + 1):
if (col, row) in temp_X:
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='crimson', linewidth=1.2)
ax.text(x, y, f"{temp_X[(col, row)]:.2f}", ha='center', va='center',
fontsize=font_size-2, color='crimson')
ax.axis('off')
sm = plt.cm.ScalarMappable(cmap=plt.cm.viridis, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array([])
plt.colorbar(sm, ax=ax, shrink=0.8).set_label('Температура (°C)', fontsize=14)
plt.title('Температура воздуха в ячейках. Вариант ' + str(variant) +
'\nСиняя сетка — сторона "Y", Красная — сторона "^"')
plt.tight_layout()
plt.show()
print(f"\nВремя моделирования {time_total:.0f} секунд; вариант {variant}\nФинальные результаты для первой пчелы:")
print(f"Торакс: \t\t{log_T_center[-1]:.2f}/{Tth_init:.2f}°C")
print(f"Средняя пчела: \t{log_T_bee_avg[-1]:.2f}/{T_bee_avg_init:.2f}°C")
print(f"Воздух: \t\t{log_T_ai_bee[-1]:.2f}/{T_ambient:.2f}°C")
print(f"Мощность: \t\t{log_Power[-1]:.2f}/{Power_init * 1e3:.2f} мВт")
Комментарии
Отправить комментарий