program solve_temperatures
implicit none
! Declaración de variables
real :: Tem_a, Tem_sky, Tem_f_in, Tem_pv, Tem_conc
real :: Tem_abs, Tem_sub, Tem_f_out, P_ele_pv, q_f
real :: tol
integer :: max_iter
! Variables de las ecuaciones internas
real :: q_solar_pv, q_solar_abs, q_rad_pv, q_rad_abs
real :: q_conv_pv, q_conv_abs, q_cond_abs_x_sub
! Valores de las constantes del problema
real :: G_tot, G_b, A_pv, A_abs, A_conc, alpha_pv, alpha_abs, alpha_conc, eta_opt, eta_pv, IAM_th, IAM_elec
real :: CR_pv, sigma, epsilon_pv, epsilon_abs, epsilon_conc, m_dot_f, C_p_f, E
real :: R_conv_pv, R_conv_abs, R_cond_abs, R_cond_sub
! Inicialización de las constantes
G_tot = 1000
G_b = 800
A_pv = 1.2
A_abs = 0.6
A_conc = 30
alpha_pv = 0.97
alpha_abs = 0.9
alpha_conc = 0.03
eta_opt = 0.83
eta_pv = 0.2998
IAM_th = 0.9
IAM_elec = 0.95
CR_pv = 10
sigma = 5.67e-8
epsilon_pv = 0.2
epsilon_abs = 0.2
epsilon_conc = 0.3
m_dot_f = 0.15
C_p_f = 4187
E = 0.9646
R_conv_pv = 0.02497
R_conv_abs = 0.03532
R_cond_abs = 0.00002439
R_cond_sub = 0.0001667
! Inicialización de las temperaturas conocidas
Tem_a = 298.2
Tem_sky = 298.2
Tem_f_in = 343.2
Tem_pv = 355.5
Tem_conc = 307
! Parámetros para el método de Newton-Raphson
tol = 1.0e-6
max_iter = 100
! Estimación inicial de las temperaturas a resolver
Tem_abs = 350.0
Tem_sub = 340.0
Tem_f_out = 320.0
! Llamar a la subrutina para resolver el sistema de ecuaciones
call newton_raphson(Tem_abs, Tem_sub, Tem_f_out, tol, max_iter)
! Imprimir los resultados
print *, 'Tem_abs: ', Tem_abs
print *, 'Tem_sub: ', Tem_sub
print *, 'Tem_f_out: ', Tem_f_out
contains
! Función para calcular el sistema de ecuaciones
subroutine equations(Tem_abs, Tem_sub, Tem_f_out, f)
real, intent(in) :: Tem_abs, Tem_sub, Tem_f_out
real, intent(out) :: f(3)
real :: q_solar_pv, q_solar_abs, q_rad_pv, q_rad_abs
real :: q_conv_pv, q_conv_abs, q_cond_abs_x_sub, q_f
q_solar_pv = G_b * A_pv * CR_pv * eta_opt * IAM_th * alpha_pv
q_solar_abs = G_tot * A_abs * alpha_abs
P_ele_pv = G_b * A_pv * CR_pv * eta_opt * IAM_elec * eta_pv
q_rad_pv = sigma * epsilon_pv * A_pv * (Tem_pv**4 - Tem_conc**4)
q_rad_abs = sigma * epsilon_abs * A_abs * (Tem_abs**4 - Tem_sky**4)
q_conv_pv = (Tem_pv - Tem_a) * R_conv_pv
q_conv_abs = (Tem_abs - Tem_a) * R_conv_abs
q_cond_abs_x_sub = (Tem_abs - Tem_sub) / (R_cond_abs + R_cond_sub)
q_f = m_dot_f * C_p_f * (Tem_f_out - Tem_f_in)
f(1) = q_solar_pv - P_ele_pv - q_conv_pv - q_rad_pv + q_solar_abs - q_conv_abs - q_rad_abs - q_f
f(2) = q_solar_abs - q_conv_abs - q_rad_abs - q_cond_abs_x_sub
f(3) = m_dot_f * C_p_f * (Tem_f_out - Tem_f_in) - E * m_dot_f * C_p_f * (Tem_sub - Tem_f_in)
end subroutine equations
! Función para calcular la jacobiana
subroutine jacobian(Tem_abs, Tem_sub, Tem_f_out, J)
real, intent(in) :: Tem_abs, Tem_sub, Tem_f_out
real, intent(out) :: J(3, 3)
real :: df_dTem_abs, df_dTem_sub, df_dTem_f_out
real :: h = 1.0e-6
real :: f1(3), f2(3), f3(3)
call equations(Tem_abs, Tem_sub, Tem_f_out, f1)
call equations(Tem_abs + h, Tem_sub, Tem_f_out, f2)
J(1, 1) = (f2(1) - f1(1)) / h
J(2, 1) = (f2(2) - f1(2)) / h
J(3, 1) = (f2(3) - f1(3)) / h
call equations(Tem_abs, Tem_sub + h, Tem_f_out, f2)
J(1, 2) = (f2(1) - f1(1)) / h
J(2, 2) = (f2(2) - f1(2)) / h
J(3, 2) = (f2(3) - f1(3)) / h
call equations(Tem_abs, Tem_sub, Tem_f_out + h, f2)
J(1, 3) = (f2(1) - f1(1)) / h
J(2, 3) = (f2(2) - f1(2)) / h
J(3, 3) = (f2(3) - f1(3)) / h
end subroutine jacobian
! Subrutina para resolver el sistema de ecuaciones usando el método de Newton-Raphson
subroutine newton_raphson(Tem_abs, Tem_sub, Tem_f_out, tol, max_iter)
real, intent(inout) :: Tem_abs, Tem_sub, Tem_f_out
real, intent(in) :: tol
integer, intent(in) :: max_iter
real :: f(3), J(3, 3), dx(3)
real :: norm_f
integer :: iter
do iter = 1, max_iter
call equations(Tem_abs, Tem_sub, Tem_f_out, f)
call jacobian(Tem_abs, Tem_sub, Tem_f_out, J)
! Resolver el sistema lineal J * dx = -f
call solve_linear_system(J, f, dx)
! Actualizar las variables
Tem_abs = Tem_abs + dx(1)
Tem_sub = Tem_sub + dx(2)
Tem_f_out = Tem_f_out + dx(3)
! Calcular la norma de f para el criterio de convergencia
norm_f = sqrt(sum(f**2))
if (norm_f < tol) then
return
end if
end do
print *, 'Advertencia: Número máximo de iteraciones alcanzado'
end subroutine newton_raphson
! Subrutina para resolver un sistema lineal A * x = b
subroutine solve_linear_system(A, b, x)
real, intent(in) :: A(3, 3), b(3)
real, intent(out) :: x(3)
real :: A_inv(3, 3)
integer :: i, j, k
! Calcular la inversa de A
call invert_matrix(A, A_inv)
! Calcular x = A_inv * b
do i = 1, 3
x(i) = 0.0
do j = 1, 3
x(i) = x(i) + A_inv(i, j) * b(j)
end do
end do
end subroutine solve_linear_system
! Subrutina para invertir una matriz 3x3
subroutine invert_matrix(A, A_inv)
real, intent(in) :: A(3, 3)
real, intent(out) :: A_inv(3, 3)
real :: det
integer :: i, j
det = A(1, 1) * (A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2)) - &
A(1, 2) * (A(2, 1) * A(3, 3) - A(2, 3) * A(3, 1)) + &
A(1, 3) * (A(2, 1) * A(3, 2) - A(2, 2) * A(3, 1))
if (abs(det) < 1.0e-10) then
print *, 'Error: Matriz singular'
stop
end if
A_inv(1, 1) = (A(2, 2) * A(3, 3) - A(2, 3) * A(3, 2)) / det
A_inv(1, 2) = -(A(1, 2) * A(3, 3) - A(1, 3) * A(3, 2)) / det
A_inv(1, 3) = (A(1, 2) * A(2, 3) - A(1, 3) * A(2, 2)) / det
A_inv(2, 1) = -(A(2, 1) * A(3, 3) - A(2, 3) * A(3, 1)) / det
A_inv(2, 2) = (A(1, 1) * A(3, 3) - A(1, 3) * A(3, 1)) / det
A_inv(2, 3) = -(A(1, 1) * A(2, 3) - A(1, 3) * A(2, 1)) / det
A_inv(3, 1) = (A(2, 1) * A(3, 2) - A(2, 2) * A(3, 1)) / det
A_inv(3, 2) = -(A(1, 1) * A(3, 2) - A(1, 2) * A(3, 1)) / det
A_inv(3, 3) = (A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1)) / det
end subroutine invert_matrix
end program solve_temperatures
To embed this program on your website, copy the following code and paste it into your website's HTML: