program main_NR_total
implicit none
Double Precision :: Tabs, Tsub, Tout, TPTC
Double Precision :: G_GHI, G_DNI, Ta, Tsky, Tin, T_PVi
Double Precision :: W_PV, W_abs, L, A_PV, A_abs, A_ap, A_PTC
Double Precision :: CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs
Double Precision :: IAM_elec, IAM_th, eta_opt, sigma
Double Precision :: mdot, Cp, U, A_hx, P_PV, T_PV
Double Precision :: Rcond_abs, Rcond_PV, Rcond_sub
Double Precision :: Rconv_abs, Rconv_PV, Rconv_PTC
Double Precision :: alpha_PTC, epsilon_PTC
integer :: iter
! ==== Valores fijos ====
! ---- ---- ---- ---- Boundary conditions
G_GHI = 1000.0d0
G_DNI = 800.0d0
Ta = 25.0d0 + 273.15d0
Tsky = 25.0d0 + 273.15d0
Tin = 70.0d0 + 273.15d0
T_PVi = 355.45d0
! ---- ---- ---- ---- Geometry of the PTC and SRC-PVT
W_PV = 0.12d0
W_abs = 0.06d0
L = 10.0d0
A_PV = W_PV*L
A_abs = W_abs*L
A_ap = 1.2d0*L
A_PTC = 3.0d0*L
CR_PTC = A_ap/A_PV
! ---- ---- ---- ---- Design parameters
alpha_PV = 0.97d0
alpha_abs = 0.90d0
alpha_PTC = 0.03d0
eps_PV = 0.2d0
eps_abs = 0.2d0
epsilon_PTC = 0.3d0
IAM_elec = 0.72d0
IAM_th = 0.86d0
eta_opt = 0.83d0
sigma = 5.67d-8 !5.670374d-8
mdot = 0.15d0
Cp = 4187.0d0
U=780.1d0
A_hx=0.9425d0
P_PV=1723.0d0
Rcond_abs=0.00002439d0
Rcond_PV=0.00005d0
Rcond_sub=0.00008472d0
Rconv_abs=0.4848d0
Rconv_PV=0.3428d0
Rconv_PTC=0.04336d0
! ==== Paso 1: Calcular T_PTC (Newton-Raphson 1 variable) ====
call calc_T_PTC(Ta, Tsky, T_PVi, TPTC, sigma, G_GHI, A_PTC, alpha_PTC, &
eps_PV, A_PV, epsilon_PTC, Rconv_PTC, iter)
print *, "--------------------------------------"
print *, "Cálculo de T_PTC"
print *, "Iteraciones: ", iter
print *, "T_PTC [K]: ", TPTC
! ==== Paso 2: Calcular Tabs, Tsub, Tout (Newton-Raphson 3 variables) ====
call solve_NR_3x3(Tabs, Tsub, Tout, &
G_GHI, G_DNI, Ta, Tsky, Tin, T_PVi, TPTC, &
A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
sigma, mdot, Cp, U, A_hx, P_PV, &
Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)
print *, "--------------------------------------"
print *, "Resultados Newton–Raphson 3x3:"
print *, "T_abs [K]: ", Tabs
print *, "T_sub [K]: ", Tsub
print *, "T_outHTF[K]: ", Tout
print *, "--------------------------------------"
! ==== Paso 3: Calcular T_PV con balance adicional ====
call calc_T_PV(Tabs, Tsub, Tout, Tin, Rcond_PV, Rcond_sub, Rcond_abs, &
mdot, Cp, T_PV)
print *, "--------------------------------------"
print *, "Cálculo de T_PV:"
print *, "T_PV [K]: ", T_PV
print *, "--------------------------------------"
end program main_NR_total
!===============================================================
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_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_NR_3x3(Tabs, Tsub, Tout, &
G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
sigma, mdot, Cp, U, A_hx, P_PV, &
Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)
implicit none
Double Precision, intent(out) :: Tabs, Tsub, Tout
Double Precision, intent(in) :: G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC
Double Precision, intent(in) :: A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs
Double Precision, intent(in) :: eps_PV, eps_abs, sigma
Double Precision, intent(in) :: mdot, Cp, U, A_hx, P_PV
Double Precision, intent(in) :: Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV
Integer :: it, maxit
Double Precision :: tol, normF
Double Precision :: F(3), dx(3), J(3,3)
tol = 1.0d-8
maxit = 100
Tabs = 355.0d0
Tsub = 355.0d0
Tout = 355.0d0
do it = 1, maxit
call residual_and_jacobian(Tabs, Tsub, Tout, F, J, &
G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
sigma, mdot, Cp, U, A_hx, P_PV, &
Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)
normF = max( max(abs(F(1)),abs(F(2))), abs(F(3)) )
if (normF < tol) exit
call solve_3x3(J, -F, dx)
Tabs = Tabs + dx(1)
Tsub = Tsub + dx(2)
Tout = Tout + dx(3)
end do
print *, "Iteraciones: ", it-1
print *, "||F||_inf : ", normF
end subroutine solve_NR_3x3
!===============================================================
subroutine residual_and_jacobian(Tabs, Tsub, Tout, F, J, &
G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC, &
A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs, eps_PV, eps_abs, &
sigma, mdot, Cp, U, A_hx, P_PV, &
Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV)
implicit none
Double Precision, intent(in) :: Tabs, Tsub, Tout
Double Precision, intent(out):: F(3), J(3,3)
Double Precision, intent(in) :: G_GHI, G_DNI, Ta, Tsky, Tin, TPV, TPTC
Double Precision, intent(in) :: A_PV, A_abs, CR_PTC, alpha_PV, alpha_abs
Double Precision, intent(in) :: eps_PV, eps_abs, sigma
Double Precision, intent(in) :: mdot, Cp, U, A_hx, P_PV
Double Precision, intent(in) :: Rcond_abs, Rcond_sub, Rconv_abs, Rconv_PV
Double Precision :: q_solar_abs, q_solar_ptc
Double Precision :: rad_abs, rad_pv, conv_pv, conv_abs
Double Precision :: DT, eps_log, denom, Lmtd, Rcond_sum
Rcond_sum = Rcond_abs + Rcond_sub
q_solar_abs = G_GHI*A_abs*alpha_abs
q_solar_ptc = G_DNI*A_PV*alpha_PV*CR_PTC*0.83d0*0.86d0
rad_abs = eps_abs*sigma*A_abs*(Tabs**4 - Tsky**4)
rad_pv = eps_PV *sigma*A_PV*(TPV**4 - TPTC**4)
conv_pv = (TPV - Ta)/Rconv_PV
conv_abs= (Tabs - Ta)/Rconv_abs
DT = Tout - Tin
eps_log = 1.0d-9
denom = max( (Tsub - Tout), eps_log )
Lmtd = log( max( (Tsub - Tin), eps_log ) / denom )
F(1) = q_solar_abs + q_solar_ptc - P_PV - rad_abs - rad_pv - conv_pv - conv_abs - mdot*Cp*DT
F(2) = q_solar_abs - rad_abs - conv_abs - (Tabs - Tsub)/Rcond_sum
F(3) = mdot*Cp*DT - U*A_hx * ( DT / Lmtd )
J(1,1) = -4.0d0*eps_abs*sigma*A_abs*(Tabs**3) - 1.0d0/Rconv_abs
J(1,2) = 0.0d0
J(1,3) = - mdot*Cp
J(2,1) = -4.0d0*eps_abs*sigma*A_abs*(Tabs**3) - 1.0d0/Rconv_abs - 1.0d0/Rcond_sum
J(2,2) = 1.0d0/Rcond_sum
J(2,3) = 0.0d0
J(3,1) = 0.0d0
J(3,2) = U*A_hx * DT / (Lmtd*Lmtd) * ( 1.0d0/(max(Tsub - Tin,eps_log)) - 1.0d0/(denom) )
J(3,3) = mdot*Cp - U*A_hx * ( ( Lmtd - DT/(denom) ) / (Lmtd*Lmtd) )
end subroutine residual_and_jacobian
!===============================================================
subroutine solve_3x3(A, b, x)
implicit none
Double Precision, intent(in) :: A(3,3), b(3)
Double Precision, intent(out) :: x(3)
Double Precision :: M(3,3), det
M = A
det = M(1,1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2)) &
- M(1,2)*(M(2,1)*M(3,3)-M(2,3)*M(3,1)) &
+ M(1,3)*(M(2,1)*M(3,2)-M(2,2)*M(3,1))
if (abs(det) < 1.0d-16) then
print *, "Jacobiano casi singular."
end if
x(1) = ( b(1)*(M(2,2)*M(3,3)-M(2,3)*M(3,2)) &
- b(2)*(M(1,2)*M(3,3)-M(1,3)*M(3,2)) &
+ b(3)*(M(1,2)*M(2,3)-M(1,3)*M(2,2)) ) / det
x(2) = ( - b(1)*(M(2,1)*M(3,3)-M(2,3)*M(3,1)) &
+ b(2)*(M(1,1)*M(3,3)-M(1,3)*M(3,1)) &
- b(3)*(M(1,1)*M(2,3)-M(1,3)*M(2,1)) ) / det
x(3) = ( b(1)*(M(2,1)*M(3,2)-M(2,2)*M(3,1)) &
- b(2)*(M(1,1)*M(3,2)-M(1,2)*M(3,1)) &
+ b(3)*(M(1,1)*M(2,2)-M(1,2)*M(2,1)) ) / det
end subroutine solve_3x3
!===============================================================
subroutine calc_T_PV(Tabs, Tsub, Tout, Tin, Rcond_PV, Rcond_sub, Rcond_abs, &
mdot, Cp, T_PV)
implicit none
Double Precision, intent(in) :: Tabs, Tsub, Tout, Tin
Double Precision, intent(in) :: Rcond_PV, Rcond_sub, Rcond_abs
Double Precision, intent(in) :: mdot, Cp
Double Precision, intent(out) :: T_PV
Double Precision :: term1, term2
term1 = mdot*Cp*(Tout - Tin) ! q_dot_HTF
term2 = (Tabs - Tsub)/(Rcond_abs + Rcond_sub) ! q_dot_cond_abs_x_sub
T_PV = Tsub + (term1 - term2)*(Rcond_PV + Rcond_sub)
end subroutine calc_T_PV
To embed this program on your website, copy the following code and paste it into your website's HTML: