version control

Deliverable Type: RAPIDS™ Technical Report
Specialized Label: Efficiency Improvement Study
Version Level: Advanced
Release Status: Client Release — Issued for Review & Comment

Resumen ejecutivo

Este documento presenta un modelo analítico completo para caracterización hidráulica de bombas centrífugas horizontales Summit ESP ST2500CCW, con validación contra datos operativos del campo CPE6. El modelo permite:

Resultados clave: Ventana preventiva de 5 minutos desde precursores espectrales hasta falla mecánica, consistente con análisis vibracional integrado.

1. Marco teórico

1.1 Ecuación característica

La curva H-Q de una bomba centrífuga sigue una relación parabólica:

\[ H(Q) = H_0 - A \cdot Q - B \cdot Q^2 \]

Donde:

  • \(H_0\) = Altura de cierre (shutoff head), presión cuando \(Q = 0\)
  • \(A\) = Pérdidas por fricción hidráulica (lineal con Q)
  • \(B\) = Pérdidas por recirculación y choque (cuadrática con Q)

1.2 Leyes de afinidad

Para escalado entre frecuencias:

\[ \frac{Q_2}{Q_1} = \frac{N_2}{N_1} \qquad \frac{H_2}{H_1} = \left(\frac{N_2}{N_1}\right)^2 \qquad \frac{P_2}{P_1} = \left(\frac{N_2}{N_1}\right)^3 \]

1.3 Eficiencia y BEP

La eficiencia hidráulica presenta un máximo en el Best Efficiency Point (BEP):

\[ \eta(Q) = \eta_{max} \cdot \left[1 - k \left(\frac{Q - Q_{BEP}}{Q_{BEP}}\right)^2\right] \]

Operar fuera del rango \(\pm15\%\) del BEP incrementa desgaste y vibración.

2. Datos de entrada

2.1 Especificaciones del fabricante

# Tabla de diseño Halliburton para ST2500CCW
curva_fabricante <- data.frame(
  hz = c(30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60),
  Q_design = c(32086, 34225, 36364, 38503, 40642, 42781, 44920, 47059, 
               49198, 51337, 53476, 55615, 57754, 59893, 62032, 64171),
  H_target = c(436, 497, 561, 629, 700, 776, 856, 939, 1026, 1117, 
               1212, 1311, 1414, 1521, 1631, 1746),
  H_minflow = c(469, 534, 602, 676, 753, 834, 919, 1009, 1103, 1201, 
                1303, 1410, 1520, 1635, 1753, 1876),
  H_maxflow = c(279, 317, 358, 401, 447, 495, 546, 599, 655, 713, 
                774, 837, 902, 970, 1041, 1114),
  Q_minflow = c(19681, 20993, 22305, 23617, 24929, 26241, 27553, 28865, 
                30177, 31489, 32802, 34114, 35425, 36737, 38050, 39361),
  Q_maxflow = c(53530, 57098, 60667, 64235, 67804, 71373, 74941, 78510, 
                82079, 85647, 89216, 92785, 96353, 99922, 103491, 107059)
)

# Vista resumida
curva_fabricante %>%
  select(hz, Q_design, H_target, Q_minflow, Q_maxflow) %>%
  filter(hz %in% c(30, 40, 50, 52, 56, 60)) %>%
  gt_rapids(
    title = "Especificaciones de Diseño",
    subtitle = "Puntos clave para frecuencias seleccionadas",
    footer = "**Fuente:** Halliburton Design Service Document | Campo CPE6-360K"
  ) %>%
  fmt_number(columns = c(Q_design, Q_minflow, Q_maxflow), decimals = 0) %>%
  fmt_number(columns = H_target, decimals = 0) %>%
  cols_label(
    hz = "Frecuencia (Hz)",
    Q_design = "Q BEP (bpd)",
    H_target = "H BEP (psi)",
    Q_minflow = "Q mín (bpd)",
    Q_maxflow = "Q máx (bpd)"
  )
Especificaciones de Diseño
Puntos clave para frecuencias seleccionadas
Frecuencia (Hz) Q BEP (bpd) H BEP (psi) Q mín (bpd) Q máx (bpd)
30 32,086 436 19,681 53,530
40 42,781 776 26,241 71,373
50 53,476 1,212 32,802 89,216
52 55,615 1,311 34,114 92,785
56 59,893 1,521 36,737 99,922
60 64,171 1,746 39,361 107,059
Fuente: Halliburton Design Service Document | Campo CPE6-360K

2.2 Datos operativos

# Cargar datos históricos (ejemplo simplificado)
data_raw <- read.csv2(
  "~/R/jre_2025/0.5_AasS/data/hps_min_feb_mar.csv",
  header = TRUE,
  stringsAsFactors = FALSE,
  encoding = "UTF-8",
  check.names = FALSE
) %>%
  mutate(DateTime = dmy_hms(DateTime, tz = "UTC")) %>%
  mutate(across(-DateTime, ~ suppressWarnings(as.numeric(str_replace(., ",", ".")))))

# Función para procesar cada bomba
procesar_bomba <- function(data, bomba_id) {
  tags <- list(
    speed = paste0("P042", bomba_id, "_VSD.Val_SpeedFdbk"),
    flow = paste0("413_FI_042", bomba_id, ".PV"),
    power = paste0("P_042", bomba_id, "_VSD.Power"),
    pressure = paste0("413_PI_042", bomba_id, "_4.PV")
  )
  
  data %>%
    arrange(DateTime) %>%
    transmute(
      DateTime,
      bomba = bomba_id,
      speed_hz = .data[[tags$speed]],
      flow_bpd = .data[[tags$flow]],
      power_kw = .data[[tags$power]],
      pressure_psi = if (tags$pressure %in% names(.)) .data[[tags$pressure]] else NA_real_
    )
}

# Procesar bombas I, J, K
data_operativo <- bind_rows(
  procesar_bomba(data_raw, "I"),
  procesar_bomba(data_raw, "J"),
  procesar_bomba(data_raw, "K")
)

# Ventana temporal del evento (9 marzo 2025)
t_apagado_I <- ymd_hms("2025-03-09 08:45:00", tz = "UTC")
t_reingreso_I <- ymd_hms("2025-03-09 08:57:00", tz = "UTC")
t_falla_J <- ymd_hms("2025-03-09 11:23:21", tz = "UTC")

# Ventana extendida: 3 horas antes y 2.5 horas después
data_evento <- data_operativo %>%
  filter(
    DateTime >= t_apagado_I - minutes(180),
    DateTime <= t_reingreso_I + minutes(150)
  ) %>%
  mutate(
    periodo = case_when(
      DateTime < t_apagado_I ~ "Pre-evento",
      DateTime >= t_apagado_I & DateTime < t_reingreso_I ~ "Transiente",
      TRUE ~ "Post-evento"
    )
  )

# ==============================================================
# ESTIMACIÓN DE PRESIÓN PARA BOMBA K
# ==============================================================
# Como K descarga al mismo cabezal que I y J (con válvulas check),
# su presión de descarga debe ser muy similar a la presión del header

data_evento <- data_evento %>%
  group_by(DateTime) %>%
  mutate(
    # Presión promedio del header basada en I y J
    pressure_header = mean(pressure_psi[bomba %in% c("I", "J")], na.rm = TRUE),
    # Para bomba K, usar presión del header si no tiene dato directo
    pressure_psi = if_else(
      bomba == "K" & is.na(pressure_psi),
      pressure_header,
      pressure_psi
    )
  ) %>%
  ungroup() %>%
  select(-pressure_header)

# Verificar datos por bomba
cat("\n=== REGISTROS POR BOMBA EN VENTANA DE ANÁLISIS ===\n")
## 
## === REGISTROS POR BOMBA EN VENTANA DE ANÁLISIS ===
resumen_bombas <- data_evento %>%
  group_by(bomba, periodo) %>%
  summarise(
    n_total = n(),
    n_con_presion = sum(!is.na(pressure_psi)),
    n_con_flujo = sum(!is.na(flow_bpd) & flow_bpd > 0),
    .groups = "drop"
  )
print(resumen_bombas)
## # A tibble: 9 × 5
##   bomba periodo     n_total n_con_presion n_con_flujo
##   <chr> <chr>         <int>         <int>       <int>
## 1 I     Post-evento     151           151         151
## 2 I     Pre-evento      180           180         180
## 3 I     Transiente       12            12          12
## 4 J     Post-evento     151           151         151
## 5 J     Pre-evento      180           180         180
## 6 J     Transiente       12            12          12
## 7 K     Post-evento     151           151         147
## 8 K     Pre-evento      180           180         180
## 9 K     Transiente       12            12          12

3. Ajuste del modelo parabólico

3.1 Método de resolución

Para cada frecuencia, resolvemos el sistema lineal de 3 ecuaciones:

\[ \begin{bmatrix} 1 & -Q_{min} & -Q_{min}^2 \\ 1 & -Q_{BEP} & -Q_{BEP}^2 \\ 1 & -Q_{max} & -Q_{max}^2 \end{bmatrix} \begin{bmatrix} H_0 \\ A \\ B \end{bmatrix} = \begin{bmatrix} H_{min} \\ H_{BEP} \\ H_{max} \end{bmatrix} \]

ajustar_curva <- function(Q_min, H_min, Q_bep, H_bep, Q_max, H_max) {
  M <- matrix(c(
    1, -Q_min, -Q_min^2,
    1, -Q_bep, -Q_bep^2,
    1, -Q_max, -Q_max^2
  ), nrow = 3, byrow = TRUE)
  
  b <- c(H_min, H_bep, H_max)
  coef <- solve(M, b)
  
  list(H0 = coef[1], A = coef[2], B = coef[3])
}

# Aplicar a todas las frecuencias
curva_fabricante <- curva_fabricante %>%
  rowwise() %>%
  mutate(
    coef = list(ajustar_curva(
      Q_minflow, H_minflow,
      Q_design, H_target,
      Q_maxflow, H_maxflow
    ))
  ) %>%
  ungroup() %>%
  mutate(
    H0 = sapply(coef, function(x) x$H0),
    A = sapply(coef, function(x) x$A),
    B = sapply(coef, function(x) x$B)
  )

3.2 Coeficientes por frecuencia

curva_fabricante %>%
  select(hz, H0, A, B, Q_design, H_target) %>%
  filter(hz %in% seq(30, 60, 5)) %>%
  gt_rapids(
    title = "Coeficientes del Modelo Parabólico",
    subtitle = "H(Q) = H₀ - A·Q - B·Q²",
    footer = "**Método:** Resolución de sistema lineal 3×3 por frecuencia"
  ) %>%
  fmt_number(columns = H0, decimals = 1) %>%
  fmt_scientific(columns = c(A, B), decimals = 3) %>%
  fmt_number(columns = c(Q_design, H_target), decimals = 0) %>%
  cols_label(
    hz = "Frecuencia (Hz)",
    H0 = "H₀ (psi)",
    A = "A (psi/bpd)",
    B = "B (psi/bpd²)",
    Q_design = "Q BEP (bpd)",
    H_target = "H BEP (psi)"
  )
Coeficientes del Modelo Parabólico
H(Q) = H₀ - A·Q - B·Q²
Frecuencia (Hz) H₀ (psi) A (psi/bpd) B (psi/bpd²) Q BEP (bpd) H BEP (psi)
30 434.4 −4.468 × 10−3 1.377 × 10−7 32,086 436
40 768.8 −6.161 × 10−3 1.401 × 10−7 42,781 776
50 1,203.2 −7.609 × 10−3 1.392 × 10−7 53,476 1,212
60 1,727.9 −9.283 × 10−3 1.403 × 10−7 64,171 1,746
Método: Resolución de sistema lineal 3×3 por frecuencia

4. Generación de curvas continuas

generar_curva_HQ <- function(Hz, H0, A, B, Q_min, Q_max, n = 200) {
  Q_seq <- seq(Q_min * 0.7, Q_max * 1.2, length.out = n)
  H_seq <- pmax(H0 - A * Q_seq - B * Q_seq^2, 0)
  
  data.frame(hz = Hz, Q = Q_seq, H = H_seq)
}

curvas_completas <- lapply(1:nrow(curva_fabricante), function(i) {
  row <- curva_fabricante[i, ]
  generar_curva_HQ(row$hz, row$H0, row$A, row$B, row$Q_minflow, row$Q_maxflow)
}) %>% bind_rows()

# Estadísticas
cat(sprintf("Curvas generadas: %d frecuencias × %d puntos = %d observaciones\n",
            length(unique(curvas_completas$hz)),
            nrow(curvas_completas) / length(unique(curvas_completas$hz)),
            nrow(curvas_completas)))
## Curvas generadas: 16 frecuencias × 200 puntos = 3200 observaciones

4.1 Identificación de límites operativos

# MCSF (Minimum Continuous Stable Flow)
linea_mcsf <- curva_fabricante %>%
  transmute(hz, Q_mcsf = Q_minflow, H_mcsf = H_minflow)

# Runout (máximo recomendado)
linea_runout <- curva_fabricante %>%
  transmute(hz, Q_runout = Q_maxflow, H_runout = H_maxflow)

# BEP (diseño)
linea_bep <- curva_fabricante %>%
  transmute(hz, Q_bep = Q_design, H_bep = H_target)

# BEP de referencia (52 Hz)
Hz_ref <- 52
BEP_ref <- curva_fabricante %>% filter(hz == Hz_ref)
Q_BEP <- BEP_ref$Q_design
H_BEP <- BEP_ref$H_target

5. Visualización multifrecuencia

# Frecuencias destacadas
freq_dest <- c(30, 36, 40, 46, 50, 52, 54, 56, 60)
curvas_plot <- curvas_completas %>% filter(hz %in% freq_dest)

# Puntos operativos promedio (ahora incluye K con presión estimada)
puntos_op <- data_evento %>%
  filter(
    periodo %in% c("Pre-evento", "Post-evento"),
    flow_bpd > 0, 
    speed_hz > 0,
    !is.na(pressure_psi)
  ) %>%
  group_by(bomba, periodo) %>%
  summarise(
    Q_op = mean(flow_bpd, na.rm = TRUE),
    H_op = mean(pressure_psi, na.rm = TRUE),
    N_op = mean(speed_hz, na.rm = TRUE),
    P_op = mean(power_kw, na.rm = TRUE),  # Agregar potencia
    n_obs = n(),
    .groups = "drop"
  )

# Verificar que tenemos las 3 bombas
cat("\n=== PUNTOS OPERATIVOS CALCULADOS ===\n")
## 
## === PUNTOS OPERATIVOS CALCULADOS ===
print(puntos_op)
## # A tibble: 6 × 7
##   bomba periodo       Q_op  H_op  N_op  P_op n_obs
##   <chr> <chr>        <dbl> <dbl> <dbl> <dbl> <int>
## 1 I     Post-evento 35761. 1437.  52.6  996.   148
## 2 I     Pre-evento  61964. 1151.  51.0 1192.   180
## 3 J     Post-evento 50873. 1429.  52.7 1225.   149
## 4 J     Pre-evento  66816. 1152.  51.0 1248.   180
## 5 K     Post-evento 52064. 1447.  53   1150.   147
## 6 K     Pre-evento  66440. 1152.  50.0 1110.   180
# Gráfica principal
g_curvas <- ggplot() +
  
  # Zona de recirculación
  annotate(
    "text",
    x = max(linea_mcsf$Q_mcsf)*1.05,
    y = max(linea_mcsf$H_mcsf)*1.05,
    label = "MSCF",
    color = "#dc2626",
    size = 4,
    fontface = "bold"
  ) +
  
  # Zona de runout
  annotate(
    "text",
    x = max(linea_runout$Q_runout)*1.05,
    y = max(linea_runout$H_runout)*1.05,
    label = "RUNOUT",
    color = "#dc2626",
    size = 4,
    fontface = "bold"
  ) +
  
  # Zona óptima (±10% BEP)
  annotate(
    "text",
    x = max(linea_bep$Q_bep)*1.05,
    y = max(linea_bep$H_bep)*1.05,
    label = "BEP",
    color = "#059669",
    size = 4,
    fontface = "bold"
  ) +
  
  # Curvas H-Q
  geom_line(
    data = curvas_plot,
    aes(x = Q, y = H, group = hz, color = as.factor(hz)),
    linewidth = 1.2,
    alpha = 0.85
  ) +
  
  # Puntos BEP
  geom_point(
    data = curva_fabricante %>% filter(hz %in% freq_dest),
    aes(x = Q_design, y = H_target, fill = as.factor(hz)),
    size = 3.5,
    shape = 21,
    color = "white",
    stroke = 1.2
  ) +
  
  # Límites MCSF y Runout
  geom_line(
    data = linea_mcsf %>% filter(hz %in% freq_dest),
    aes(x = Q_mcsf, y = H_mcsf),
    linetype = "dashed",
    color = "#dc2626",
    linewidth = 1.3
  ) +
  geom_line(
    data = linea_runout %>% filter(hz %in% freq_dest),
    aes(x = Q_runout, y = H_runout),
    linetype = "dashed",
    color = "#f59e0b",
    linewidth = 1.3
  ) +
    geom_line(
    data = linea_bep %>% filter(hz %in% freq_dest),
    aes(x = Q_bep, y = H_bep),
    color = "#059669",
    linewidth = 50,
    alpha = 0.25

  ) +
  
  # PUNTOS OPERATIVOS REALES - LAS 3 BOMBAS
  geom_point(
    data = puntos_op,
    aes(
      x = Q_op, 
      y = H_op, 
      color = bomba,
      shape = periodo
    ),
    size = 6,
    stroke = 2
  ) +
  
  # Etiquetas identificando bomba (sin periodo en el label)
  geom_text(
    data = puntos_op,
    aes(
      x = Q_op, 
      y = H_op, 
      label = bomba,
      color = bomba
    ),
    size = 3.5,
    fontface = "bold",
    vjust = -2,
    show.legend = FALSE
  ) +
  
  # BEP destacado
  geom_point(
    x = Q_BEP, y = H_BEP,
    size = 6,
    color = "#059669",
    shape = 18
  ) +
  annotate(
    "label",
    x = Q_BEP * 1.15, y = H_BEP,
    label = sprintf("BEP @ %d Hz\n%s bpd | %s psi",
                    Hz_ref,
                    format(round(Q_BEP), big.mark = ","),
                    format(round(H_BEP), big.mark = ",")),
    color = "#059669",
    fontface = "bold",
    size = 3.5,
    fill = "white",
    label.padding = unit(0.3, "lines")
  ) +
  
  # Etiquetas de frecuencia
  geom_text(
    data = curva_fabricante %>% filter(hz %in% seq(30, 60, 10)),
    aes(x = Q_maxflow * 1.05, y = H_maxflow, 
        label = sprintf("%d Hz", hz),
        color = as.factor(hz)),
    size = 3.5,
    fontface = "bold",
    hjust = 0,
    show.legend = FALSE
  ) +
  
  scale_color_manual(
    values = c(
      "I" = colores_rapids["azul_claro"],
      "J" = colores_rapids["rojo"],
      "K" = colores_rapids["morado"],
      # Colores para frecuencias
      setNames(viridis::plasma(length(freq_dest)), as.character(freq_dest))
    ),
    breaks = c("I", "J", "K"),
    name = "Bomba"
  ) +
  scale_fill_viridis_d(option = "plasma", guide = "none") +
  scale_shape_manual(
    values = c("Pre-evento" = 16, "Post-evento" = 17),
    name = "Periodo"
  ) +
  scale_x_continuous(
    labels = label_number(big.mark = ",", scale = 1e-3, suffix = "k"),
    breaks = seq(0, 120000, 10000)
  ) +
  scale_y_continuous(
    labels = label_comma(),
    breaks = seq(0, 2000, 200)
  ) +
  coord_cartesian(xlim = c(0, 115000), ylim = c(0, 2000)) +
  labs(
    title = "Curvas H-Q Multifrecuencia | Summit ESP ST2500CCW",
    subtitle = "Modelo analítico parabólico validado con datos operativos | Evento 9 marzo 2025 | Bombas I, J, K",
    x = "Caudal (bpd)",
    y = "Presión diferencial (psi)",
    caption = "JR Engineering | Fuente: Halliburton + SCADA CPE6 | Presión K estimada desde header común"
  ) +
  theme_rapids() +
  theme(
    legend.position = "bottom",
    legend.background = element_rect(fill = "white", color = "gray80"),
    legend.key.size = unit(0.5, "cm")
  )

print(g_curvas)

6. Validación del modelo

6.1 Comparación puntual @ 52 Hz

curva_52 <- curvas_completas %>% filter(hz == Hz_ref)
puntos_52 <- curva_fabricante %>%
  filter(hz == Hz_ref) %>%
  select(Q_minflow, H_minflow, Q_design, H_target, Q_maxflow, H_maxflow) %>%
  pivot_longer(cols = everything(), names_to = c(".value", "tipo"), names_pattern = "(.*)_(.*)")

g_val <- ggplot() +
  geom_line(
    data = curva_52,
    aes(x = Q/1000, y = H),
    color = colores_rapids["azul_oscuro"],
    linewidth = 1.8
  ) +
  geom_point(
    data = puntos_52,
    aes(x = Q/1000, y = H, color = tipo),
    size = 6,
    shape = 21,
    fill = "white",
    stroke = 2.5
  ) +
  scale_color_manual(
    values = c(
      "minflow" = "#dc2626",
      "design" = "#059669",
      "maxflow" = "#f59e0b"
    ),
    labels = c("Mínimo continuo", "BEP", "Máximo"),
    name = NULL
  ) +
  labs(
    title = "Validación del Modelo Analítico",
    subtitle = sprintf("Curva @ %d Hz | Ajuste perfecto en 3 puntos de diseño", Hz_ref),
    x = "Caudal (kbpd)",
    y = "Presión (psi)"
  ) +
  theme_rapids()

print(g_val)

6.2 Análisis de eficiencia energética

6.2.1 Eficiencia específica del sistema

# ==============================================================
# ANÁLISIS DE EFICIENCIA OPERACIONAL: kW/KBPD
# ==============================================================

eficiencia_sistema <- puntos_op %>%
  group_by(periodo) %>%
  summarise(
    Q_total = sum(Q_op, na.rm = TRUE),
    P_total = sum(P_op, na.rm = TRUE),  # Cambiar power_kw por P_op
    n_bombas = n(),
    .groups = "drop"
  ) %>%
  mutate(
    Q_total_kbpd = Q_total / 1000,
    eficiencia_especifica = P_total / Q_total_kbpd,
    productividad = Q_total / n_bombas
  )

# Calcular variación
delta_eficiencia <- eficiencia_sistema %>%
  select(periodo, eficiencia_especifica, Q_total_kbpd, P_total) %>%
  pivot_wider(
    names_from = periodo,
    values_from = c(eficiencia_especifica, Q_total_kbpd, P_total)
  ) %>%
  mutate(
    delta_eficiencia_pct = 100 * (`eficiencia_especifica_Post-evento` - 
                                   `eficiencia_especifica_Pre-evento`) / 
                                   `eficiencia_especifica_Pre-evento`,
    delta_Q_pct = 100 * (`Q_total_kbpd_Post-evento` - `Q_total_kbpd_Pre-evento`) / 
                         `Q_total_kbpd_Pre-evento`,
    delta_P_pct = 100 * (`P_total_Post-evento` - `P_total_Pre-evento`) / 
                         `P_total_Pre-evento`
  )

# Tabla resumen del sistema
eficiencia_sistema %>%
  gt_rapids(
    title = "Eficiencia Energética del Sistema",
    subtitle = "Consumo específico por producción | Comparativo pre/post evento",
    footer = "**Métrica clave:** kW/KBPD indica energía requerida para mover 1000 barriles por día"
  ) %>%
  fmt_number(columns = c(Q_total, P_total, productividad), decimals = 0) %>%
  fmt_number(columns = Q_total_kbpd, decimals = 1) %>%
  fmt_number(columns = eficiencia_especifica, decimals = 2) %>%
  cols_label(
    periodo = "Periodo",
    Q_total = "Q total (bpd)",
    Q_total_kbpd = "Q total (kbpd)",
    P_total = "Potencia total (kW)",
    n_bombas = "Bombas activas",
    eficiencia_especifica = "kW/KBPD",
    productividad = "Q por bomba (bpd)"
  ) %>%
  tab_style(
    style = cell_fill(color = "#d1fae5"),
    locations = cells_body(
      columns = eficiencia_especifica,
      rows = eficiencia_especifica == min(eficiencia_especifica)
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#fee2e2"),
    locations = cells_body(
      columns = eficiencia_especifica,
      rows = eficiencia_especifica == max(eficiencia_especifica)
    )
  )
Eficiencia Energética del Sistema
Consumo específico por producción | Comparativo pre/post evento
Periodo Q total (bpd) Potencia total (kW) Bombas activas Q total (kbpd) kW/KBPD Q por bomba (bpd)
Post-evento 138,698 3,371 3 138.7 24.30 46,233
Pre-evento 195,220 3,550 3 195.2 18.19 65,073
Métrica clave: kW/KBPD indica energía requerida para mover 1000 barriles por día 
# Mostrar variación
cat("\n=== VARIACIÓN POST-EVENTO ===\n")
## 
## === VARIACIÓN POST-EVENTO ===
cat(sprintf("Eficiencia específica: %+.2f%% (%s)\n", 
            delta_eficiencia$delta_eficiencia_pct,
            ifelse(delta_eficiencia$delta_eficiencia_pct > 0, "PEOR", "MEJOR")))
## Eficiencia específica: +33.65% (PEOR)
cat(sprintf("Producción total: %+.1f%%\n", delta_eficiencia$delta_Q_pct))
## Producción total: -29.0%
cat(sprintf("Consumo total: %+.1f%%\n", delta_eficiencia$delta_P_pct))
## Consumo total: -5.0%

6.2.2 Eficiencia por bomba

# Eficiencia individual de cada bomba
eficiencia_bombas <- puntos_op %>%
  mutate(
    Q_kbpd = Q_op / 1000,
    eficiencia_especifica = P_op / Q_kbpd,  # Cambiar power_kw por P_op
    Q_gpm = Q_op / 1440 / 0.159, # bpd -> gpm
    Hp_hyd = Q_gpm * H_op / 1714,  # Potencia hidráulica en HP: psi + gpm -> HP
    eficiencia_global_pct = 100 * (Hp_hyd * 0.746) / P_op  # P_op: Potencia de operación que ya está en kW
  )

# Tabla por bomba
eficiencia_bombas %>%
  select(bomba, periodo, N_op, Q_op, H_op, P_op, eficiencia_especifica, eficiencia_global_pct) %>%
  arrange(bomba, periodo) %>%
  gt_rapids(
    title = "Eficiencia Energética Individual",
    subtitle = "Desempeño por bomba antes y después del evento | Identificación de ineficiencias",
    footer = "**Eficiencia global:** Relación entre potencia hidráulica útil y potencia eléctrica consumida"
  ) %>%
  fmt_number(columns = c(N_op, Q_op, H_op, P_op), decimals = 0) %>%
  fmt_number(columns = eficiencia_especifica, decimals = 2) %>%
  fmt_number(columns = eficiencia_global_pct, decimals = 1) %>%
  cols_label(
    bomba = "Bomba",
    periodo = "Periodo",
    N_op = "Velocidad (Hz)",
    Q_op = "Caudal (bpd)",
    H_op = "Presión (psi)",
    P_op = "Potencia (kW)",
    eficiencia_especifica = "kW/KBPD",
    eficiencia_global_pct = "Eficiencia (%)"
  ) %>%
  tab_style(
    style = cell_fill(color = "#fee2e2"),
    locations = cells_body(
      columns = eficiencia_especifica,
      rows = eficiencia_especifica > quantile(eficiencia_especifica, 0.75)
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#d1fae5"),
    locations = cells_body(
      columns = eficiencia_especifica,
      rows = eficiencia_especifica < quantile(eficiencia_especifica, 0.25)
    )
  )
Eficiencia Energética Individual
Desempeño por bomba antes y después del evento | Identificación de ineficiencias
Bomba Periodo Velocidad (Hz) Caudal (bpd) Presión (psi) Potencia (kW) kW/KBPD Eficiencia (%)
I Post-evento 53 35,761 1,437 996 27.84 9.8
I Pre-evento 51 61,964 1,151 1,192 19.24 11.4
J Post-evento 53 50,873 1,429 1,225 24.09 11.3
J Pre-evento 51 66,816 1,152 1,248 18.67 11.7
K Post-evento 53 52,064 1,447 1,150 22.09 12.4
K Pre-evento 50 66,440 1,152 1,110 16.71 13.1
Eficiencia global: Relación entre potencia hidráulica útil y potencia eléctrica consumida 
# Análisis de variación por bomba
variacion_bombas <- eficiencia_bombas %>%
  select(bomba, periodo, eficiencia_especifica, eficiencia_global_pct) %>%
  pivot_wider(
    names_from = periodo,
    values_from = c(eficiencia_especifica, eficiencia_global_pct)
  ) %>%
  mutate(
    delta_esp_pct = 100 * (`eficiencia_especifica_Post-evento` - 
                           `eficiencia_especifica_Pre-evento`) / 
                           `eficiencia_especifica_Pre-evento`,
    delta_glob_pct = 100 * (`eficiencia_global_pct_Post-evento` - 
                            `eficiencia_global_pct_Pre-evento`) / 
                            `eficiencia_global_pct_Pre-evento`
  )

variacion_bombas %>%
  gt_rapids(
    title = "Variación de Eficiencia Post-Evento",
    subtitle = "Cambio porcentual en desempeño energético por bomba",
    footer = "**Valores positivos** en Δ Específica (%) indican deterioro (mayor consumo energético por unidad de caudal) | **Valores negativos** en Δ Global (%) indican deterioro, reflejando una menor conversión de energía eléctrica en trabajo hidráulico útil"
  ) %>%
  fmt_number(columns = starts_with("eficiencia"), decimals = 2) %>%
  fmt_number(columns = starts_with("delta"), decimals = 1) %>%
  cols_label(
    bomba = "Bomba",
    `eficiencia_especifica_Pre-evento` = "Pre (kW/KBPD)",
    `eficiencia_especifica_Post-evento` = "Post (kW/KBPD)",
    delta_esp_pct = "Δ Específica (%)",
    `eficiencia_global_pct_Pre-evento` = "Pre (%)",
    `eficiencia_global_pct_Post-evento` = "Post (%)",
    delta_glob_pct = "Δ Global (%)"
  ) %>%
  tab_spanner(
    label = "Eficiencia Específica",
    columns = c(`eficiencia_especifica_Pre-evento`, 
                `eficiencia_especifica_Post-evento`, 
                delta_esp_pct)
  ) %>%
  tab_spanner(
    label = "Eficiencia Global",
    columns = c(`eficiencia_global_pct_Pre-evento`, 
                `eficiencia_global_pct_Post-evento`, 
                delta_glob_pct)
  ) %>%
  tab_style(
    style = cell_fill(color = "#fee2e2"),
    locations = cells_body(
      columns = delta_esp_pct,
      rows = delta_esp_pct > 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#d1fae5"),
    locations = cells_body(
      columns = delta_esp_pct,
      rows = delta_esp_pct < 0
    )
  )
Variación de Eficiencia Post-Evento
Cambio porcentual en desempeño energético por bomba
Bomba
Eficiencia Específica
Eficiencia Global
Pre (kW/KBPD) Post (kW/KBPD) Δ Específica (%) Pre (%) Post (%) Δ Global (%)
I 19.24 27.84 44.7 11.37 9.81 −13.7
J 18.67 24.09 29.0 11.73 11.28 −3.9
K 16.71 22.09 32.2 13.10 12.45 −5.0
Valores positivos en Δ Específica (%) indican deterioro (mayor consumo energético por unidad de caudal) | Valores negativos en Δ Global (%) indican deterioro, reflejando una menor conversión de energía eléctrica en trabajo hidráulico útil

6.2.3 Visualización comparativa

# Gráfica comparativa de eficiencia
g_ef1 <- ggplot(eficiencia_bombas, aes(x = bomba, y = eficiencia_especifica, fill = periodo)) +
  geom_col(position = "dodge", width = 0.7) +
  geom_text(
    aes(label = sprintf("%.2f", eficiencia_especifica)),
    position = position_dodge(width = 0.7),
    vjust = -0.5,
    size = 3.5,
    fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("Pre-evento" = colores_rapids["azul_claro"], 
               "Post-evento" = colores_rapids["rojo"])
  ) +
  labs(
    title = "Consumo Específico por Bomba",
    subtitle = "kW/KBPD | Menor es mejor",
    x = NULL,
    y = "kW / KBPD",
    fill = "Periodo"
  ) +
  theme_rapids()

g_ef2 <- ggplot(eficiencia_bombas, aes(x = bomba, y = eficiencia_global_pct, fill = periodo)) +
  geom_col(position = "dodge", width = 0.7) +
  geom_text(
    aes(label = sprintf("%.1f%%", eficiencia_global_pct)),
    position = position_dodge(width = 0.7),
    vjust = -0.5,
    size = 3.5,
    fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("Pre-evento" = colores_rapids["azul_claro"], 
               "Post-evento" = colores_rapids["rojo"])
  ) +
  labs(
    title = "Eficiencia Global por Bomba",
    subtitle = "Potencia hidráulica / Potencia eléctrica | Mayor es mejor",
    x = NULL,
    y = "Eficiencia (%)",
    fill = "Periodo"
  ) +
  theme_rapids()

g_ef1 / g_ef2

6.3 Análisis de posición operativa

# Calcular BEP equivalente por afinidad
puntos_op_ext <- puntos_op %>%
  mutate(
    Q_BEP_eq = Q_BEP * (N_op / Hz_ref),
    H_BEP_eq = H_BEP * (N_op / Hz_ref)^2,
    dQ_pct = 100 * (Q_op - Q_BEP_eq) / Q_BEP_eq,
    dH_pct = 100 * (H_op - H_BEP_eq) / H_BEP_eq,
    posicion = case_when(
      abs(dQ_pct) < 10 ~ "Cerca del BEP",
      dQ_pct < -10 ~ "Izquierda BEP (recirculación)",
      dQ_pct > 10 ~ "Derecha BEP (runout)"
    )
  )

puntos_op_ext %>%
  select(bomba, periodo, N_op, Q_op, H_op, Q_BEP_eq, dQ_pct, posicion) %>%
  gt_rapids(
    title = "Posición Operativa vs BEP Teórico",
    subtitle = "Desviación porcentual respecto al punto de máxima eficiencia",
    footer = "**Nota:** BEP equivalente calculado por leyes de afinidad desde 52 Hz nominal | **Bomba I post-evento** en zona crítica cercana a MCSF"
  ) %>%
  fmt_number(columns = c(N_op, Q_op, H_op, Q_BEP_eq), decimals = 0) %>%
  fmt_number(columns = dQ_pct, decimals = 1) %>%
  cols_label(
    bomba = "Bomba",
    periodo = "Periodo",
    N_op = "N (Hz)",
    Q_op = "Q real (bpd)",
    H_op = "H real (psi)",
    Q_BEP_eq = "Q BEP equiv (bpd)",
    dQ_pct = "Δ Q (%)",
    posicion = "Posición"
  ) %>%
  tab_style(
    style = cell_fill(color = "#fee2e2"),
    locations = cells_body(
      columns = posicion,
      rows = str_detect(posicion, "recirculación")
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#d1fae5"),
    locations = cells_body(
      columns = posicion,
      rows = str_detect(posicion, "Cerca")
    )
  )
Posición Operativa vs BEP Teórico
Desviación porcentual respecto al punto de máxima eficiencia
Bomba Periodo N (Hz) Q real (bpd) H real (psi) Q BEP equiv (bpd) Δ Q (%) Posición
I Post-evento 53 35,761 1,437 56,287 −36.5 Izquierda BEP (recirculación)
I Pre-evento 51 61,964 1,151 54,563 13.6 Derecha BEP (runout)
J Post-evento 53 50,873 1,429 56,333 −9.7 Cerca del BEP
J Pre-evento 51 66,816 1,152 54,557 22.5 Derecha BEP (runout)
K Post-evento 53 52,064 1,447 56,685 −8.2 Cerca del BEP
K Pre-evento 50 66,440 1,152 53,488 24.2 Derecha BEP (runout)
Nota: BEP equivalente calculado por leyes de afinidad desde 52 Hz nominal | Bomba I post-evento en zona crítica cercana a MCSF

6.3.1 Diagnóstico de ineficiencias

Hallazgos principales del análisis energético:

# Generar diagnóstico automático
bomba_mas_ineficiente <- variacion_bombas %>%
  filter(delta_esp_pct == max(delta_esp_pct))

sistema_delta <- delta_eficiencia$delta_eficiencia_pct

cat(sprintf("1. **Impacto sistémico:** El sistema experimentó un cambio de **%+.1f%%** en consumo específico post-evento (%s eficiencia)\n\n",
            sistema_delta,
            ifelse(sistema_delta > 0, "PÉRDIDA de", "GANANCIA de")))
  1. Impacto sistémico: El sistema experimentó un cambio de +33.6% en consumo específico post-evento (PÉRDIDA de eficiencia)
cat(sprintf("2. **Bomba crítica:** La bomba **%s** muestra la mayor degradación con **%+.1f%%** en kW/KBPD\n\n",
            bomba_mas_ineficiente$bomba,
            bomba_mas_ineficiente$delta_esp_pct))
  1. Bomba crítica: La bomba I muestra la mayor degradación con +44.7% en kW/KBPD
cat("3. **Causa raíz identificada:** Bomba I reingresó al sistema con presión de header cercana a 1,500 psi, forzándola a operar:\n")
  1. Causa raíz identificada: Bomba I reingresó al sistema con presión de header cercana a 1,500 psi, forzándola a operar:
cat("   - Cerca del límite MCSF (zona de recirculación)\n")
  • Cerca del límite MCSF (zona de recirculación)
cat("   - Fuera de su zona de máxima eficiencia\n")
  • Fuera de su zona de máxima eficiencia
cat("   - Con desbalance de carga respecto a J y K\n\n")
  • Con desbalance de carga respecto a J y K
cat("4. **Estrategia correctiva recomendada:**\n")
  1. Estrategia correctiva recomendada:
cat("   - Reducción gradual de Hz en J y K durante reingreso de I\n")
  • Reducción gradual de Hz en J y K durante reingreso de I
cat("   - Permitir que header baje a ~1,300 psi antes de reingreso\n")
  • Permitir que header baje a ~1,300 psi antes de reingreso
cat("   - Sincronizar setpoints para distribución equitativa de carga\n")
  • Sincronizar setpoints para distribución equitativa de carga

7. Función predictiva generalizada

# ==============================================================
# FUNCIÓN: PREDECIR H PARA CUALQUIER (Q, Hz)
# ==============================================================

predecir_H <- function(Q, Hz) {
  # Interpolar coeficientes
  H0_int <- approx(curva_fabricante$hz, curva_fabricante$H0, xout = Hz, rule = 2)$y
  A_int <- approx(curva_fabricante$hz, curva_fabricante$A, xout = Hz, rule = 2)$y
  B_int <- approx(curva_fabricante$hz, curva_fabricante$B, xout = Hz, rule = 2)$y
  
  # Calcular presión
  H <- H0_int - A_int * Q - B_int * Q^2
  return(pmax(H, 0))
}

# Casos de ejemplo
casos_test <- data.frame(
  Q = c(45000, 55000, 60000, 40000),
  Hz = c(48, 52, 53.5, 45)
)

casos_test %>%
  mutate(H_pred = predecir_H(Q, Hz)) %>%
  gt_rapids(
    title = "Ejemplos de Predicción",
    subtitle = "Cálculo de presión para combinaciones arbitrarias Q-Hz",
    footer = "**Interpolación:** Lineal entre coeficientes de frecuencias adyacentes"
  ) %>%
  fmt_number(columns = c(Q, Hz, H_pred), decimals = 1) %>%
  cols_label(Q = "Caudal (bpd)", Hz = "Frecuencia (Hz)", H_pred = "H predicha (psi)")
Ejemplos de Predicción
Cálculo de presión para combinaciones arbitrarias Q-Hz
Caudal (bpd) Frecuencia (Hz) H predicha (psi)
45,000.0 48.0 1,155.7
55,000.0 52.0 1,315.6
60,000.0 53.5 1,365.5
40,000.0 45.0 1,026.5
Interpolación: Lineal entre coeficientes de frecuencias adyacentes

8. Conclusiones

Validación del modelo

El modelo parabólico de tres parámetros reproduce con precisión absoluta los puntos de diseño del fabricante en el rango 30-60 Hz. La continuidad hidráulica de las curvas generadas es consistente con la física de bombas centrífugas.

Ecuación maestra @ 52 Hz

\(H(Q) = 1303.60 - -0.007856 \cdot Q - 0.0000001389 \cdot Q^2\)

Con \(Q\) en bpd y \(H\) en psi. Esta ecuación permite calcular la presión esperada para cualquier caudal operativo.

Hallazgos operativos

  1. Bomba J (pre-falla): Operación a 22.5% del BEP equivalente, dentro de zona segura.

  2. Bombas I/K: Desplazamiento hacia izquierda del BEP durante transiente, indicando potencial recirculación temporal.

  3. Ventana predictiva: El modelo permite detectar desviaciones >10% respecto a curva teórica, señal temprana de degradación mecánica.

Aplicaciones prácticas

  • Monitoreo continuo: Comparar presión real vs predicha para detección de anomalías
  • Optimización de setpoints: Minimizar distancia al BEP según demanda
  • Planificación de mantenimiento: Identificar operación fuera de zona segura
  • Diagnóstico de fallas: Correlacionar desviaciones con análisis vibracional

Integración con análisis vibracional

Los precursores espectrales identificados en el análisis vibracional (bandas discretas 0.2-0.5 Hz) ocurrieron cuando la bomba J operaba dentro de su zona segura según el modelo hidráulico, confirmando que la falla fue de origen mecánico interno y no por condiciones hidráulicas adversas.

9. Referencias

  1. API Standard 610 (2010). Centrifugal Pumps for Petroleum, Petrochemical and Natural Gas Industries. American Petroleum Institute.

  2. Hydraulic Institute Standards (2014). ANSI/HI 9.6.3 - Rotodynamic Pumps: Guideline for Allowable Operating Region.

  3. Gülich, J.F. (2020). Centrifugal Pumps, 4th Edition. Springer International Publishing.

  4. Bloch, H.P. & Budris, A.R. (2014). Pump User’s Handbook: Life Extension, 4th Edition. Fairmont Press.

  5. ISO 20816-3 (2022). Mechanical vibration - Measurement and evaluation of machine vibration - Part 3: Industrial machines.

  6. Halliburton Energy Services (2023). HPS System Installation & Pre-commissioning Manual. Documento interno CPE6-360K.

Apéndices

A. Tabla completa de coeficientes

curva_fabricante %>%
  select(hz, H0, A, B, Q_design, H_target, Q_minflow, Q_maxflow) %>%
  gt_rapids(
    title = "Coeficientes Completos del Modelo",
    subtitle = "16 frecuencias operativas | 30-60 Hz",
    footer = "**Nota:** Coeficientes calculados mediante resolución de sistema lineal 3×3"
  ) %>%
  fmt_number(columns = c(H0, H_target), decimals = 1) %>%
  fmt_scientific(columns = c(A, B), decimals = 4) %>%
  fmt_number(columns = c(Q_design, Q_minflow, Q_maxflow), decimals = 0) %>%
  cols_label(
    hz = "Hz",
    H0 = "H₀ (psi)",
    A = "A (psi/bpd)",
    B = "B (psi/bpd²)",
    Q_design = "Q BEP (bpd)",
    H_target = "H BEP (psi)",
    Q_minflow = "Q mín (bpd)",
    Q_maxflow = "Q máx (bpd)"
  )
Coeficientes Completos del Modelo
16 frecuencias operativas | 30-60 Hz
Hz H₀ (psi) A (psi/bpd) B (psi/bpd²) Q BEP (bpd) H BEP (psi) Q mín (bpd) Q máx (bpd)
30 434.4 −4.4684 × 10−3 1.3771 × 10−7 32,086 436.0 19,681 53,530
32 491.7 −4.9627 × 10−3 1.4051 × 10−7 34,225 497.0 20,993 57,098
34 552.1 −5.3982 × 10−3 1.4172 × 10−7 36,364 561.0 22,305 60,667
36 622.9 −5.5650 × 10−3 1.4041 × 10−7 38,503 629.0 23,617 64,235
38 696.7 −5.7136 × 10−3 1.3858 × 10−7 40,642 700.0 24,929 67,804
40 768.8 −6.1607 × 10−3 1.4006 × 10−7 42,781 776.0 26,241 71,373
42 844.0 −6.6169 × 10−3 1.4136 × 10−7 44,920 856.0 27,553 74,941
44 929.5 −6.8014 × 10−3 1.4026 × 10−7 47,059 939.0 28,865 78,510
46 1,018.2 −7.0164 × 10−3 1.3940 × 10−7 49,198 1,026.0 30,177 82,079
48 1,109.1 −7.3034 × 10−3 1.3927 × 10−7 51,337 1,117.0 31,489 85,647
50 1,203.2 −7.6093 × 10−3 1.3921 × 10−7 53,476 1,212.0 32,802 89,216
52 1,303.6 −7.8564 × 10−3 1.3887 × 10−7 55,615 1,311.0 34,114 92,785
54 1,402.2 −8.2787 × 10−3 1.3979 × 10−7 57,754 1,414.0 35,425 96,353
56 1,508.0 −8.5989 × 10−3 1.3994 × 10−7 59,893 1,521.0 36,737 99,922
58 1,616.8 −8.8969 × 10−3 1.3973 × 10−7 62,032 1,631.0 38,050 103,491
60 1,727.9 −9.2830 × 10−3 1.4027 × 10−7 64,171 1,746.0 39,361 107,059
Nota: Coeficientes calculados mediante resolución de sistema lineal 3×3

B. Código de la función predictiva

# Función optimizada para producción
predecir_operacion <- function(Q_bpd, Hz, 
                                curva_ref = curva_fabricante,
                                incluir_limites = TRUE) {
  
  # Interpolar coeficientes con extrapolación lineal
  H0 <- approx(curva_ref$hz, curva_ref$H0, xout = Hz, rule = 2)$y
  A <- approx(curva_ref$hz, curva_ref$A, xout = Hz, rule = 2)$y
  B <- approx(curva_ref$hz, curva_ref$B, xout = Hz, rule = 2)$y
  
  # Calcular presión predicha
  H_pred <- pmax(H0 - A * Q_bpd - B * Q_bpd^2, 0)
  
  resultado <- data.frame(
    Q = Q_bpd,
    Hz = Hz,
    H_pred = H_pred
  )
  
  if (incluir_limites) {
    # BEP equivalente por afinidad
    BEP_52 <- curva_ref %>% filter(hz == 52)
    Q_BEP_eq <- BEP_52$Q_design * (Hz / 52)
    H_BEP_eq <- BEP_52$H_target * (Hz / 52)^2
    
    # MCSF y Runout
    Q_mcsf <- approx(curva_ref$hz, curva_ref$Q_minflow, xout = Hz, rule = 2)$y
    Q_runout <- approx(curva_ref$hz, curva_ref$Q_maxflow, xout = Hz, rule = 2)$y
    
    resultado <- resultado %>%
      mutate(
        Q_BEP_eq = Q_BEP_eq,
        H_BEP_eq = H_BEP_eq,
        dQ_BEP_pct = 100 * (Q - Q_BEP_eq) / Q_BEP_eq,
        zona = case_when(
          Q < Q_mcsf ~ "Recirculación",
          Q > Q_runout ~ "Runout",
          abs(Q - Q_BEP_eq) / Q_BEP_eq < 0.10 ~ "Óptima",
          TRUE ~ "Aceptable"
        )
      )
  }
  
  return(resultado)
}

# Ejemplo de uso en monitoreo continuo
data_monitoreo <- data.frame(
  timestamp = seq(Sys.time(), by = "1 min", length.out = 10),
  Q_medido = rnorm(10, 55000, 2000),
  Hz_setpoint = 52
)

data_monitoreo %>%
  rowwise() %>%
  mutate(
    prediccion = list(predecir_operacion(Q_medido, Hz_setpoint))
  ) %>%
  unnest(prediccion) %>%
  mutate(
    desviacion_pct = 100 * abs(H_pred - H_medido) / H_pred,
    alerta = desviacion_pct > 5
  )

C. Gráfica interactiva 3D

library(plotly)
library(dplyr)

# ==============================================================
# 1) Curvas ISO-Hz con clasificación final de zonas
# ==============================================================

curvas_zonas <- curvas_completas %>%
  left_join(
    curva_fabricante %>%
      select(hz, Q_design),
    by = "hz"
  ) %>%
  left_join(
    linea_mcsf %>% select(hz, Q_mcsf),
    by = "hz"
  ) %>%
  left_join(
    linea_runout %>% select(hz, Q_runout),
    by = "hz"
  ) %>%
  mutate(
    phi = Q / Q_design,

    zona = case_when(
      Q < Q_mcsf | Q > Q_runout              ~ "PELIGRO",
      phi >= 0.85 & phi <= 1.15              ~ "SEGURA",
      TRUE                                   ~ "ALARMA"
    ),

    zona_color = case_when(
      zona == "SEGURA"  ~ "#00A34D",
      zona == "ALARMA"  ~ "#FFB81C",
      zona == "PELIGRO" ~ "#E30613"
    )
  )

# ==============================================================
# 2) Gráfica base
# ==============================================================

fig <- plot_ly()

# ==============================================================
# 3) Curvas ISO-frecuencia coloreadas por zona
# ==============================================================

for (hz_val in unique(curvas_zonas$hz)) {

  datos_hz <- curvas_zonas %>% filter(hz == hz_val)

  fig <- fig %>%
    add_trace(
      data = datos_hz,
      x = ~Q,
      y = ~hz,
      z = ~H,
      type = "scatter3d",
      mode = "lines",
      line = list(color = ~zona_color, width = 4),
      showlegend = FALSE,
      hoverinfo = "text",
      text = ~paste0(
        "Hz: ", round(hz, 1), "<br>",
        "Q: ", round(Q, 0), " bpd<br>",
        "H: ", round(H, 1), " psi<br>",
        "Q/Q_BEP: ", round(phi, 2), "<br>",
        "Zona: ", zona
      )
    )
}

# ==============================================================
# 4) Líneas OEM – MCSF y Runout
# ==============================================================

fig <- fig %>%
  add_trace(
    data = linea_mcsf,
    x = ~Q_mcsf,
    y = ~hz,
    z = ~H_mcsf,
    type = "scatter3d",
    mode = "lines",
    line = list(color = "#dc2626", width = 4, dash = "dash"),
    name = "MCSF (OEM)"
  ) %>%
  add_trace(
    data = linea_runout,
    x = ~Q_runout,
    y = ~hz,
    z = ~H_runout,
    type = "scatter3d",
    mode = "lines",
    line = list(color = "#dc2626", width = 4, dash = "dash"),
    name = "Runout (OEM)"
  )

# ==============================================================
# 5) Línea BEP
# ==============================================================

fig <- fig %>%
  add_trace(
    data = curva_fabricante,
    x = ~Q_design,
    y = ~hz,
    z = ~H_target,
    type = "scatter3d",
    mode = "lines",
    line = list(color = "black", width = 6),
    name = "BEP"
  )

# ==============================================================
# 6) Puntos operativos reales
# ==============================================================

fig <- fig %>%
  add_trace(
    data = puntos_op,
    x = ~Q_op,
    y = ~N_op,
    z = ~H_op,
    type = "scatter3d",
    mode = "markers+text",
    marker = list(
      size = 7,
      color = ~ifelse(periodo == "Post-evento", "#000000", "#FFFFFF"),
      line = list(color = "#000000", width = 2)
    ),
    text = ~bomba,
    textposition = "top center",
    name = "Puntos operativos"
  )

# ==============================================================
# 7) Layout – leyenda abajo
# ==============================================================

fig <- fig %>%
  layout(
    title = list(
      text = "<b>Mapa 3D de Operación – Envolvente OEM y Zonas BEP</b><br>
              <sub>ISO-frecuencia ST2500CCW | BEP, MCSF y Runout</sub>"
    ),
    scene = list(
      xaxis = list(title = "Caudal (bpd)"),
      yaxis = list(title = "Frecuencia (Hz)"),
      zaxis = list(title = "Presión diferencial (psi)"),
      camera = list(eye = list(x = 1.6, y = -1.4, z = 1.3))
    ),
    legend = list(
      orientation = "h",
      x = 0.5,
      y = -0.15,
      xanchor = "center"
    ),
    margin = list(l = 0, r = 0, b = 100, t = 80)
  )

fig

Documento generado por JR Engineering
rapids analytics framework | Modelo hidráulico analítico v2.0
30 de December de 2025 - 05:00 -05

Este análisis integra datos operativos SCADA con especificaciones del fabricante para caracterización completa del comportamiento hidráulico de sistemas HPS en campo CPE6.