Гидродинамика гемосорбционной колонки

От одной гранулы к слою

Параметры модели

Зададим свойства сорбента и геометрию колонки с явным указанием единиц измерения. Пакет units переводит все величины в систему СИ и автоматически контролирует размерности в ходе симуляции.

# 1. Свойства сорбента
sorbents <- tibble(
  sorbent_id             = "efferon_bead",
  particle_density       = set_units(1.2, "g/cm^3"), # Плотность влажной гранулы
  bed_packing_factor     = 0.60,                      # Доля объема сорбента (packing factor)
  bed_porosity           = 0.40,                      # Пористость слоя (поровое пространство)
  particle_diameter_mean = set_units(350, "um"),      # Средний диаметр гранул
  particle_diameter_sd   = set_units(80, "um"),       # Стандартное отклонение
  particle_diameter_min = set_units(200, "um"),      # Нижняя граница фракции
  particle_diameter_max = set_units(600, "um")       # Верхняя граница фракции
)

# 2. Геометрия колонки (расчетные параметры вычисляются автоматически)
columns <- tibble(
  column_id       = "column_200ml",
  column_volume   = set_units(200, "mL"),     # Полный объем корпуса
  column_diameter = set_units(4.0, "cm"),     # Внутренний диаметр
  column_height   = set_units(16.0, "cm"),    # Полная высота корпуса
  working_height  = set_units(15.2, "cm"),    # Высота рабочей камеры
  bed_height      = set_units(12.85, "cm"),   # Высота насыпного слоя сорбента
  end_cap_volume  = set_units(5, "mL")        # Объем концевых распределителей
) %>%
  mutate(
    # Площадь поперечного сечения колонки: A = π * D^2 / 4
    column_cross_section_area = (pi / 4) * column_diameter^2,
    
    # Полезный объем рабочей камеры между сетками
    working_volume = column_cross_section_area * working_height,
    
    # Насыпной объем слоя сорбента
    bed_volume = column_cross_section_area * bed_height
  )

# Выбор активного сорбента и колонки для текущего расчета
sorbent <- sorbents %>%
  filter(sorbent_id == "efferon_bead") %>%
  slice(1) %>%
  as.list()

column <- columns %>%
  filter(column_id == "column_200ml") %>%
  slice(1) %>%
  as.list()

Сводная таблица расчетных параметров

Model Parameters Summary

Parameter	Value	Unit
particle_density	1.20	g cm-3
bed_packing_factor	0.60	-
bed_porosity	0.40	-
particle_diameter_mean	350.00	um
column_diameter	4.00	cm
bed_height	12.85	cm
column_cross_section_area	12.57	cm2
working_volume	191.01	cm3
bed_volume	161.48	cm3

Одна частица в потоке

Рассмотрим изолированную сферическую гранулу с диаметром \(d_{\text{particle}}\). Если плотность гранулы \(\rho_{\text{particle}}\) превышает плотность рабочей жидкости \(\rho_{\text{fluid}}\), частица в покоящейся жидкости оседает под действием силы тяжести с учетом выталкивающей силы Архимеда. В восходящем потоке жидкости на гранулу дополнительно действует гидродинамическая сила лобового сопротивления \(F_D\) (Figure 1).

Figure 1: Силы, действующие на изолированную гранулу в восходящем потоке.

Сила сопротивления \(F_D\) зависит от площади поперечного сечения частицы \(A_{\text{particle}}\), относительной скорости движения частицы и жидкости \(v_{\text{terminal}}\), а также коэффициента лобового сопротивления \(C_D\).

NoteМоделирование осаждения с помощью корреляции Шиллера–Наумана

На Figure 1 показан случай, в котором гранула оседает относительно восходящего потока. Сила сопротивления \(F_{\text{drag}}\) направлена вверх: она всегда противостоит относительному движению гранулы и жидкости.

В установившемся режиме оседающего движения погруженный вес частицы сбалансирован гидродинамическим сопротивлением:

\[(\rho_{\text{particle}} - \rho_{\text{fluid}}) \cdot V_{\text{particle}} \cdot g = \frac{1}{2} C_D \cdot \rho_{\text{fluid}} \cdot A_{\text{particle}} \cdot v_{\text{terminal}}^2\]

Где:
\(\rho_{\text{particle}}\) — particle_density (плотность частицы)
\(\rho_{\text{fluid}}\) — fluid_density (плотность жидкости)
\(V_{\text{particle}}\) — particle_volume (объем частицы, \(V_{\text{particle}} = \pi \cdot d_{\text{particle}}^3 / 6\))
\(g\) — gravity (ускорение свободного падения)
\(C_D\) — drag_coefficient (коэффициент лобового сопротивления)
\(A_{\text{particle}}\) — particle_cross_section_area (площадь сечения частицы, \(A_{\text{particle}} = \pi \cdot d_{\text{particle}}^2 / 4\))
\(v_{\text{terminal}}\) — terminal_velocity (установившаяся скорость свободного осаждения)

Для сферической частицы коэффициент сопротивления вычислим по корреляции Шиллера–Наумана (Schiller and Naumann 1933):

\[Re_{\text{particle}} = \frac{\rho_{\text{fluid}} \cdot v_{\text{terminal}} \cdot d_{\text{particle}}}{\mu_{\text{fluid}}}\]

\[C_D = \frac{24}{Re_{\text{particle}}} \left(1 + 0.15 \cdot Re_{\text{particle}}^{0.687}\right)\]

Где:
\(Re_{\text{particle}}\) — particle_reynolds_number (число Рейнольдса для частицы)
\(d_{\text{particle}}\) — particle_diameter (диаметр частицы)
\(\mu_{\text{fluid}}\) — fluid_viscosity (динамическая вязкость жидкости)

Уравнения задают скорость осаждения \(v_{\text{terminal}}(d_{\text{particle}})\) неявно.

NoteЧисленное решение уравнения осаждения (Солвер)

Поскольку коэффициент сопротивления \(C_D\) зависит от числа Рейнольдса \(Re\), которое само содержит искомую скорость \(v_{\text{terminal}}\), система решается методом простой итерации (неподвижной точки). В качестве начального приближения используется скорость Стокса.

calc_terminal_velocity <- function(diameter, 
                                   particle_density, 
                                   fluid_density, 
                                   fluid_viscosity, 
                                   gravity = set_units(9.81, "m/s^2")) {
  
  # 1. Снимаем единицы, приводя всё к СИ (m, kg/m^3, Pa*s, m/s^2)
  d_m    <- drop_units(set_units(diameter, "m"))
  rho_p  <- drop_units(set_units(particle_density, "kg/m^3"))
  rho_f  <- drop_units(set_units(fluid_density, "kg/m^3"))
  mu_val <- drop_units(set_units(fluid_viscosity, "Pa*s"))
  g_val  <- drop_units(set_units(gravity, "m/s^2"))
  
  delta_rho <- rho_p - rho_f
  
  # 2. Векторное начальное приближение по Стоксу
  v <- (delta_rho * g_val * d_m^2) / (18 * mu_val)
  
  # 3. Векторизованный цикл (считает ВСЕ 40 000 элементов одновременно)
  for (i in seq_len(100)) {
    re <- (rho_f * v * d_m) / mu_val
    re <- pmax(re, 1e-12) # pmax векторизован
    
    cd <- (24 / re) * (1 + 0.15 * re^0.687)
    v_next <- sqrt((4 * g_val * delta_rho * d_m) / (3 * rho_f * cd))
    
    # Проверка сходимости для всего вектора
    if (max(abs(v_next - v)) < 1e-10) break
    v <- (v + v_next) / 2
  }
  
  # 4. Навешиваем единицы СИ 1 раз на готовый вектор
  set_units(v, "m/s")
}

Чем частица меньше и проще, тем легче она уносится вверх восходящим потоком. На Figure 2 представлена карта областей уноса и оседания гранул.

fluid_density   <- set_units(1025, "kg/m^3")
fluid_viscosity <- set_units(3.5, "mPa*s")

# 1. Функция трансформации с АВТО-ГЕНЕРАТОРОМ идеальных брейксов
sqrt_mono <- trans_new(
  name      = "sqrt_mono",
  transform = function(x) sqrt(pmax(0, x)),
  inverse   = function(x) ifelse(x < 0, x, x^2),
  # Генерирует 'pretty' числа в пространстве корня и возводит в квадрат.
  # Дает идеально равные физические интервалы на экране под ЛЮБЫЕ данные!
  breaks    = function(limits) pretty(sqrt(pmax(0, limits)), n = 5)^2
)

# Сетка параметров (можно менять 10 на 5, 25, 100 — график сам подстроится)
particle_map <- crossing(
  liquid_velocity = seq(0.01, sqrt(10), length.out = 500)^2 |> set_units("mm/s"),
  diameter        = seq(100, 1000, length.out = 500) |> set_units("um")
) |>
  mutate(
    settling_velocity   = calc_terminal_velocity(diameter, sorbent$particle_density, fluid_density, fluid_viscosity),
    particle_velocity   = liquid_velocity - settling_velocity,
    liquid_velocity_num = drop_units(liquid_velocity),
    diameter_num        = drop_units(diameter),
    particle_v_num      = particle_velocity |> set_units("mm/s") |> drop_units()
  )

# Коэффициент пересчета в объемный расход Q (мл/мин)
volume_flow_per_mm_s <- (set_units(1, "mm/s") * column$column_cross_section_area) |>
  set_units("mL/min") |>
  drop_units()

# Динамический расчет координат подписей в % от максимальной скорости
x_max <- max(particle_map$liquid_velocity_num)

# Построение 100% адаптивного графика
ggplot(particle_map, aes(liquid_velocity_num, diameter_num)) +
  geom_tile(aes(fill = particle_v_num)) +
  geom_contour(aes(z = particle_v_num), breaks = 0, color = "black", linewidth = 0.8) +
  
  # Аннотации динамически привязаны к % от длины оси
  annotate("text", x = x_max * 0.1, y = 820, label = "Settling\n(particle velocity < 0)") +
  annotate("text", x = x_max * 0.55, y = 260, label = "Upward transport\n(particle velocity > 0)") +
  
  scale_fill_steps2(
    low      = "#024eab", 
    mid      = "white", 
    high     = "#ab0202", 
    midpoint = 0,
    limits   = c(-16, 9),
    #oob = scales::squish,
    n.breaks = 9,  
    name     = "Particle velocity,\nmm/s",
    guide    = guide_colorsteps(barheight = unit(6, "cm"))
  ) +
  
  # Полностью автонастраиваемые оси
  scale_x_continuous(
    trans    = sqrt_mono,
    sec.axis = sec_axis(
      ~ . * volume_flow_per_mm_s, 
      name   = "Volumetric flow rate Q, mL/min"
    )
  ) +
  
  labs(
    x = expression("Liquid superficial velocity " * U * ", mm/s"),
    y = "Particle diameter, µm"
  )

Figure 2: Направление движения изолированной гранулы в восходящем потоке. Черная линия соответствует нулевой скорости частицы (\(v_{\text{particle}} = 0\)).

References

Schiller, L., and A. Naumann. 1933. “Über Die Grundlegenden Berechnungen Bei Der Schwerkraftaufbereitung.” Zeitschrift Des Vereines Deutscher Ingenieure 77: 318–20. https://cir.nii.ac.jp/crid/1573950399224895872.