Estou com um codigo de calculo numerico para entrega sobre um sistema de aquecimento solar de uma piscina, porem o valor para o COP do sistema esta dando um valor absurdo, alem de que a temperatura para os graficos está da mesma forma. Creio que o erro está em alguma das equações ou então na variável Qa. Tenho 13 dados de temperatura para usar no código que também não foram usados. Alguém consegue me dar um norte para a questão?

import numpy as np
import matplotlib.pyplot as plt
# Dados do problema
# Definição dos parâmetros do sistema e do tempo de simulação
T0_pool = 25 # Temperatura inicial da piscina (°C)
T0_collector = 20 # Temperatura inicial da água no coletor (°C)
delta_T_on = 6 # Delta de temperatura para ligar a bomba (°C)
delta_T_off = 0 # Delta de temperatura para desligar a bomba (°C)
dt = 3600 # Passo de tempo em segundos (1 hora)
t0 = 6 * 3600 # Hora de início da simulação em segundos (6 horas)
tf = 19 * 3600 # Hora do fim da simulação em segundos (19 horas)

# Constantes físicas e características do sistema
rho_w = 998 # Densidade da água (kg/m^3)
cp_w = 4190 # Calor específico da água (J/kg·K)
V_pool = 15.85 # Volume da piscina (m^3)
K_capa = 0.05 # Condutividade térmica da capa (W/m·K)
delta_x = 0.00326 # Espessura da capa (m)
A_pool = 20.4 # Área da piscina coberta (m^2)
FR = 0.6 # Fator de remoção solar do coletor
U_L = 11.45 # Coeficiente global de perdas do coletor (W/m^2·K)
alpha_collector = 0.8 # Absorvidade do material do coletor
A_collector = 24 # Área do coletor solar (m^2)
V_water_collector = 0.12 # Volume de água no coletor (m^3)
v_bomb = 0.72 / 3600 # Vazão volumétrica da bomba (m^3/s)
Eb = 1.25 * 3600 # Consumo de energia da bomba (W·h/h)
# Inicialização das variáveis
T_pool = np.zeros(int((tf - t0) / dt) + 1)
T_collector = np.zeros(int((tf - t0) / dt) + 1)
Qa = np.zeros(int((tf - t0) / dt) + 1)
Qcd_capa = np.zeros(int((tf - t0) / dt) + 1)
# Simulação do sistema
T_pool[0] = T0_pool
T_collector[0] = T0_collector
# Loop de simulação para cada passo de tempo
for i in range(len(T_collector) - 1):
Tcd_capa = (T_pool[i] + T_collector[i]) / 2
Qcd_capa[i] = (K_capa * A_pool * (Tcd_capa - T0_pool)) / delta_x

if T_collector[i] - T_pool[i] > delta_T_on:
Qa[i] = v_bomb * rho_w * cp_w * (T_collector[i] - T_pool[i])
if v_bomb * V_water_collector < dt:
T_collector[i + 1] = (
(Qa[i] + FR * alpha_collector * rho_w * cp_w * (T_collector[i] - Tcd_capa) * A_collector)
/ (rho_w * cp_w * V_water_collector) + T_collector[i]
)
else:
T_collector[i + 1] = (
(Qa[i] + FR * alpha_collector * rho_w * cp_w * (T_collector[i] - Tcd_capa) * A_collector)
/ (rho_w * cp_w * v_bomb) + T_collector[i]
)
T_pool[i + 1] = (Qcd_capa[i] * dt + T_pool[i] * V_pool * rho_w * cp_w) / (V_pool * rho_w * cp_w)
else:
T_collector[i + 1] = (
(Qa[i] + FR * alpha_collector * rho_w * cp_w * (T_collector[i] - Tcd_capa) * A_collector)
/ (rho_w * cp_w * v_bomb) + T_collector[i]
)
T_pool[i + 1] = (Qcd_capa[i] * dt + T_pool[i] * V_pool * rho_w * cp_w) / (V_pool * rho_w * cp_w)
# Cálculo do COP
Qa_total = np.sum(Qa) * dt
COP = Qa_total / (Eb * (tf - t0) / 3600)
# Resultados
print("a) Temperatura da água da piscina no final do dia simulado:", T_pool[-1])
plt.figure(figsize=(12, 6))
plt.subplot(2, 1, 1)
plt.plot(np.arange(t0, tf + dt, dt) / 3600, T_pool, label="Temperatura da piscina")
plt.plot(np.arange(t0, tf + dt, dt) / 3600, T_collector, label="Temperatura do coletor")
plt.xlabel("Tempo (h)")
plt.ylabel("Temperatura (°C)")
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(np.arange(t0, tf + dt, dt) / 3600, Qa, label="Taxa de transferência de calor Qa")
plt.plot(np.arange(t0, tf + dt, dt) / 3600, Qcd_capa, label="Taxa de transferência de calor Qcd, capa")
plt.xlabel("Tempo (h)")
plt.ylabel("Taxa de calor (W)")
plt.legend()
plt.tight_layout()
plt.show()

print("d) Custo-benefício (COP) do sistema projetado:", COP)

https://docs.google.com/spreadsheets/d/1SZiYiDTN3kYUohvsHGzOkk8snQ14u83XhqQ4trQCZ-Y/edit?usp=sharing