program main_ET_PTC_HTF
implicit none
! ======= Constantes físicas y de geometría =======
Double Precision :: sigma, G_GHI
Double Precision :: A_PTC, alpha_PTC, epsilon_PTC
Double Precision :: A_PV, epsilon_PV, R_conv_PTC, R_conv_PV
Double Precision :: A_abs, A_sub, alpha_abs, epsilon_abs, R_conv_abs
Double Precision :: R_cond_PV, R_cond_sub, R_cond_abs
Double Precision :: m_dot_HTF, C_p_HTF, h_HTF
Double Precision :: A_hx, UA, NTU, E
Double Precision :: G_DNI, CR_PTC, eta_opt, IAM_th, alpha_PV, P_PV
Double Precision :: T_a, T_sky, T_PVi, T_in_HTF
Double Precision :: hh_abs, hh_PTC, hh_PV, hh_HTF
Double Precision :: th_PV,th_abs,th_sub,k_PV,k_abs,k_sub
! ======= Variables de iteración y resultados =======
Double Precision :: T_PTC, T_abs, T_sub, T_out_HTF, T_PV
integer :: iter1, iter2
! ===== Inicialización de constantes =====
sigma = 5.670374d-8
G_GHI = 1000.0d0
T_a = 298.2d0
T_sky = 298.2d0
T_in_HTF = 343.2d0
T_PVi = 355.5d0
A_PV = 1.2d0
A_abs = 0.6d0
A_sub = A_PV + A_abs !OJO aquí
A_PTC = 30.0d0
alpha_abs = 0.9d0
alpha_PTC = 0.03d0
epsilon_PV = 0.2d0
epsilon_abs = 0.2d0
epsilon_PTC = 0.3d0
m_dot_HTF = 0.15d0
C_p_HTF = 4187.0d0
h_HTF = 1665.0d0
G_DNI = 800.d0
CR_PTC = 10.d0
eta_opt = 0.83d0
IAM_th = 0.86d0
alpha_PV = 0.97d0
!R_conv_PTC = 0.04336d0
!R_conv_PV = 0.3428d0
!R_conv_abs = 0.0464d0
!R_cond_PV = 0.00005d0
!R_cond_sub = 0.00003082d0
!R_cond_abs = 0.00002439d0
hh_abs = 35.22d0
hh_PTC = 7.87d0
hh_PV = 24.90d0
hh_HTF = 1684.90d0
th_PV = 0.003d0 ! PV thickness
th_abs = 0.003d0 ! Absorber thickness
th_sub = 0.02541d0 ! Substrate thickness
k_PV = 50.0d0 ! Coefficient of thermal conductivity of PV
k_abs = 205.0d0 ! Coefficient of thermal conductivity of absorber
k_sub = 250.0d0 ! Coefficient of thermal conductivity of substrate
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 / ( hh_PV * A_PV )
R_conv_abs = 1 / ( hh_abs * A_abs )
R_conv_PTC = 1 / ( hh_PTC * A_PTC )
! ========= NTU
A_hx = 3.1416d0 * 0.03d0 * 10.0d0
UA = 1.0d0 / (1.0d0/h_HTF + R_cond_sub)
NTU = (UA * A_hx) / (m_dot_HTF * C_p_HTF)
E = 1.0d0 - exp(-NTU)
P_PV = 1723.0d0
! ===== 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)
print *, "-------------------------------------------"
print *, "Temperatura PTC calculada:", T_PTC, "K", T_PTC-273.15d0, "°C"
print *, "Iteraciones método 1:", iter1
print *, "-------------------------------------------"
! ===== Paso 2: Calcular T_abs, T_sub, T_out_HTF y E =====
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)
print *, "-------------------------------------------"
print *, "T_abs (K) =", T_abs, " (°C) =", T_abs-273.15d0
print *, "T_sub (K) =", T_sub, " (°C) =", T_sub-273.15d0
print *, "T_out_HTF (K) =", T_out_HTF, " (°C) =", T_out_HTF-273.15d0
print *, "E =", E
print *, "Iteraciones método 2:", iter2
print *, "-------------------------------------------"
! ===== Paso 3: Calcular nueva T_PV =====
call calc_T_PV(T_abs, T_sub, T_in_HTF, R_cond_PV, R_cond_sub, R_cond_abs, &
E, m_dot_HTF, C_p_HTF, T_PV)
print *, "-------------------------------------------"
print *, "Temperatura PV calculada:", T_PV, "K", T_PV-273.15d0, "°C"
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
!T_PVi = 355.5d0
!T_a = 298.2d0
!T_sky = 298.2d0
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 & ! q_dot_sol_PTC
- (T_PTC**4 - T_sky**4) * A_PTC * epsilon_PTC * sigma & ! q_dot_rad_PTC
+ (T_PVi**4 - T_PTC**4) * A_PV * epsilon_PV * sigma & ! q_dot_rad_PV
- ((T_PTC - T_a)/(R_conv_PTC)) ! q_dot_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 :: P_PV
Double Precision :: f1, f2, f3, J(3,3), rhs(3), dx(3)
Double Precision :: UA, A_hx, NTU, tol
integer :: max_iter
!T_a = 298.2d0
!T_sky = 298.2d0
!T_PVi = 355.5d0
!q_dot_sol_PV = 6647.0d0
!P_PV = 1723.0d0
T_abs = 350.0d0
T_sub = 350.0d0
T_out_HTF = 350.0d0
!A_hx = 3.1416d0 * 0.03d0 * 10.0d0
!UA = 1.0d0 / (1.0d0/h_HTF + R_cond_sub)
!NTU = (UA * A_hx) / (m_dot_HTF * C_p_HTF)
!E = 1.0d0 - exp(-NTU)
tol = 1.0d-8
max_iter = 100
do iter = 1, max_iter
f1 = G_GHI * A_abs * alpha_abs & ! q_dot_sol_abs
+ G_DNI * A_PV * alpha_PV * CR_PTC * eta_opt * IAM_th & ! q_dot_sol_PV
- P_PV & ! P_PV
- (T_abs**4 - T_sky**4) * A_abs * epsilon_abs * sigma & ! q_dot_rad_abs
- (T_PVi**4 - T_PTC**4) * A_PV * epsilon_PV * sigma & ! q_dot_rad_PV
- ((T_abs - T_a)/(R_conv_abs)) & ! q_dot_conv_abs
- ((T_PVi - T_a)/(R_conv_PV)) & ! q_dot_conv_PV
!- ((T_abs - T_sub)/(R_cond_abs + R_cond_sub)) & ! q_dot_cond_abs_x_sub
- (T_sub - T_in_HTF) * E * m_dot_HTF * C_p_HTF ! q_dot_HTF, ojo
f2 = G_GHI * A_abs * alpha_abs & ! q_dot_sol_abs
- (T_abs**4 - T_sky**4) * A_abs * epsilon_abs * sigma & ! q_dot_rad_abs
- ((T_abs - T_a)/(R_conv_abs)) & ! q_dot_conv_abs
- ((T_abs - T_sub)/(R_cond_abs + R_cond_sub)) ! q_dot_cond_abs_x_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) = -1.0d0/R_conv_abs - 4.0d0*epsilon_abs*sigma*A_abs*T_abs**3 - 1.0d0/(R_cond_abs + R_cond_sub)
J(1,2) = 1.0d0/(R_cond_abs + R_cond_sub) - E*m_dot_HTF*C_p_HTF
J(1,3) = 0.0d0
J(2,1) = -1.0d0/R_conv_abs - 4.0d0*epsilon_abs*sigma*A_abs*T_abs**3 - 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_in_HTF, R_cond_PV, R_cond_sub, R_cond_abs, &
E, m_dot_HTF, C_p_HTF, T_PV)
implicit none
Double Precision, intent(in) :: T_abs, T_sub, T_in_HTF
Double Precision, intent(in) :: R_cond_PV, R_cond_sub, R_cond_abs
Double Precision, intent(in) :: E, m_dot_HTF, C_p_HTF
Double Precision, intent(out) :: T_PV
Double Precision :: term1, term2
term1 = E*m_dot_HTF*C_p_HTF*(T_sub - T_in_HTF) ! q_dot_HTF, ojo
term2 = (T_abs - T_sub)/(R_cond_abs + R_cond_sub) ! q_dot_cond_abs_x_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: