program main_ET_PTC_HTF
implicit none
! ======= Constantes físicas y de geometría =======
Double Precision :: G_GHI, u_Air, rho_Water
Double Precision :: R_conv_PTC, R_conv_PV, R_conv_abs
Double Precision :: R_cond_PV, R_cond_sub, R_cond_abs
Double Precision :: C_p_HTF
Double Precision :: G_DNI
Double Precision :: T_a, T_sky, T_PVi, T_in_HTF
Double Precision :: hh_abs, hh_PTC, hh_PV, hh_HTF
Double Precision :: h_abs, h_PTC, h_PV, h_HTF
Double Precision :: D,W_PV,W_abs,W_PTC, L_PTC,L_SRC_PVT
Double Precision :: A_PV,A_abs,A_sub,A_ap,A_PTC, CR_PTC
Double Precision :: alpha_PV,alpha_abs,alpha_PTC
Double Precision :: epsilon_PV,epsilon_abs,epsilon_PTC
Double Precision :: IAM_elec,IAM_th,eta_opt,sigma
Double Precision :: P_Air,P_HTF,m_dot_HTF, th_PV,th_abs,th_sub,k_PV,k_abs,k_sub
Double Precision :: r_pipe,L_ver, L_edg, th_min,th_max
Double Precision :: A_cs,u_Water, A_hx,NTU,E
Double Precision :: eta_PV, P_PV, q_dot_HTF, eta_elec, eta_th
! ======= Variables de iteración y resultados =======
Double Precision :: T_PTC, T_abs, T_sub, T_out_HTF, T_PV_calc
integer :: iter1, iter2
! !!! NUEVAS VARIABLES PARA EL BUCLE GLOBAL !!!
integer :: iter_global
integer, parameter :: max_iter_global = 100
Double Precision :: tol_global, err_global, relax_factor
! ===== Inicialización de constantes =====
! ---- Boundary conditions
G_GHI = 1000.0d0
G_DNI = 910.0d0
u_Air = 5.0d0
T_a = 18.0d0 + 273.15d0
T_sky = T_a - 5.0d0
T_in_HTF = 20.0d0 + 273.15d0
! --- SEMILLA INICIAL ---
T_PVi = 355.45d0
! ---- Geometry
D = 0.04d0
W_PV = 0.06d0
W_abs = 0.06d0
W_PTC = 1.97d0
L_PTC = 5.08d0
L_SRC_PVT = L_PTC
A_PV = W_PV * L_SRC_PVT
A_abs = A_PV
A_sub = A_PV !+ A_abs
A_ap = W_PTC * L_PTC
A_PTC = 3.0d0 * L_PTC ! Verificar ese "3" para el área del PTC.
CR_PTC = A_ap / A_PV
! ---- Design parameters
alpha_PV = 0.97d0
alpha_abs = 0.90d0
alpha_PTC = 0.03d0
epsilon_PV = 0.2d0
epsilon_abs = 0.2d0
epsilon_PTC = 0.3d0
IAM_elec = 0.72d0 !0.28
IAM_th = 0.86d0 !0.14
eta_opt = 0.645d0
sigma = 5.67d-8
P_Air = 1.01325
P_HTF = 0.3119
m_dot_HTF = 0.0722d0
th_PV = 0.003d0
th_abs = 0.003d0
! th_sub = Calculo del triangulo inscrito
k_PV = 50.0d0
k_abs = 205.0d0
k_sub = 250.0d0
! ---- Coefficients (Constant assumption for now)
C_p_HTF = 4187.0d0
rho_Water = 977.7d0 ! Asumido standard si no estaba definido
hh_abs = 33.73d0 !3.438d0
hh_PTC = 7.54d0 !0.7687d0
hh_PV = 23.85d0 !2.431d0
hh_HTF = 1681.02d0 !835.3d0
h_abs = hh_abs
h_PTC = hh_PTC
h_PV = hh_PV
h_HTF = hh_HTF
! ======================= Thickness of the substrate
r_pipe = D/2.0d0
L_ver = (sqrt(3.0d0)/3.0d0) * W_abs
L_edg = (sqrt(3.0d0)/6.0d0) * W_abs
th_min = L_ver - r_pipe
th_max = L_edg - r_pipe
th_sub = (th_min + 2.0d0 * th_max)/3.0d0
! ======================= Thermal resistances
R_cond_PV = th_PV/(k_PV * A_PV)
R_cond_abs = th_abs/(k_abs * A_abs)
R_cond_sub = th_sub/(k_sub * A_sub)
R_conv_PV = 1.0d0*(h_PV * A_PV) ! 1.0d0/(h_PV * A_PV)
R_conv_abs = 1.0d0*(h_abs * A_abs) ! 1.0d0/(h_abs * A_abs)
R_conv_PTC = 1.0d0*(h_PTC * A_PTC) ! 1.0d0/(h_PTC * A_PTC)
! ======================= Water (HTF) velocity
A_cs = (3.1416d0 * (D/2.0d0)**(2))
u_Water = (m_dot_HTF)/(rho_Water * A_cs)
! ======================= Heat exchanger efficiency (NTU)
A_hx = 3.1416d0 * D * L_SRC_PVT
NTU = ((1.0d0/((1.0d0/h_HTF) + R_cond_sub)) * A_hx)/(m_dot_HTF * (C_p_HTF))
E = 1.0d0 - exp(- NTU)
! =========================================================================
! COMIENZO DEL BUCLE ITERATIVO GLOBAL
! =========================================================================
tol_global = 1.0d-4
relax_factor = 0.5d0 ! 50% valor nuevo, 50% valor viejo para estabilidad
print *, ">>> INICIO BUCLE CONVERGENCIA T_PV <<<"
print *, "Iter | T_semilla (K) | T_calc (K) | Error"
do iter_global = 1, max_iter_global
! 1. Calcular Potencia Electrica usando la SEMILLA actual (T_PVi)
! IMPORTANTE: Usamos T_PVi porque T_PV_calc aún no existe en esta vuelta
eta_PV = 0.298d0 + 0.0142d0 * (log(CR_PTC)) + &
(- 0.000715d0 + 0.0000697d0 * (log(CR_PTC))) * (T_PVi - 298.0d0)
P_PV = G_DNI * A_PV * CR_PTC * eta_opt * IAM_elec * eta_PV
! 2. Paso 1: Calcular T_PTC
call calc_T_PTC(T_a, T_sky, T_PVi, T_PTC, sigma, G_GHI, A_PTC, alpha_PTC, &
epsilon_PV, A_PV, epsilon_PTC, R_conv_PTC, iter1)
! 3. Paso 2: Calcular Sistema HTF (T_abs, T_sub, T_out_HTF)
call solve_HTF_system(T_a, T_sky, T_PVi, T_in_HTF, T_PTC, T_abs, T_sub, T_out_HTF, &
E, sigma, G_GHI, &
A_abs, alpha_abs, epsilon_abs, A_PV, epsilon_PV, &
R_conv_abs, R_conv_PV, R_cond_abs, R_cond_sub, &
G_DNI, CR_PTC, eta_opt, IAM_th, alpha_PV, P_PV, &
m_dot_HTF, C_p_HTF, h_HTF, iter2)
! 4. Paso 3: Calcular nueva T_PV (Resultado del balance)
call calc_T_PV(T_abs, T_sub, T_out_HTF, T_in_HTF, R_cond_PV, R_cond_sub, R_cond_abs, &
m_dot_HTF, C_p_HTF, T_PV_calc)
! 5. Verificar Convergencia
err_global = abs(T_PV_calc - T_PVi)
write(*, '(I4, 2X, F10.4, 2X, F10.4, 2X, E10.3)') iter_global, T_PVi, T_PV_calc, err_global
if (err_global < tol_global) then
print *, ">>> CONVERGENCIA ALCANZADA <<<"
exit
end if
! 6. Actualizar semilla para la siguiente vuelta (con relajación)
T_PVi = (1.0d0 - relax_factor)*T_PVi + (relax_factor*T_PV_calc)
end do
! =========================================================================
! Resultados Finales
q_dot_HTF = (T_out_HTF - T_in_HTF) * m_dot_HTF * C_p_HTF
eta_elec = (P_PV)/(G_DNI * A_ap)
eta_th = (q_dot_HTF)/(G_DNI * A_ap)
print *, " "
print *, "====== RESULTADOS FINALES ======"
print *, "T_PV Inicial (C):", 355.45 - 273.15d0 !T_PVi
print *, "T_PV Final (C) :", T_PV_calc - 273.15d0
print *, "T_abs (C) :", T_abs - 273.15d0
print *, "T_sub (C) :", T_sub - 273.15d0
print *, "T_out_HTF (C) :", T_out_HTF - 273.15d0
print *, "T_PTC (C) :", T_PTC - 273.15d0
print *, "Potencia PV (W) :", P_PV
print *, "Eta Elec (-) :", eta_elec
print *, "Eta Termica (-) :", eta_th
print *, "================================"
contains
!-------------------------------------------
subroutine calc_T_PTC(T_a, T_sky, T_PVi, T_PTC, sigma, G_GHI, A_PTC, alpha_PTC, epsilon_PV, A_PV, &
epsilon_PTC, R_conv_PTC, iter)
implicit none
Double Precision, intent(out) :: T_PTC
Double Precision, intent(in) :: T_a, T_sky, T_PVi
Double Precision, intent(in) :: sigma, G_GHI, A_PTC, alpha_PTC
Double Precision, intent(in) :: epsilon_PV, A_PV, epsilon_PTC, R_conv_PTC
integer, intent(out) :: iter
Double Precision :: f_T_PTC, df_T_PTC, T_PTCe, tol
integer :: max_iter
! Reiniciar si viene con basura, o mantener valor anterior para velocidad
if (T_PTC < 100.0d0) T_PTC = 320.0d0
tol = 1.0d-6
max_iter = 100
do iter = 1, max_iter
f_T_PTC = G_GHI * A_PTC * alpha_PTC &
- (T_PTC**4 - T_sky**4) * A_PTC * epsilon_PTC * sigma &
+ (T_PVi**4 - T_PTC**4) * A_PV * epsilon_PV * sigma &
- ((T_PTC - T_a)/(R_conv_PTC))
df_T_PTC = - 4.0d0 * T_PTC**3 * A_PTC * epsilon_PTC * sigma &
- 4.0d0 * T_PTC**3 * A_PV * epsilon_PV * sigma &
- 1.0d0/R_conv_PTC
T_PTCe = T_PTC - f_T_PTC/df_T_PTC
if (abs(T_PTCe - T_PTC) < tol) then
T_PTC = T_PTCe
exit
endif
T_PTC = T_PTCe
end do
end subroutine calc_T_PTC
!-------------------------------------------
subroutine solve_HTF_system(T_a, T_sky, T_PVi, T_in_HTF, T_PTC, T_abs, T_sub, T_out_HTF, E, sigma, G_GHI, &
A_abs, alpha_abs, epsilon_abs, A_PV, epsilon_PV, &
R_conv_abs, R_conv_PV, R_cond_abs, R_cond_sub, &
G_DNI, CR_PTC, eta_opt, IAM_th, alpha_PV, P_PV, &
m_dot_HTF, C_p_HTF, h_HTF, iter)
implicit none
Double Precision, intent(in) :: T_a, T_sky, T_PVi, T_in_HTF, T_PTC, sigma, G_GHI
Double Precision, intent(in) :: A_abs, alpha_abs, epsilon_abs, A_PV, epsilon_PV
Double Precision, intent(in) :: R_conv_abs, R_conv_PV, R_cond_abs, R_cond_sub
Double Precision, intent(in) :: E, m_dot_HTF, C_p_HTF, h_HTF
Double Precision, intent(in) :: G_DNI, CR_PTC, eta_opt, IAM_th, alpha_PV, P_PV
integer, intent(out) :: iter
Double Precision, intent(out) :: T_abs, T_sub, T_out_HTF
Double Precision :: f1, f2, f3, J(3,3), rhs(3), dx(3)
Double Precision :: tol
integer :: max_iter
! Inicializar guesses razonables si es la primera vez
if (T_abs < 100.d0) T_abs = 350.0d0
if (T_sub < 100.d0) T_sub = 350.0d0
if (T_out_HTF < 100.d0) T_out_HTF = 350.0d0
tol = 1.0d-8
max_iter = 100
do iter = 1, max_iter
f1 = G_GHI * A_abs * alpha_abs &
+ G_DNI * A_PV * alpha_PV * CR_PTC * eta_opt * IAM_th &
- P_PV &
- (T_abs**4 - T_sky**4) * A_abs * epsilon_abs * sigma &
- (T_PVi**4 - T_PTC**4) * A_PV * epsilon_PV * sigma &
- ((T_abs - T_a)/(R_conv_abs)) &
- ((T_PVi - T_a)/(R_conv_PV)) &
- (T_out_HTF - T_in_HTF) * m_dot_HTF * C_p_HTF
f2 = G_GHI * A_abs * alpha_abs &
- (T_abs**4 - T_sky**4) * A_abs * epsilon_abs * sigma &
- ((T_abs - T_a)/(R_conv_abs)) &
- ((T_abs - T_sub)/(R_cond_abs + R_cond_sub))
f3 = (T_out_HTF - T_in_HTF) * m_dot_HTF * C_p_HTF &
- (T_sub - T_in_HTF) * E * m_dot_HTF * C_p_HTF
! Jacobiano
J(1,1) = - 4.0d0 * T_abs**3 * A_abs * epsilon_abs * sigma - 1.0d0/R_conv_abs
J(1,2) = 0.0d0
J(1,3) = - m_dot_HTF * C_p_HTF
J(2,1) = - 4.0d0 * T_abs**3 * A_abs * epsilon_abs * sigma - 1.0d0/R_conv_abs - 1.0d0/(R_cond_abs + R_cond_sub)
J(2,2) = - 1.0d0/(R_cond_abs + R_cond_sub)
J(2,3) = 0.0d0
J(3,1) = 0.0d0
J(3,2) = - E * m_dot_HTF * C_p_HTF
J(3,3) = m_dot_HTF * C_p_HTF
rhs(1) = -f1
rhs(2) = -f2
rhs(3) = -f3
call solve3x3(J, rhs, dx)
T_abs = T_abs + dx(1)
T_sub = T_sub + dx(2)
T_out_HTF = T_out_HTF + dx(3)
if (max(abs(dx(1)), max(abs(dx(2)), abs(dx(3)))) < tol) exit
end do
end subroutine solve_HTF_system
!-------------------------------------------
subroutine solve3x3(A, b, x)
implicit none
Double Precision, intent(inout) :: A(3,3)
Double Precision, intent(in) :: b(3)
Double Precision, intent(out) :: x(3)
Double Precision :: M(3,3), rhs(3)
integer :: i,j,k
Double Precision :: factor
M = A
rhs = b
do k = 1,2
if (abs(M(k,k)) < 1.0d-14) stop "Pivote casi cero"
do i = k+1, 3
factor = M(i,k)/M(k,k)
do j = k,3
M(i,j) = M(i,j) - factor*M(k,j)
end do
rhs(i) = rhs(i) - factor*rhs(k)
end do
end do
x(3) = rhs(3)/M(3,3)
x(2) = (rhs(2) - M(2,3)*x(3))/M(2,2)
x(1) = (rhs(1) - M(1,2)*x(2) - M(1,3)*x(3))/M(1,1)
end subroutine solve3x3
!-------------------------------------------
subroutine calc_T_PV(T_abs, T_sub, T_out_HTF, T_in_HTF, R_cond_PV, R_cond_sub, R_cond_abs, &
m_dot_HTF, C_p_HTF, T_PV)
implicit none
Double Precision, intent(in) :: T_abs, T_sub,T_out_HTF, T_in_HTF
Double Precision, intent(in) :: R_cond_PV, R_cond_sub, R_cond_abs
Double Precision, intent(in) :: m_dot_HTF, C_p_HTF
Double Precision, intent(out) :: T_PV
Double Precision :: term1, term2
term1 = (T_out_HTF - T_in_HTF) * m_dot_HTF * C_p_HTF
term2 = (T_abs - T_sub)/(R_cond_abs + R_cond_sub)
T_PV = T_sub + (term1 - term2)*(R_cond_PV + R_cond_sub)
end subroutine calc_T_PV
end program main_ET_PTC_HTF
To embed this program on your website, copy the following code and paste it into your website's HTML: