diff --git a/app/main.f90 b/app/main.f90 index d10d718..4192e02 100644 --- a/app/main.f90 +++ b/app/main.f90 @@ -5,16 +5,13 @@ program main implicit none type(FluidDataOut) :: fluid_to_characterize - - fluid_to_characterize=characterize(file='YPF2.nml',mw_source="experimental",& - method = "plus_mw", pho_method = 1 , fix_C=.true., eos='PR') - + + fluid_to_characterize=characterize(file='YPF2.nml', mw_source="experimental",& + method = "plus_mw", rho_method = 1 , fix_C=.TRUE., eos='PR') + + call general_result( fluid_to_characterize) write(*, *) fluid_to_characterize - print*, fluid_to_characterize%C, fluid_to_characterize%a, fluid_to_characterize%b, & - fluid_to_characterize%plus_mw, fluid_to_characterize%plus_z - print*, fluid_to_characterize%c_max - end program main diff --git a/src/bfr_routines/bfr_routines.f90 b/src/bfr_routines/bfr_routines.f90 index 5010f9d..3f41a50 100644 --- a/src/bfr_routines/bfr_routines.f90 +++ b/src/bfr_routines/bfr_routines.f90 @@ -1,224 +1,293 @@ module bfr_routines + !! This module calculates the best possible linear regression for the + !! distribution of mole fractions as a function of carbon number. use constants contains subroutine Linear_Regression(x, y, a, b, r2) - !! This subroutine computes the regression line for a data set of x, y variables. + !! This subroutine computes the linear regression parameters + !! (slope, intercept, and R²) for a given dataset of independent (x) + !! and dependent (y) variables. implicit none + ! Variable declarations + integer :: i ! Loop index + integer :: n ! Total number of data points + real(pr), intent(in) :: x(:) + !! Input array containing independent variable values (size n) + real(pr), intent(in) :: y(:) + !! Input array containing dependent variable values (size n) + real(pr), intent(out) :: a + !! Output variable: Slope of the regression line + real(pr), intent(out) :: b + !! Output variable: Intercept of the regression line + real(pr), intent(out) :: r2 + !! Output variable: Coefficient of determination (R²) + + ! Internal variables for summations + real(pr) :: t1, t2, t3, t4 + ! Summation terms for slope and intercept calculation + real(pr) :: aux1, aux2, aux3, aux4, aux5 + ! Auxiliary variables for R² calculation + + ! Get the number of data points + n = size(x) - integer :: i - integer :: n !! total number of data set. - real(pr), intent(in) :: x(:) !! x: input array of length n which contains the set of independent variable - real(pr), intent(in) :: y(:) !! y: input array of length n which contains the set of dependent variable - real(pr) :: t1,t2,t3,t4, aux1,aux2,aux3,aux4,aux5 ! internal variables - real(pr), intent(out) :: a !! a: output real variable. Slope of the regression line - real(pr), intent(out) :: b !! b: output real variable. Intercept of the regression line - real(pr), intent(out) :: r2 !! r2: output real variable. Square correlation coefficient + ! Initialize summation variables + t1 = 0.0; t2 = 0.0; t3 = 0.0; t4 = 0.0 - n = size(x) - !Calculation of a y b --> y=a*x+b - t1=0.;t2=0.;t3=0.;t4=0.; - - do i=1,n - t1=t1+x(i)*y(i) - t2=t2+x(i) - t3=t3+y(i) - t4=t4+x(i)**2 + ! Compute summations needed for slope (a) and intercept (b) + do i = 1, n + t1 = t1 + x(i) * y(i) + t2 = t2 + x(i) + t3 = t3 + y(i) + t4 = t4 + x(i) ** 2 end do - a=(n*t1-t2*t3)/(n*t4-t2**2) - b=(t3-a*t2)/n - !coefficient calculation of correlations r2 - aux1=0.;aux2=0.;aux3=0.;aux4=0.;aux5=0.; - - do i=1, n - aux1= aux1 + x(i)*y(i) - aux2= aux2 + x(i) - aux3= aux3 + y(i) - aux4= aux4 + x(i)**2 - aux5= aux5 + y(i)**2 + ! Compute slope (a) and intercept (b) using least squares formula + a = (n * t1 - t2 * t3) / (n * t4 - t2 ** 2) + b = (t3 - a * t2) / n + + ! Initialize auxiliary variables for correlation coefficient (R²) + aux1 = 0.0; aux2 = 0.0; aux3 = 0.0; aux4 = 0.0; aux5 = 0.0 + + ! Compute summations needed for R² calculation + do i = 1, n + aux1 = aux1 + x(i) * y(i) + aux2 = aux2 + x(i) + aux3 = aux3 + y(i) + aux4 = aux4 + x(i) ** 2 + aux5 = aux5 + y(i) ** 2 end do - r2=(aux1-aux2*aux3/n)**2 /((aux4-aux2**2/n)*(aux5-aux3**2/n)) + ! Compute coefficient of determination (R²) + r2 = (aux1 - aux2 * aux3 / n) ** 2 / & + ((aux4 - aux2 ** 2 / n) * (aux5 - aux3 ** 2 / n)) end subroutine Linear_Regression subroutine Best_Linear_Regression(scn_nc, scn, scn_z, plus_z, a, b, r2, & n_init, c_max_blr) - !! This subroutine calculates the best regression line for an fluid. + !! This subroutine calculates the best linear regression for a + !! fluid dataset. It determines the optimal regression line parameters + !! and the carbon number at which plus_z is reached. + implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_blr !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), intent(in) :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a !! output real variable. Slope of the best regression line. - real(pr), intent(out) :: b !! output real variable. Intercept of the best regression line. - real(pr), intent(out) :: r2 !! output real variable. Square correlation coefficient. - integer, intent(out) :: n_init !! minimum carbon number obtained from the best linear regression - integer :: i, j, k , k_old, n_best, x_aux ! internal variables - real(pr), dimension(scn_nc) :: x_blr, y_blr ! x_blr: carbon number vector, y_blr: logarithm of mole fraction vector; for each step - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux ! auxiliary variables - - k=5 - r2=0.0001_pr - r2_old=0.00001_pr - r2_best=0.0001_pr - - do while (r2.gt.r2_old.or.r2_old.lt.0.9) - k_old=k - r2_old=r2 - a_old=a - b_old=b - - if (r2.gt.r2_best) then - r2_best=r2 - a_best=a - b_best=b - n_best=scn(scn_nc-k+2) + ! Input variables + integer, intent(in) :: scn_nc + !! Total number of single cuts considered in the fluid + integer, intent(in) :: scn(:) + !! Array of single cuts considered in the fluid + real(pr), intent(in) :: scn_z(:) + !! Corresponding mole fractions of SCN cuts + real(pr), intent(in) :: plus_z + !! Composition of the residual fraction from input file + + ! Output variables + real(pr), intent(out) :: a !! Slope of the best regression line + real(pr), intent(out) :: b !! Intercept of the best regression line + real(pr), intent(out) :: r2 !! Coefficient of determination (R²) + integer, intent(out) :: n_init + !! Minimum carbon number obtained from the best linear regression + integer, intent(out) :: c_max_blr + !! CN at which plus_z is reached, computed from the best linear + !! distribution (BLR) + + ! Internal variables + integer :: i, j, k, k_old, n_best, x_aux + real(pr), dimension(scn_nc) :: x_blr, y_blr + !! x_blr: Carbon number vector, y_blr: Logarithm of mole fraction vector + real(pr) :: r2_old, r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux + !! Auxiliary variables + + ! Initialization + k = 5 ! This value sets the carbon number from which the regressions start + r2 = 0.0001_pr + r2_old = 0.00001_pr + r2_best = 0.0001_pr + + ! Iterative regression refinement + do while (r2 > r2_old .or. r2_old < 0.9) + k_old = k + r2_old = r2 + a_old = a + b_old = b + + if (r2 > r2_best) then + r2_best = r2 + a_best = a + b_best = b + n_best = scn(scn_nc - k + 2) end if - if (k.gt.scn_nc) then - if (r2.gt.r2_best)then - n_init=scn(scn_nc-k+2) + if (k > scn_nc) then + if (r2 > r2_best) then + n_init = scn(scn_nc - k + 2) go to 22 else - r2=r2_best - a=a_best - b=b_best - n_init=n_best + r2 = r2_best + a = a_best + b = b_best + n_init = n_best go to 22 end if end if - j=1 - x_blr=0. - y_blr=0. + j = 1 + x_blr = 0.0 + y_blr = 0.0 - do i=scn_nc-k+1, scn_nc - x_blr(j)=scn(i) - y_blr(j)=scn_z(i) - j=j+1 + do i = scn_nc - k + 1, scn_nc + x_blr(j) = scn(i) + y_blr(j) = scn_z(i) + j = j + 1 end do + ! Perform linear regression on selected data call Linear_Regression(x_blr(:k), y_blr(:k), a, b, r2) - k=k+1 + k = k + 1 end do - r2=r2_old - a=a_old - b=b_old - n_init=scn(scn_nc-k_old+2) + ! Store best values + r2 = r2_old + a = a_old + b = b_old + n_init = scn(scn_nc - k_old + 2) 22 continue - z_sum = 0d0 - x_aux=scn(scn_nc) + ! Compute CN at which plus_z is reached + z_sum = 0.0d0 + x_aux = scn(scn_nc) - do while (z_sum.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a*x_aux+b) - z_sum=z_sum+z_aux + do while (z_sum < plus_z .and. x_aux < 300) + x_aux = x_aux + 1 + z_aux = exp(a * x_aux + b) + z_sum = z_sum + z_aux end do - c_max_blr=x_aux - - !print*, a, b, c_max_blr, r2, n_init + c_max_blr = x_aux end subroutine Best_Linear_Regression subroutine LimitLine(scn_nc, scn, plus_z, a_blr, b_blr, a_lim, b_lim, & c_max_lim, half) - !! This subroutine obtains the limit line constants for an fluid. + !! This subroutine obtains the limit line constants for a fluid dataset. + !! It determines the parameters of the limit line and the carbon number + !! at which plus_z is reached. implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid. - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid. - integer, intent(out) :: c_max_lim !! maximum carbon number obtained for limit line. - !! output CN at which plus_z is reached, as the summation of single z(i) from the limit distribution (blr). - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression. - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file. - real(pr), intent(out) :: a_lim !! output real variable. Slope of the limit line. - real(pr), intent(out) :: b_lim !! output real variable. Intercept of the limit line. - real(pr), intent(out) :: half ! variable defined in the numerical method to converge the area under the curve of the function. + ! Input variables + integer, intent(in) :: scn_nc + !! Total number of single cuts considered in the fluid + integer, intent(in) :: scn(:) + !! Array of single cuts considered in the fluid + real(pr), intent(in) :: a_blr + !! Slope from the best linear regression + real(pr), intent(in) :: b_blr + !! Intercept from the best linear regression + real(pr), intent(in) :: plus_z + !! Composition of the residual fraction from input file + + ! Output variables + real(pr), intent(out) :: a_lim !! Slope of the limit line + real(pr), intent(out) :: b_lim !! Intercept of the limit line + integer, intent(out) :: c_max_lim + !! Maximum carbon number obtained for the limit line + real(pr), intent(out) :: half + !! Parameter used in the numerical method to converge the area under the curve + + ! Internal variables real(pr) :: z_lim, cross_cn, z_cross, z_aux - integer :: x_aux ! internal variable + integer :: x_aux - !A and B limit calculation. + ! A and B limit calculation z_lim = 0.0_pr - half = 0.506_pr ! modified by oscar, the variable half was removed from the argument + half = 0.506_pr !! Initial half-value - do while (z_lim.lt.plus_z) + do while (z_lim < plus_z) half = half + 0.002d0 - cross_cn = scn(scn_nc)+half ! typically 19.508 in first try - z_cross = exp((a_blr*cross_cn) + b_blr) - a_lim = -(z_cross)/plus_z - b_lim = log(z_cross)-(a_lim*cross_cn) - x_aux=scn(scn_nc) + cross_cn = scn(scn_nc) + half !! Typically 19.508 in first iteration + z_cross = exp((a_blr * cross_cn) + b_blr) + a_lim = -(z_cross) / plus_z + b_lim = log(z_cross) - (a_lim * cross_cn) + x_aux = scn(scn_nc) z_lim = 0._pr - do while (z_lim.lt.plus_z.and.x_aux<300) - x_aux=x_aux+1 - z_aux= exp(a_lim*x_aux+b_lim) - z_lim=z_lim+z_aux + + do while (z_lim < plus_z .and. x_aux < 300) + x_aux = x_aux + 1 + z_aux = exp(a_lim * x_aux + b_lim) + z_lim = z_lim + z_aux end do - continue end do c_max_lim = x_aux end subroutine LimitLine - subroutine Line_C60_max (scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & + subroutine Line_C60_max(scn_nc, scn, plus_z, a_blr, b_blr, half, a_60, & b_60, c_max_60) - !! This subroutine obtains the Cmax60 line constants for an fluid. - + !! This subroutine obtains the Cmax60 line constants for a fluid. implicit none - integer, intent(in) :: scn_nc !! integer input variable set to the total number of single cuts being considered in the fluid - integer, intent(in) :: scn(:) !! set of singles cuts being considered in the fluid - integer, intent(out) :: c_max_60 !!output CN at which Zp is reached, as the summation of single z(i) from the Cmax60 distribution. - real(pr), intent(in) :: a_blr !! input constant from the Best linear regression - real(pr), intent(in) :: b_blr !! input constant from the Best linear regression - real(pr), intent(in) :: plus_z !! composition of residual fraction from input file - real(pr), intent(out) :: a_60 !! output real variable. Slope of the limit line - real(pr), intent(out) :: b_60 !! output real variable. Intercept of the C60max line. - real(pr), intent(in) :: half !variable defined in the numerical method to converge the area under the curve of the function - integer :: x_aux ! auxiliary variable - real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range ! internal variables - real(pr) :: a_old, Zp_60, F, dF_dA ! internal variables - - cross_cn = scn(scn_nc)+half ! typically 19.51 - z_cross= exp((a_blr*cross_cn)+b_blr) - F_tol=1_pr - a_tol=1_pr - a_60=a_blr + + ! Input variables + integer, intent(in) :: scn_nc + !! Total number of single cuts considered in the fluid + integer, intent(in) :: scn(:) + !! Array of single cuts considered in the fluid + real(pr), intent(in) :: plus_z + !! Composition of residual fraction from input file + real(pr), intent(in) :: a_blr + !! Slope from the best linear regression + real(pr), intent(in) :: b_blr + !! Intercept from the best linear regression + real(pr), intent(in) :: half + !! Numerical method variable to adjust the integration range + + ! Output variables + real(pr), intent(out) :: a_60 !! Slope of the Cmax60 line + real(pr), intent(out) :: b_60 !! Intercept of the Cmax60 line + integer, intent(out) :: c_max_60 + !! Carbon number at which plus_z is reached from the Cmax60 distribution + + ! Internal variables + integer :: x_aux + real(pr) :: z_sum, cross_cn, z_cross, z_aux, F_tol, a_tol, var_range + real(pr) :: a_old, Zp_60, F, dF_dA + + ! Compute initial values + cross_cn = scn(scn_nc) + half + z_cross = exp((a_blr * cross_cn) + b_blr) + F_tol = 1_pr + a_tol = 1_pr + a_60 = a_blr var_range = 60.5_pr - cross_cn - do while (a_tol.gt.1e-6_pr.or.F_tol.gt.1e-6_pr) - a_old=a_60 - Zp_60 = (exp(a_60*var_range+log(z_cross))-z_cross)/a_60 + ! Newton-Raphson iteration for a_60 + do while (a_tol > 1e-6_pr .or. F_tol > 1e-6_pr) + a_old = a_60 + Zp_60 = (exp(a_60 * var_range + log(z_cross)) - z_cross) / a_60 F = plus_z - Zp_60 - dF_dA = - (var_range*(a_60*Zp_60+z_cross)-Zp_60) / a_60 ! modified by oscar, error in the derivative in previous version is corrected - a_60=a_old-(F/dF_dA) - a_tol=abs(a_old-a_60) - F_tol=abs(F) + dF_dA = - (var_range * (a_60 * Zp_60 + z_cross) - Zp_60) / a_60 + a_60 = a_old - (F / dF_dA) + a_tol = abs(a_old - a_60) + F_tol = abs(F) end do - b_60=log(z_cross)-a_60*cross_cn - z_sum=0.0_pr - x_aux=scn(scn_nc) + b_60 = log(z_cross) - a_60 * cross_cn + z_sum = 0.0_pr + x_aux = scn(scn_nc) - do while (z_sum.lt.plus_z) - x_aux=x_aux+1 - z_aux= exp(a_60*x_aux+b_60) - z_sum=z_sum+z_aux + ! Compute c_max_60 + do while (z_sum < plus_z) + x_aux = x_aux + 1 + z_aux = exp(a_60 * x_aux + b_60) + z_sum = z_sum + z_aux end do - c_max_60=x_aux + c_max_60 = x_aux - !print*, a_60,b_60,c_max_60 end subroutine Line_C60_max diff --git a/src/characterization/characterization.f90 b/src/characterization/characterization.f90 index 9d42ae0..2839c67 100644 --- a/src/characterization/characterization.f90 +++ b/src/characterization/characterization.f90 @@ -1,58 +1,136 @@ module characterization - - use constants - use dtypes, only: FluidData, FluidDataOut - use data_from_input, only: data_from_file, FluidData - use molecular_weigth - use density - use lumping - use critical_parameters - - + !! This module provides functions to characterize a fluid using input data. + !! The characterization includes: + !! - Molecular weight computation + !! - Density estimation using various methods + !! - Lumping of components for compositional modeling + !! - Calculation of critical parameters (Tc, Pc, ω) + !! The module relies on external modules for data handling + !!and property calculations. + + use constants + !! Module containing physical constants + use dtypes, only: FluidData, FluidDataOut + !! Data structures for fluid properties + use data_from_input, only: data_from_file, FluidData + !! Function to read fluid data from a file + use molecular_weigth + !! Module for molecular weight calculations + use density + !! Module for density calculations + use lumping + !! Module for lumping methodologies + use critical_parameters + !! Module for critical parameter estimation + contains - type(FluidDataOut) function characterize(file, mw_source, method, pho_method, fix_C, eos) & - result(characterization) + type(FluidDataOut) function characterize(file, mw_source, method, & + rho_method, fix_C, eos) result(characterization) + !! This function performs the full characterization of a reservoir + !! fluid based on input data. The characterization process includes + !! molecular weight assignment, density computation, lumping of + !! components, and estimation of critical parameters. + !! The function reads fluid composition from a file and applies different + !! calculation methodologies depending on the specified input parameters. - type(FluidData) :: fluid - character(len=*), intent(in) :: file !! file name - character(len=*), intent(in), optional :: method !! plus_mw or global_mw + !input variables + character(len=*), intent(in) :: file + !! Name of the input file containing the fluid composition data. character(len=*), intent(in) :: mw_source - logical, intent(in), optional :: fix_C + !! Determines whether molecular weights from predefined SCNs in the input + !! file are used or whether correlations are used to estimate these values. + !! Indicates data source ("experimental" or "calculated") + character(len=*), intent(in), optional :: method + !! Specifies the method for molecular weight calculation. + !! Available options: + !! - "plus_mw": reproduces the molecular weight of the residual fraction. + !! - "global_mw": reproduces the global molecular weight of fluid. character(len=*), intent(in) :: eos - integer, intent(in) :: pho_method !! selector method of density funtion. - - - - fluid = data_from_file(file) - characterization%input_data = fluid - allocate(characterization%scn_z(fluid%scn_nc)) - allocate(characterization%log_scn_z(fluid%scn_nc)) - allocate(characterization%scn_mw(fluid%scn_nc)) - allocate(characterization%carbon_number_plus(0)) - allocate(characterization%plus_z_i(0)) - allocate(characterization%product_z_mw_plus_i(0)) - allocate(characterization%scn_zm(fluid%scn_nc)) - allocate(characterization%scn_i(0)) - allocate(characterization%plus6_density(0)) - allocate(characterization%mol_fraction(0)) - allocate(characterization%lumped_z(0)) - allocate(characterization%lumped_mw(0)) - allocate(characterization%lumped_densities(0)) - allocate(characterization%critical_temperature(0)) - allocate(characterization%critical_pressure(0)) - allocate(characterization%acentric_factor(0)) - allocate(characterization%m_funtion(0)) - - - call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, fix_C=fix_C, characterization= characterization) - call density_funtion(fluid=fluid, mw_source=mw_source, pho_method= pho_method, characterization=characterization) + !! Specifies the equation of state (EOS) used for + !! calculating critical parameters. + logical, intent(in), optional :: fix_C + !! If present and set to .TRUE., it fixes the value of the constant 'C' + !! instead of estimating it dynamically. + integer, intent(in) :: rho_method + !! Selector for the density computation method. + !! Different methods can be used for different fluid + !! characterization approaches. + ! + !! method "1":use experimental scn's densities reported to calculate + !! experimental volume for 6plus. The densities values are calculated + !! since CN plus + ! + !! method "2": This densities doesn't experimentals, they can be + !! reference values. the densities values are calculated since 6 plus + !! by correlations. + ! + !! method "3": use experimental value of density of C6+ to calculate + !! the experimental volume + + ! Internal variables + type(FluidData) :: fluid + !! Structure to store the fluid composition and properties read + !! from the input file. + + ! Read input fluid data + fluid = data_from_file(file) + !! Read the fluid composition from the input file + characterization%input_data = fluid + !! Store input data in the output structure + + ! Dynamic memory allocation + allocate(characterization%scn_z(fluid%scn_nc)) + ! Mole fractions of single carbon numbers (SCNs) + allocate(characterization%log_scn_z(fluid%scn_nc)) + ! Logarithm of SCN mole fractions + allocate(characterization%scn_mw(fluid%scn_nc)) + ! Molecular weights of SCNs + allocate(characterization%carbon_number_plus(0)) + ! Carbon number for the heavy fraction + allocate(characterization%plus_z_i(0)) + ! Mole fraction of the heavy fraction + allocate(characterization%product_z_mw_plus_i(0)) + ! Product of mole fraction and molecular weight for the heavy fraction + allocate(characterization%scn_zm(fluid%scn_nc)) + ! Mole fraction adjusted for molecular weight calculations + allocate(characterization%scn_i(0)) + ! Indexing for SCN components + allocate(characterization%plus6_density(0)) + ! Density of the heavy fraction + allocate(characterization%mol_fraction(0)) + ! Mole fractions after characterization + allocate(characterization%lumped_z(0)) + ! Lumped mole fractions + allocate(characterization%lumped_mw(0)) + ! Lumped molecular weights + allocate(characterization%lumped_densities(0)) + ! Lumped component densities + allocate(characterization%critical_temperature(0)) + ! Estimated critical temperature (Tc) + allocate(characterization%critical_pressure(0)) + ! Estimated critical pressure (Pc) + allocate(characterization%acentric_factor(0)) + ! Acentric factor (ω) + allocate(characterization%m_funtion(0)) + ! Function for molecular weight adjustment + + + ! Computation steps + ! Step 1: Compute molecular weight for the heavy fraction + call get_c_or_m_plus(fluid=fluid, mw_source=mw_source, method=method, & + fix_C=fix_C, characterization=characterization) + + ! Step 2: Compute density using the selected method + call density_funtion(fluid=fluid, rho_method=rho_method, & + characterization=characterization) + + ! Step 3: Apply lumping methodology to reduce component count call lump(fluid=fluid, characterization=characterization) - call get_critical_constants(fluid= fluid, characterization=characterization, eos=eos) + !! Step 4: Compute critical parameters (Tc, Pc, ω) based on EOS + call get_critical_constants( characterization=characterization, eos=eos) end function characterize end module characterization - - diff --git a/src/constans.f90 b/src/constans.f90 index 47ee4a7..4c4b608 100644 --- a/src/constans.f90 +++ b/src/constans.f90 @@ -1,4 +1,5 @@ module constants + !! This module set float precision use iso_fortran_env, only: real64 implicit none integer, parameter :: pr=real64 diff --git a/src/critical_parameters/critical_parameters.f90 b/src/critical_parameters/critical_parameters.f90 index 39bfcf4..d3ab19a 100644 --- a/src/critical_parameters/critical_parameters.f90 +++ b/src/critical_parameters/critical_parameters.f90 @@ -1,60 +1,105 @@ module critical_parameters - - use constants - use dtypes + !! This module is responsible for calculating the critical + !! parameters of a fluid, namely: + !! - Critical Temperature (Tc) + !! - Critical Pressure (Pc) + !! - Acentric Factor (omega) + !! + !! The calculations are based on Pedersen correlations + !! (adapted from the code "C7plusPedersenTcPcOm") + !! and on the evaluation of EOS parameters taken from CubicParam. + + use constants ! Contains physical and mathematical constants + use dtypes !Defines data types, including FluidDataOut use defined_critical_parameters + ! Defines critical parameters and related constants contains - subroutine get_critical_constants(fluid, characterization, eos) - !! this subroutine doing: - !! - Calculation of Tc, Pc, omega from Pedersen correlations (adapted from code "C7plusPedersenTcPcOm") - !! - Calculation of EoS parameters (taken from CubicParam) - !use critical_parameters + subroutine get_critical_constants(characterization, eos) + !! This subroutine calculates the critical parameters of a fluid, namely: + !! - Critical Temperature (Tc) + !! - Critical Pressure (Pc) + !! - Acentric Factor (omega) + !! + !! Additionally, it computes a modification function (m_funtion) + !! used in the correlations. + !! + !! The equations used in the calculations are: + !! tc = c1 * rho + c2 * log(mw) + c3 * mw + (c4 / mw) + !! pc = exp(d1_p + d2 * (rho ** d5) + (d3 / mw) + (d4 / (mw ** 2))) + !! om = (-b + sqrt(b ** 2 - 4 * a * modified_c)) / (2 * a) + !! where: + !! modified_c = c0 - m_funtion + !! + !! Depending on the EOS specified (e.g., "SRK", "PR", or "RKPR"), + !! an appropriate expression for m_funtion is selected. + !! + !! Input: + !! characterization : + !! FluidDataOut type (in/out) containing fluid information. + !! It will be updated with the calculated critical parameters. + !! + !! eos : + !! Character string specifying the equation of state to be used. + !! + !! Internal variables: + !! Allocatable arrays for storing intermediate values for Tc, Pc, omega, + !! and m_funtion. implicit none - type(FluidData), intent(in) :: fluid type(FluidDataOut), intent(inout) :: characterization character(len=*), intent(in) :: eos + ! Allocatable arrays for intermediate calculations real(pr), allocatable :: modified_c(:) - real(pr), allocatable :: tc(:) - real(pr), allocatable :: pc(:) - real(pr), allocatable :: om(:) + !! Array for the modified constant: modified_c = c0 - m_funtion + real(pr), allocatable :: tc(:) !! Array for the critical temperature (Tc) + real(pr), allocatable :: pc(:) !! Array for the critical pressure (Pc) + real(pr), allocatable :: om(:) !! Array for the acentric factor (omega) real(pr), allocatable :: m_funtion(:) - integer :: i + !! Array for the m function used in the correlations associate(& rho => characterization%lumped_densities, & + ! rho: Lumped densities of the fluid mw => characterization%lumped_mw & + ! mw: Lumped molecular weights of the fluid ) - + ! Allocate memory for the intermediate arrays allocate(modified_c(0)) allocate(tc(0)) allocate(pc(0)) allocate(om(0)) allocate(m_funtion(0)) + ! Call the subroutine that sets parameters for critical calculations + ! based on the EOS call get_parameteres_for_critical(eos) - tc = c1*rho + c2*log(mw) + c3*mw + (c4/mw) - pc = exp(d1_p + d2*(rho**(d5)) + (d3/mw) + (d4/(mw**2))) + ! Calculate the critical temperature (Tc) using the Pedersen correlation: + tc = c1 * rho + c2 * log(mw) + c3 * mw + (c4 / mw) + + ! Calculate the critical pressure (Pc) using an exponential correlation: + pc = exp(d1_p + d2 * (rho ** d5) + (d3 / mw) + (d4 / (mw ** 2))) + !Determine the modification function (m_funtion) based on the specified EOS if (eos == "SRK") then - m_funtion = e1 + e2*mw + e3*rho + e4*(mw**(2)) + ! For the Soave-Redlich-Kwong (SRK) EOS: + m_funtion = e1 + e2 * mw + e3 * rho + e4 * (mw ** 2) else if (eos == "PR" .or. eos =="RKPR") then - m_funtion = e3*rho + 0.2 + 1.8*(1-exp(-mw/220)) + ! For the Peng-Robinson (PR) or RKPR EOS: + m_funtion = e3 * rho + 0.2 + 1.8 * (1 - exp(-mw / 220)) ! Alternative expression without max with MW - May 2019 endif endif + ! Calculate the modified constant: modified_c = c0 - m_funtion - om = (-b + sqrt(b**2 - 4*a*modified_c))/(2*a) - - !do i = 1, size(mw) - ! print*, mw(i), rho(i), tc(i), pc(i), om(i), m_funtion(i) - !end do - + ! Calculate the acentric factor (omega) using the quadratic formula: + om = (-b + sqrt(b ** 2 - 4 * a * modified_c))/(2 * a) + ! Update the characterization structure with the calculated + ! critical parameters. characterization%critical_temperature = tc characterization%critical_pressure = pc characterization%acentric_factor = om diff --git a/src/critical_parameters/defined_critical_param.f90 b/src/critical_parameters/defined_critical_param.f90 index b00e1d1..1d4941d 100644 --- a/src/critical_parameters/defined_critical_param.f90 +++ b/src/critical_parameters/defined_critical_param.f90 @@ -1,25 +1,69 @@ module defined_critical_parameters + !! This module defines global arrays and coefficients used for calculating + !! the critical properties of a fluid. These properties include the critical + !! temperature (Tc), critical pressure (Pc), and acentric factor (omega). + !! The module provides default values for these properties as well as + !! correlation coefficients used in Pedersen correlations and + !! EOS parameter calculations. use constants implicit none - real(pr), allocatable :: tc_def(:), pc_def(:), om_def(:) - real(pr):: c1, c2, c3, c4, d1_p, d2, d3, d4, d5, e1, e2, e3, e4 - real(pr):: a, b, c0, del_1 + ! Allocatable arrays for default critical properties for various components: + real(pr), allocatable :: tc_def(:) !! Define critical temperatures [K] + real(pr), allocatable :: pc_def(:) !! Define critical pressures [bar] + real(pr), allocatable :: om_def(:) !! Define acentric factors + + ! Global coefficients for critical property correlations: + ! Coefficients for critical temperature (Tc) correlation: + real(pr) :: c1, c2, c3, c4 + ! Coefficients for critical pressure (Pc) correlation: + real(pr) :: d1_p, d2, d3, d4, d5 + ! Coefficients for the modification function (m_funtion) used in + ! the correlations: + real(pr) :: e1, e2, e3, e4 + ! Additional coefficients for acentric factor (omega) calculation via a + ! quadratic expression: + real(pr) :: a, b, c0, del_1 + ! 'del_1' is an additional parameter (currently not used). contains subroutine get_parameteres_for_critical(eos) + !! This subroutine initializes the define componets critical property + !! arrays (tc_def, pc_def,om_def) and sets the global correlation + !! coefficients based on the specified equation of state (EOS). + !! + !! The default critical arrays are defined element-wise for different + !! components. + !! + !! Depending on the input string 'eos', the subroutine assigns + !! different sets of coefficients: + !! - "SRK" : Soave-Redlich-Kwong EOS coefficients. + !! - "PR" : Peng-Robinson EOS coefficients (original Pedersen values). + !! - "RKPR" : Modified Peng-Robinson (RKPR) EOS coefficients. + !! + !! Input: + !! eos : A character string specifying the equation of state to be used. + !! Supported values: "SRK", "PR", "RKPR". + !! + !! The coefficients set in this subroutine are used by other modules to + !! calculate the fluid's critical temperature (Tc), pressure (Pc), and + !! acentric factor (omega). + implicit none character(len=*) :: eos + ! Define default critical properties for various components: tc_def = [126.2, 304.21, 190.564, 305.32, 369.83, 408.14, 425.12, 460.43, 469.7] pc_def = [34.0, 73.83, 45.99, 48.72, 42.48, 36.48, 37.96, 33.81, 33.7] om_def = [0.038, 0.224, 0.012, 0.099, 0.152, 0.181, 0.20, 0.228, 0.252] + ! Set correlation coefficients based on the specified EOS. select case(eos) case("SRK") + ! For the Soave-Redlich-Kwong (SRK) EOS: c1 = 1.6312d2 c2 = 8.6052d1 c3 = 4.3475d-1 diff --git a/src/density/density_funtion.f90 b/src/density/density_funtion.f90 index 32eb4b1..8ac78ce 100644 --- a/src/density/density_funtion.f90 +++ b/src/density/density_funtion.f90 @@ -1,154 +1,223 @@ module density + !! This module is responsible for calculating the density function for a + !! fluid. It contains routines to: + !! - Compute the experimental volume using density data (density_method). + !! - Iteratively adjust the density function parameters so that the + !! calculated volume matches the experimental volume (density_funtion). + !! - Calculate the volume based on the current density function parameters + !! (calculate_volume) + use constants use dtypes, only: FluidData, FluidDataOut contains - subroutine density_funtion(fluid, mw_source, pho_method, characterization) - !! this subroutine ... + subroutine density_funtion(fluid, rho_method, characterization) + !! This subroutine adjust the density function constants (a_d and b_d) + !! until the calculated volume matches the experimental volume. This is + !! done using an iterative (secant-like) approach. + !! + !! Process: + !! 1. Call density_method to compute the experimental volume + !! (volume_exp) based on the fluid data. + !! 2. Initialize the density function constants a_d and b_d with an + !! initial guess. + !! 3. Call calculate_volume to compute the calculated volume + !! (volume_cal) using the current density function parameters. + !! 4. Compute the difference between volume_cal and volume_exp. + !! 5. Iteratively adjust a_d (and update b_d accordingly) until the + !! absolute difference is less than the specified tolerance (0.001). + implicit none - type(FluidData) :: fluid + type(FluidData) :: fluid !! (FluidData) Contains fluid properties. type(FluidDataOut) :: characterization - character(len=*), intent(in) :: mw_source - integer, intent(in) :: pho_method - !real(pr), allocatable :: volume_6plus_cal(:) - real(pr) :: volume_exp + !! Structure to store density function constants and volumes. + integer, intent(in) :: rho_method + !! Selects the density calculation method real(pr) :: a_d_old, aux real(pr) :: difference, difference_old - !volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & - ! (fluid%product_z_mw_plus/fluid%plus_density) - - - - - - associate( & a_d => characterization%a_d , & + ! Slope constant of the density function b_d => characterization%b_d, & + ! Intercept constant of the density function volume_cal => characterization%volume_cal, & + ! Calculated volume from the density function. volume_exp => characterization%volume_exp & + ! Experimental volume obtained from the fluid data ) - - call density_method (fluid, pho_method,characterization) + + ! Compute the experimental volume using the selected density method. + call density_method (fluid, rho_method,characterization) ! Now find constants for the Density function - a_d = -0.50_pr ! initial guess - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume(fluid, mw_source, pho_method, characterization) + ! Initialize the density function constants with initial guesses + a_d = -0.50_pr ! Initial guess for a_d. + b_d = 0.685_pr - a_d * exp(-0.60_pr) ! Initial guess for b_d + + ! Calculate the volume using the current density function. + call calculate_volume(fluid, rho_method, characterization) difference_old = volume_cal - volume_exp a_d_old = a_d + ! Update the initial guess slightly. a_d = -0.49_pr - b_d = 0.685_pr - a_d*exp(-0.60_pr) - call calculate_volume(fluid, mw_source, pho_method, characterization) + b_d = 0.685_pr - a_d * exp(-0.60_pr) + call calculate_volume(fluid, rho_method, characterization) difference = volume_cal - volume_exp - + + ! Iteratively adjust a_d (and update b_d) until the difference is + ! within tolerance. do while (abs(difference) > 0.001_pr) aux = a_d - a_d = a_d - difference*(a_d - a_d_old)/(difference - difference_old) - b_d = 0.685_pr - a_d*exp(-0.60_pr) + a_d = a_d - difference * (a_d - a_d_old) / (difference - difference_old) + b_d = 0.685_pr - a_d * exp(-0.60_pr) a_d_old = aux difference_old = difference - call calculate_volume(fluid, mw_source, pho_method, characterization) + call calculate_volume(fluid, rho_method, characterization) difference = volume_cal - volume_exp end do end associate end subroutine density_funtion - subroutine density_method (fluid, pho_method,characterization) - !! this subroutine ... + subroutine density_method (fluid, rho_method,characterization) + !! This subroutine compute the experimental volume (volume_exp) using the + !! fluid's density data. + !! + !! Details: + !! - For rho_method 1 and 2: + !! The experimental volume is calculated using SCN densities + !! provided in the fluid data. The volume for the 6+ fraction is + !! computed from the product of mole fraction and molecular weight + !! divided by the corresponding density. method 1 use experimental + !! densities for scn, while method 2 used reference values of + !! densities for scn cuts. + !! - For rho_method 3: + !! A single experimental value is used to compute the density + !! distribution. + implicit none type(FluidData) :: fluid + !! Contains fluid properties and density data type(FluidDataOut) :: characterization - integer, intent(in) :: pho_method + !! Structure to store density function constants and volumes. + integer, intent(in) :: rho_method + !! Selector for the density computation method real(pr) :: volume_exp + !! Experimental volume obtained from the fluid data - select case (pho_method) ! cambiar a rho + select case (rho_method) case (1) - ! case 1 use experimental scn's densities reported to calculate - ! experimental volume for 6plus. The densities values are calculated - ! since CN plus + ! Case 1: + ! Use experimental SCN densities reported to calculate the + ! experimental volume for the 6+ fraction. The densities are based + ! on values calculated from the carbon number plus fraction. volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & (fluid%product_z_mw_plus/fluid%plus_density) characterization%volume_exp = volume_exp case (2) - !This case use scn densities reported to calculate experimental volume - ! this densities doesn't experimental, they can be reference values. - ! the densities values are calculated since 6 plus - + ! Case 2: + ! Use reported (reference) SCN densities to calculate the + ! experimental volume. The approach is similar to case 1. volume_exp = sum((fluid%product_z_mw_scn)/(fluid%scn_density)) + & (fluid%product_z_mw_plus/fluid%plus_density) characterization%volume_exp = volume_exp case (3) - ! This case use one experimental value (6plus) to calculate the density - ! distribution - + ! Case 3: + ! Use a single experimental value (6plus) to compute the density + ! distribution. volume_exp = (fluid%plus6_z_exp*fluid%plus6_mw_exp)/(fluid%plus6_density_exp) characterization%volume_exp = volume_exp end select - - - - - - - end subroutine density_method - subroutine calculate_volume(fluid, mw_source, pho_method, characterization) - !! this subroutine ... + subroutine calculate_volume(fluid, rho_method, characterization) + !! This subroutine calculate the volume (volume_cal) based on the current + !! density function parameters. + !! + !! Process: + !! 1. Compute estimated densities for the SCN fraction (scn_density_cal) + !! and the plus fraction (plus_density_cal) using the density + !! function constants a_d and b_d. + !! 2. Depending on the selected rho_method, calculate the total volume + !! by summing contributions: + !! - The SCN fraction volume is calculated as the sum over + !! adjusted mole fractions (scn_zm) divided by the corresponding + !! density. + !! - The plus fraction volume is calculated as the sum over the + !! product of mole fraction and molecular weight for the plus + !! fraction divided by the computed plus density. + !! 3. Update the characterization structure with the calculated volume + !! and the density array for the 6+ fraction. + implicit none - type(FluidData), intent(in) :: fluid + type(FluidData), intent(in) :: fluid + !! Contains fluid composition and density data. type(FluidDataOut), intent(inout) :: characterization - character(len=*), intent(in) :: mw_source - integer, intent(in) :: pho_method - real(pr) :: volume_cal - real(pr), allocatable :: plus_density_cal(:), scn_density_cal(:), plus6_density(:) - - !allocate(density_plus_cal(0), density_scn_cal(0), plus6_density(0)) - + !! Structure to update with calculated volume and density arrays + integer, intent(in) :: rho_method + !! Density computation method selector. + real(pr) :: volume_cal + ! Calculated volume based on the current density function + real(pr), allocatable :: plus_density_cal(:) + real(pr), allocatable :: plus6_density(:) + ! Density values for the 6+ fraction + real(pr), allocatable :: scn_density_cal(:) associate(& - a_d => characterization%a_d, & - b_d => characterization%b_d, & - cn_plus => characterization%carbon_number_plus, & - z_m_plus_i => characterization%product_z_mw_plus_i, & - scn_z_i => characterization%scn_z, & - scn_mw_i => characterization%scn_mw, & - scn_zm => characterization%scn_zm & + a_d => characterization%a_d, & ! Density function slope constant. + b_d => characterization%b_d, & ! Density function intercept constant. + cn_plus => characterization%carbon_number_plus, & + ! Carbon number for the plus fraction. + z_m_plus_i => characterization%product_z_mw_plus_i, & + ! Product of mole fraction and MW for the plus fraction. + scn_z_i => characterization%scn_z, & ! SCN mole fractions. + scn_mw_i => characterization%scn_mw, & ! SCN molecular weights. + scn_zm => characterization%scn_zm & + ! Adjusted SCN mole fractions for volume calculation. ) - - scn_density_cal = ((a_d)*(exp(- (real(fluid%scn, pr))/10._pr))) + b_d - plus_density_cal =((a_d)*(exp(- (real(cn_plus, pr) )/10._pr))) + b_d + + ! Compute estimated density for the SCN fraction using the density function: + scn_density_cal = ((a_d) * (exp(- (real(fluid%scn, pr))/10._pr))) + b_d + ! Compute estimated density for the plus fraction using the carbon number plus: + plus_density_cal =((a_d) * (exp(- (real(cn_plus, pr) )/10._pr))) + b_d + - select case (pho_method) + select case (rho_method) case(1) + ! Use the fluid's provided SCN densities and the computed + ! plus_density_cal. volume_cal = sum((scn_zm)/(fluid%scn_density)) + & sum(z_m_plus_i/plus_density_cal) plus6_density = [fluid%scn_density, plus_density_cal] case(2) + ! Use the computed SCN densities (scn_density_cal) and + ! plus_density_cal. volume_cal = sum((scn_zm)/(scn_density_cal)) + & sum(z_m_plus_i/plus_density_cal) plus6_density = [scn_density_cal, plus_density_cal] case(3) + ! Similar to case 2: use computed densities for both SCN and + ! plus fractions volume_cal = sum((scn_zm)/(scn_density_cal)) + & sum(z_m_plus_i/plus_density_cal) plus6_density = [scn_density_cal, plus_density_cal] end select + ! Update the characterization structure with the calculated volume and + ! density values characterization%volume_cal = volume_cal characterization%plus6_density = plus6_density @@ -156,5 +225,4 @@ subroutine calculate_volume(fluid, mw_source, pho_method, characterization) end subroutine calculate_volume - end module density diff --git a/src/input_data/data_input.f90 b/src/input_data/data_input.f90 index ed37eaa..1ed8a7a 100644 --- a/src/input_data/data_input.f90 +++ b/src/input_data/data_input.f90 @@ -3,7 +3,7 @@ module data_from_input !! an input file and organize it into a data structure. use dtypes, only: FluidData ! The FluidData type is imported, which stores all the fluid information. - use constants, only: pr + use constants, only: pr ! Pr precision is used to define floating point variables. implicit none @@ -36,7 +36,7 @@ end subroutine read_setup subroutine read_components(file, def_components, scn, scn_plus) !! This subroutine reads the component names from input file - character(len=*), intent(in) :: file + character(len=*), intent(in) :: file !! File name of fluid to characterize character(len=15), allocatable, intent(out) :: def_components(:) !! Names of defined components @@ -69,7 +69,7 @@ subroutine read_composition(file, def_comp_z, scn_z, plus_z, plus6_z_exp, & !! This subroutine reads the molar compositions of each component !! from input file - character(len=*), intent(in) :: file + character(len=*), intent(in) :: file !! File name of fluid to characterize real(pr), allocatable, intent(out) :: def_comp_z(:) !! Set of corresponding mole fractions of defined components @@ -134,7 +134,7 @@ subroutine read_molecular_weight(file, def_comp_mw, scn_mw, plus_mw, & integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified !! number of pseudos - integer :: numbers_ps + integer :: numbers_ps !! Number of pseudos in which the scn fractions grouped integer :: funit @@ -213,7 +213,7 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& !! molecular weight of residual fraction !! C30+ is considered the residual fraction real(pr), allocatable :: product_z_mw_i(:) - ! product of mole fraction and molecular weight of each + ! product of mole fraction and molecular weight of each ! component of the fluid real(pr), allocatable, intent(out) :: def_comp_w(:) !! mass fractions of defined components @@ -244,7 +244,8 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& product_z_mw_def_comp = (def_comp_z)*(def_comp_mw) product_z_mw_scn = (scn_z)*(scn_mw) product_z_mw_plus = (plus_z)*(plus_mw) - product_z_mw_i = [product_z_mw_def_comp, product_z_mw_scn, product_z_mw_plus] + product_z_mw_i = [product_z_mw_def_comp, product_z_mw_scn, & + product_z_mw_plus] sum_z_mw_i = sum(product_z_mw_i) w = (product_z_mw_i) / (sum_z_mw_i) def_comp_w = w(1:def_comp_nc) @@ -254,31 +255,34 @@ subroutine mass_fractions(file, w, product_z_mw_def_comp, product_z_mw_scn,& end subroutine mass_fractions subroutine read_density(file,scn_density,plus_density, plus6_density_exp, & - plus7_density_exp, plus12_density_exp, plus30_density_exp & - ) - !! This subroutine reads the density of each component from the input - !! file and calculated molar volume since C6 fraction. + plus7_density_exp, plus12_density_exp, plus30_density_exp) + !! Reads the density of each component from the input file and calculated + !! molar volume since C6 fraction. - character(len=*), intent(in) :: file !! file name - real(pr), allocatable, intent(out) :: scn_density(:) + character(len=*), intent(in) :: file !! File name + real(pr), allocatable, intent(out) :: scn_density(:) !! Set of corresponding densities of scn cuts - real(pr), intent(out):: plus_density + real(pr), intent(out):: plus_density !! Density of residual fraction from input file - real(pr), intent(out):: plus6_density_exp !! Experimental density for C6+ - real(pr), intent(out):: plus7_density_exp !! Experimental density for C7+ - real(pr), intent(out):: plus12_density_exp !! Experimental density for C12+ + real(pr), intent(out):: plus6_density_exp + !! Experimental density for C6+ + real(pr), intent(out):: plus7_density_exp + !! Experimental density for C7+ + real(pr), intent(out):: plus12_density_exp + !! Experimental density for C12+ real(pr), intent(out):: plus30_density_exp - integer :: def_comp_nc + !! Experimental density for C30+ + integer :: def_comp_nc !! Number of defined components being considered in the oil - integer :: scn_nc + integer :: scn_nc !! Number of single cuts being considered in the oil - integer :: scn_nc_ps - !! CN from which all SCN fractions will be lumped into the specified + integer :: scn_nc_ps + !! CN from which all SCN fractions will be lumped into the specified !! number of pseudos - integer :: numbers_ps - !! Number of pseudos in which the scn fractions grouped + integer :: numbers_ps + !! number of pseudos in which the scn fractions grouped integer :: funit - + namelist /nml_density/ scn_density, plus_density, plus6_density_exp, & plus7_density_exp, plus12_density_exp, plus30_density_exp ! Namelist for reading density values from the input file @@ -286,18 +290,19 @@ subroutine read_density(file,scn_density,plus_density, plus6_density_exp, & call read_setup(file, def_comp_nc, scn_nc, scn_nc_ps, numbers_ps) allocate(scn_density(scn_nc)) - open(newunit=funit, file=file) + open(newunit=funit, file=file) ! Open the input file and read densities using the namelist + read(funit, nml=nml_density) close(funit) end subroutine read_density type(FluidData) function data_from_file(file) result(data) - !! This funtion allows to obtain experimental data from data imput + !! This funtion allows to obtain experimental data from imput file + character(len=*), intent(in) :: file !! file name data%filename = file - call read_setup(file, data%def_comp_nc, data%scn_nc, data%scn_nc_ps, & data%numbers_ps & ) diff --git a/src/lumping/lumping_method.f90 b/src/lumping/lumping_method.f90 index e082987..357e4fe 100644 --- a/src/lumping/lumping_method.f90 +++ b/src/lumping/lumping_method.f90 @@ -4,55 +4,65 @@ module lumping contains subroutine lump(fluid, characterization) - !! This sobroutine ... + !! This subroutine performs lumping of fluid components based on their + !! molecular weight and composition to create pseudo-components for + !! reservoir simulations. implicit none - type(FluidData), intent(in) :: fluid + type(FluidData), intent(in) :: fluid !! Input fluid characterization data type(FluidDataOut), intent(inout) :: characterization - integer :: last_C ! CN of last single cut: typically 19 - integer :: i_last ! Number of elements in the distribution of CN+ - integer :: last ! Total elements starting after defined components. - integer :: scn_nc_input !! inicial number of single cuts being considered in the oil from data input - integer :: scn_nc_new !! new number of single cuts being considered in the oil defined as from scn_nc_ps variable. - real(pr) :: plus_w_new + !! Output lumped fluid characterization + integer :: last_C ! Carbon number of the last single cut, typically 19 + integer :: i_last ! Number of elements in the distribution of carbon numbers (CN+) + integer :: last + ! Total number of elements including defined components and lumped components + integer :: scn_nc_input + ! Initial number of single cuts considered in the oil from data input + real(pr) :: plus_w_new ! Adjusted weight fraction of the plus fraction real(pr), allocatable :: z_m_plus_i(:) - real(pr), allocatable :: plus_z_i(:) - real(pr), allocatable :: var_aux_1(:) - real(pr), allocatable :: var_aux_2(:) - integer :: i, prev_i + ! Product of molar fraction and molecular weights distribution for plus fractions + real(pr), allocatable :: plus_z_i(:) ! Molar fractions for plus fractions + real(pr), allocatable :: var_aux_1(:) ! Auxiliary variable array + real(pr), allocatable :: var_aux_2(:) ! Auxiliary variable array + integer :: i ! Iterators for loops integer, dimension(15) :: j_ps - integer :: i_ps ! auxiliary variables - real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps ! auxiliary variables - real(pr), allocatable :: plus_z_ps(:) - real(pr), allocatable :: plus_mw_ps(:) - integer :: numbers_ps - real(pr), allocatable :: density_ps(:) - real(pr), allocatable :: w_ps(:) - real(pr) :: sum_zm_last_ps - real(pr), allocatable :: moles(:) - real(pr), allocatable :: lumped_z(:) - real(pr), allocatable :: lumped_mw(:) - real(pr), allocatable :: lumped_densities(:) + ! Array for storing the number of components per pseudo-fraction + integer :: i_ps ! Index for pseudo-fractions + real(pr) :: rec_zm, remain_plus_zm, sum_z, sum_zm, sum_volume_ps + ! Auxiliary summation variables + real(pr), allocatable :: plus_z_ps(:) ! Molar fractions for pseudo-components + real(pr), allocatable :: plus_mw_ps(:) ! Molecular weight of pseudo-components + integer :: numbers_ps ! Number of pseudo-components + real(pr), allocatable :: density_ps(:) ! Densities of pseudo-components + real(pr), allocatable :: w_ps(:) ! Weight fractions of pseudo-components + real(pr) :: sum_zm_last_ps ! Summation variable for last pseudo-component + real(pr), allocatable :: moles(:) ! Molar distribution of the lumped components + real(pr), allocatable :: lumped_z(:) ! Lumped molar fractions + real(pr), allocatable :: lumped_mw(:) ! Lumped molecular weights + real(pr), allocatable :: lumped_densities(:) ! Lumped densities associate (& - def_nc => fluid%def_comp_nc , & - scn_mw => characterization%scn_mw, & - plus_zm => characterization%plus_zm, & - plus_density => characterization%plus_density, & + def_nc => fluid%def_comp_nc , & ! Number of defined components + scn_mw => characterization%scn_mw, & ! Molecular weights of SCN components + plus_zm => characterization%plus_zm, & ! Summed molar fraction of the plus fraction + plus_density => characterization%plus_density, & ! Density of plus fraction scn_zm => characterization%scn_zm, & - plus_z => characterization%plus_z, & - scn_z => characterization%scn_z, & - plus_w => fluid%plus_w, & - w => fluid%w, & - plus_mw => characterization%plus_mw, & + ! Product of molar fractions and molecular weight of SCN components + plus_z => characterization%plus_z, & ! Molar fraction of plus fraction + scn_z => characterization%scn_z, & ! Molar fraction of SCN components + plus_w => fluid%plus_w, & ! Weight fraction of plus fraction + w => fluid%w, & ! Weight fraction of all components + plus_mw => characterization%plus_mw, & ! Molecular weight of plus fraction carbon_number_plus => characterization%carbon_number_plus, & - plus6_density => characterization%plus6_density, & - C => characterization%C, & + ! Carbon numbers of plus fractions + plus6_density => characterization%plus6_density, & ! Densities of SCN components + C => characterization%C, & ! Constant for molecular weight distribution scn_nc_new => characterization%scn_nc_new & + ! Updated number of single carbon number (SCN) cuts ) - last_C = fluid%scn(fluid%scn_nc) !19 - i_last = characterization%nc_plus !124 - last = fluid%scn_nc + i_last !138 + last_C = fluid%scn(fluid%scn_nc) ! Last defined SCN component, typically 19 + i_last = characterization%nc_plus ! Number of plus components, typically 124 + last = fluid%scn_nc + i_last ! Total number of elements after lumping allocate (z_m_plus_i(0)) allocate (plus_z_i(0)) @@ -67,23 +77,23 @@ subroutine lump(fluid, characterization) allocate(lumped_mw(0)) allocate(lumped_densities(0)) - - scn_nc_input = fluid%scn_nc - scn_nc_new = fluid%scn_nc_ps - 6 - ! scn_nc_ps : CN from which all SCN fractions will be lumped + scn_nc_input = fluid%scn_nc ! Store initial SCN count + scn_nc_new = fluid%scn_nc_ps - 6 ! Define new SCN count after lumping adjustment + ! scn_nc_ps : CN from which all SCN fractions will be lumped (from input data) ! into the specified number of pseudos - characterization%plus_w = plus_w + characterization%plus_w = plus_w ! Assign weight fraction to output z_m_plus_i = characterization%product_z_mw_plus_i - plus_z_i = characterization%plus_z_i + ! Assign plus product fraction molar and molecular weight distribution + plus_z_i = characterization%plus_z_i ! Assign plus fraction molar fractions if (scn_nc_new < scn_nc_input) then - ! Plus Fraction needs to be extended to include lower CN's + ! Extend the plus fraction to include additional SCN components (lower CN's) plus_zm = plus_zm + sum(scn_zm(scn_nc_new + 1 : scn_nc_input)) plus_z = plus_z + sum(scn_z(scn_nc_new + 1 : scn_nc_input)) !extended zp plus_w_new = plus_w + sum(w(fluid%def_comp_nc+scn_nc_new + 1 & : fluid%def_comp_nc+scn_nc_input)) - !! use new variable for plus_w because type fluidata can't be modified. + ! use new variable for plus_w because type fluidata can't be modified. plus_mw = plus_zm / plus_z z_m_plus_i = scn_zm(scn_nc_new + 1 : scn_nc_input ) z_m_plus_i = [z_m_plus_i, characterization%product_z_mw_plus_i] @@ -96,17 +106,15 @@ subroutine lump(fluid, characterization) characterization%plus_z_i = plus_z_i end if - !do i= 1, scn_nc_new - ! print*, scn_z(i), scn_mw(i) - !end do + ! Lumping into pseudo-components numbers_ps = fluid%numbers_ps j_ps = 0 ! Lumping into Nps pseudos rec_zm = plus_zm / numbers_ps ! Recommended value for the product z*M (proportional to weight) ! for each pseudo according to pedersen - remain_plus_zm = plus_zm + remain_plus_zm = plus_zm ! Initialize summation variables i_ps = 1._pr sum_z = 0.0_pr sum_zm = 0.0_pr @@ -114,20 +122,23 @@ subroutine lump(fluid, characterization) i = 0.0_pr + !! Loop for pseudo-component assignment do while (i_ps < numbers_ps) i = i + 1 j_ps(i_ps) = j_ps(i_ps) + 1 sum_z = sum_z + plus_z_i(i) sum_zm = sum_zm + z_m_plus_i(i) - sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / plus6_density(scn_nc_new+i) - if (z_m_plus_i(i+1) > 2*(rec_zm - sum_zm)) then + sum_volume_ps = sum_volume_ps + z_m_plus_i(i) / & + plus6_density(scn_nc_new + i) + if (z_m_plus_i(i + 1) > 2 * (rec_zm - sum_zm)) then ! when adding one more would go too far plus_z_ps = [plus_z_ps, sum_z] plus_mw_ps = [plus_mw_ps, sum_zm/sum_z] - w_ps = [w_ps, (characterization%plus_w*(sum_zm/plus_zm))] + w_ps = [w_ps, (characterization%plus_w * (sum_zm/plus_zm))] remain_plus_zm = remain_plus_zm - sum_zm if (remain_plus_zm < rec_zm) numbers_ps = i_ps + 1 - carbon_number_plus(scn_nc_new + i_ps) = 6+((plus_mw_ps(i_ps)-84)/C) + carbon_number_plus(scn_nc_new + i_ps) = 6 + & + ((plus_mw_ps(i_ps) - 84)/C) density_ps = [density_ps, (sum_zm/sum_volume_ps)] i_ps = i_ps + 1 sum_z = 0._pr @@ -136,27 +147,25 @@ subroutine lump(fluid, characterization) end if end do - plus_z_ps = [plus_z_ps, sum(plus_z_i(i+1:i_last))] + ! Final lumping calculations + plus_z_ps = [plus_z_ps, sum(plus_z_i(i + 1 : i_last))] ! at this point, Nps is the order for the last NON-Asphaltenic pseudo comp. ! (e.g. 4 if 5 is for Asphaltenes) - plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i+1:i_last))/plus_z_ps(numbers_ps)] - w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps))* & + plus_mw_ps =[plus_mw_ps, sum(z_m_plus_i(i + 1 : i_last))/plus_z_ps(numbers_ps)] + w_ps = [w_ps, characterization%plus_w * ((plus_z_ps(numbers_ps)) * & (plus_mw_ps(numbers_ps))/(plus_zm))] carbon_number_plus(scn_nc_new + numbers_ps) = 6 + & ((plus_mw_ps(numbers_ps) - 84) / C) - sum_zm_last_ps = (sum(z_m_plus_i(i+1:i_last))) - density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / (plus6_density(scn_nc_new+i+1:last)))] - j_ps(numbers_ps) = i_last - sum(j_ps(1: numbers_ps-1)) !! revisar manana - - !do i = 1, numbers_ps - ! print*, 'ps',i, plus_z_ps(i), plus_mw_ps(i), density_ps(i), j_ps(i), carbon_number_plus(scn_nc_new+i) - !end do + sum_zm_last_ps = (sum(z_m_plus_i(i + 1 : i_last))) + density_ps = [density_ps, (sum_zm_last_ps)/sum((z_m_plus_i(i+1:i_last)) / & + (plus6_density(scn_nc_new + i + 1 : last)))] + j_ps(numbers_ps) = i_last - sum(j_ps(1 : numbers_ps - 1)) !! revisar manana - plus_density = plus_zm / sum(plus_z_ps(1:numbers_ps)* & - plus_mw_ps(1:numbers_ps)/(density_ps(1:numbers_ps))) + plus_density = plus_zm / sum(plus_z_ps(1 : numbers_ps) * & + plus_mw_ps(1 : numbers_ps)/(density_ps(1 : numbers_ps))) moles = [(fluid%def_comp_w)/(fluid%def_comp_mw)] - moles = [moles, fluid%w(def_nc+1:def_nc+scn_nc_new)/ scn_mw(1:scn_nc_new)] + moles = [moles, fluid%w(def_nc + 1:def_nc + scn_nc_new)/ scn_mw(1 : scn_nc_new)] moles = [moles, w_ps/plus_mw_ps ] characterization%mol_fraction = moles/sum(moles) ! mole fractions normalize @@ -167,11 +176,10 @@ subroutine lump(fluid, characterization) ! In this point we create a new array to content final Mis and densities ! whit pseudos. - lumped_z = [scn_z(1:scn_nc_new), plus_z_ps] - lumped_mw = [scn_mw(1:scn_nc_new), plus_mw_ps] - lumped_densities = [plus6_density(1:scn_nc_new), density_ps ] - - + lumped_z = [scn_z(1 : scn_nc_new), plus_z_ps] + lumped_mw = [scn_mw(1 : scn_nc_new), plus_mw_ps] + lumped_densities = [plus6_density(1 : scn_nc_new), density_ps ] + ! Store lumped results characterization%lumped_z = lumped_z characterization%lumped_mw = lumped_mw characterization%lumped_densities = lumped_densities diff --git a/src/molecular_weigth/mw_calculations.f90 b/src/molecular_weigth/mw_calculations.f90 index de34db4..e579e18 100644 --- a/src/molecular_weigth/mw_calculations.f90 +++ b/src/molecular_weigth/mw_calculations.f90 @@ -1,4 +1,10 @@ module molecular_weigth + !! This module contains subroutines for determining and calculating the + !! molecular weight of reservoir fluid fractions, particularly the heaviest + !! components. It allows selecting experimental data or estimated values for + !! the molecular weight and recalculates mole fractions accordingly. The + !! module is essential for compositional modeling in reservoir characterization. + use constants use dtypes, only: FluidData, FluidDataOut use bfr_routines @@ -6,23 +12,34 @@ module molecular_weigth contains subroutine select_method(fluid, mw_source, method, characterization) - !! This subroutine defines the calculation method to perform the characterization,& - !! based on the available experimental data. If the molecular weight data is experimental, - !! define the mole fractions, molecular weights and densities directly from the input file. - !! On the other hand, if the molecular weights are assumed, they are recalculated according - !! to the methodology described by Martin et al and the molar fractions of the fluid are recalculated. + !! Determines the molecular weight calculation method based on available + !! experimental data. + !! - If experimental data is available, mole fractions, molecular weights, + !! and densities are directly assigned from input. + !! - If molecular weights are assumed, they are recalculated using Martin + !! et al.'s methodology, adjusting molar fractions accordingly. implicit none - type(FluidData), intent(inout) :: fluid !! Derived type to save inlet or experimental information of fluid to be characterized. - type(FluidDataOut), intent(inout) :: characterization !! Derive type variable to save results of characterization methodology. + type(FluidData), intent(inout) :: fluid + !! Input fluid data including experimental values + type(FluidDataOut), intent(inout) :: characterization + !! Output characterization results character(len=*), intent(in), optional :: method - !! this option allows choose beetwen obtaining plus_mw or global_mw reported in the experimental information . - character(len=*), intent(in) :: mw_source !! This variable indicates the source of the input data, that is, whether they are experimental or not. - real(pr) :: sum_def_comp_z_plus !! compositions sum from (define component +1) to (plus component) - real(pr), allocatable :: def_comp_moles(:) !! - real(pr), allocatable :: scn_moles(:) !! - real(pr) :: plus_moles - real(pr) :: total_moles + !! Specifies the method for molecular weight calculation. + !! Available options: + !! - "plus_mw": reproduces the molecular weight of the residual fraction. + !! - "global_mw": reproduces the global molecular weight of fluid. + character(len=*), intent(in) :: mw_source + !! Determines whether molecular weights from predefined SCNs in the input + !! file are used or whether correlations are used to estimate these values. + !! Indicates data source ("experimental" or "calculated") + real(pr) :: sum_def_comp_z_plus + !! Sum of compositions for defined components and plus fraction + real(pr), allocatable :: def_comp_moles(:) !! Moles of defined components + real(pr), allocatable :: scn_moles(:) + !! Moles of single carbon number (SCN) fractions. + real(pr) :: plus_moles ! Moles of the plus fraction + real(pr) :: total_moles ! Total moles of the mixture associate (& scn_nc => fluid%scn_nc, def_comp_nc => fluid%def_comp_nc, & @@ -43,6 +60,8 @@ subroutine select_method(fluid, mw_source, method, characterization) select case (mw_source) case("experimental") + ! Use experimental molecular weights and calculate mole fractions + ! accordingly. scn_mw = fluid%scn_mw scn_moles = (scn_w)/(scn_mw) plus_moles = (plus_w)/(plus_mw) @@ -51,17 +70,22 @@ subroutine select_method(fluid, mw_source, method, characterization) plus_z = plus_moles / total_moles case("calculated") + ! Estimate molecular weights using correlation and recalculate + ! mole fractions if (method=="global_mw")then - scn_mw = 84 + C*(scn-6) + scn_mw = 84 + C*(scn-6) + ! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] scn_z = fluid%product_z_mw_scn / scn_mw sum_def_comp_z_plus = sum(fluid%scn_z) + (fluid%plus_z) - plus_z = sum_def_comp_z_plus - sum(scn_z) ! new molar fraction for residual cut, based on C + plus_z = sum_def_comp_z_plus - sum(scn_z) + ! Adjusted molar fraction for residual cut plus_mw = fluid%product_z_mw_plus / plus_z end if if (method=="plus_mw")then scn_mw = 84 + C*(scn-6) + ! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] scn_moles = (scn_w)/(scn_mw) plus_moles = (plus_w)/(plus_mw) total_moles = sum(def_comp_moles)+sum(scn_moles)+(plus_moles) @@ -70,6 +94,8 @@ subroutine select_method(fluid, mw_source, method, characterization) endif end select + + ! Compute adjusted mole fractions for mass consistency scn_zm = scn_z*scn_mw plus_zm = plus_z*plus_mw log_scn_z = log(scn_z) @@ -79,31 +105,39 @@ end subroutine select_method subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & characterization) - - !! this subroutine ... + !! This subroutine computes the difference between estimated and + !! experimental molecular weights for the plus fraction. implicit none - type(FluidData), intent(inout) :: fluid - type(FluidDataOut), intent(inout) :: characterization - logical :: start + type(FluidData), intent(inout) :: fluid !! Input fluid data + type(FluidDataOut), intent(inout) :: characterization + !! Output characterization data + logical :: start + !! Indicates if this is the first iteration of the calculation character(len=*), intent(in), optional :: method + !! Specifies the method for molecular weight calculation. + !! Available options: + !! - "plus_mw": reproduces the molecular weight of the residual fraction. + !! - "global_mw": reproduces the global molecular weight of fluid. character(len=*), intent(in) :: mw_source - real(pr) :: plus_mw_cal !! calculated molecular weight of residual fraction + !! Determines whether molecular weights from predefined SCNs in the input + !! file are used or whether correlations are used to estimate these values. + !! Indicates data source ("experimental" or "calculated") + real(pr) :: plus_mw_cal + !! calculated molecular weight of residual fraction real(pr), intent(out) :: difference - real(pr) :: half - real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux - real(pr):: sum_z - real(pr) :: denom - !! set to molecular weights calculated by \[M = 84-characterization%C(i-6)\] - real(pr) :: sum_def_comp_z_plus - !! compositions sum from (define component +1) to (plus component) - integer :: j, k , k_old, n_best, x_aux, i_0! internal variables - integer :: i + !! Difference between computed and expected plus molecular weight + real(pr) :: half ! Auxiliary variable + !real(pr) :: r2_old,r2_best, a_old, b_old, a_best, b_best, z_sum, z_aux + real(pr):: sum_z ! Sum moles fractions of plus fraction + real(pr) :: denom ! Denominator for C value calculation + integer :: i, i_0 ! Iteration variables integer, dimension(300) :: carbon_number_plus real(pr), dimension(300) :: plus_z_i real(pr), dimension(300) :: product_z_mw_plus_i real(pr), dimension(300) :: scn_i + ! Associate block for easier variable reference associate (& scn_nc => fluid%scn_nc, scn => fluid%scn,& plus_z => characterization%plus_z, plus_mw => characterization%plus_mw, & @@ -124,11 +158,13 @@ subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & call Best_Linear_Regression(scn_nc, scn, log_scn_z, & plus_z, a_blr, b_blr, r2, n_init, c_max_blr) call LimitLine(scn_nc, scn, plus_z, a_blr, b_blr,a_lim,b_lim,c_max_lim,half) - call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) ! add line by oscar. + call Line_C60_max (scn_nc,scn,plus_z,a_blr,b_blr,half,a_60,b_60,c_max_60) + ! add line by oscar. + ! selection best feasible regression line if(a_blr < a_lim)then - !Added by oscar 05/12/2023.se elimina por que se agregan las restricciones - !para characterization%C>12 y characterization%C<12. + ! Added by oscar 05/12/2023.se elimina por que se agregan las + ! restricciones para characterization%C>12 y characterization%C<14. a = a_lim b = b_lim else ! this line was removed call Line_C60_max @@ -141,6 +177,7 @@ subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & end if end if + ! Compute the maximun carbon number when reach molar fraction value sum_z = 0.0_pr i = 0 do while (sum_z < plus_z .and. i< 300) @@ -150,51 +187,63 @@ subroutine difference_mw_plus(fluid, mw_source, method, start, difference, & product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) carbon_number_plus(i) = i+i_0 end do - + ! Adjust C constant value if (i == 300.and.sum_z < plus_z)then if(start) C = C - 0.07 if(.not.start) C = C - 0.01 go to 1 end if + ! cmax: output CN at which characterization%plus_z is reached, as the + ! summation of single z(i) from the best linear distribution (blr) c_max = i+i_0 plus_z_i(i) = plus_z_i(i) - (sum_z-plus_z) ! Adjustment to Zp (z20+) product_z_mw_plus_i(i) = plus_z_i(i)*(84+C*(i+i_0-6)) + ! Compute directly C constant value if (mw_source == "experimental" .and. C>12 .and. C<14 ) then - ! mayor igual a 12 o menor igual 14. denom = sum((plus_z_i(1:i))*((carbon_number_plus(1:i)-6))) ! denominator of equation to obtain C value directly C = (plus_z*(plus_mw-84))/denom endif - + + ! Compute diferrence between experimental and calculated molecular weight plus_mw_cal = sum(product_z_mw_plus_i(1:i))/plus_z difference = plus_mw_cal-plus_mw scn_i = [scn, carbon_number_plus(1:i)] end associate + ! Save outputs values into characterization structure characterization%plus_z_i = plus_z_i(1:i) characterization%product_z_mw_plus_i = product_z_mw_plus_i(1:i) characterization%carbon_number_plus = carbon_number_plus(1:i) characterization%nc_plus = i characterization%scn_i = scn_i - - end subroutine difference_mw_plus subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) - !! This subroutine... + !! This subroutine determines the correlation parameter C or adjusts the + !! plus fraction molecular weight iteratively. + !! Ensures consistency between estimated and experimental values. + implicit none - type(FluidData) :: fluid - type(FluidDataOut) :: characterization - logical :: start - character(len=*), intent(in) :: mw_source + type(FluidData) :: fluid !! Input fluid data. + type(FluidDataOut) :: characterization !! Output characterization data + logical :: start + !! Indicates if this is the first iteration of the calculation + character(len=*), intent(in) :: mw_source + !! Determines whether molecular weights from predefined SCNs in the input + !! file are used or whether correlations are used to estimate these values. + !! Indicates data source ("experimental" or "calculated") character(len=*), intent(in), optional :: method - logical, intent(in) :: fix_C !! If C is < 12 or > 14 then fix to the - !! closest value. - !integer :: c_max !! output CN at which characterization%plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: log_scn_z(:) !! logarithm of corresponding mole fractions of scn cutsn + !! Specifies the method for molecular weight calculation. + !! Available options: + !! - "plus_mw": reproduces the molecular weight of the residual fraction. + !! - "global_mw": reproduces the global molecular weight of fluid. + logical, intent(in) :: fix_C + !! If present and set to .TRUE., it fixes the value of the constant 'C' + !! instead of estimating it dynamically. real(pr) :: difference !! real(pr) :: difference_old, plus_mw_old real(pr) :: C_old @@ -236,25 +285,27 @@ subroutine get_c_or_m_plus(fluid, mw_source, method, fix_C, characterization) end do end select - if(fix_C .and. C>14 .or. C<12)then + if (fix_C .and. C>14 .or. C<12) then - if(C>14) C=14 - if(C<12) C=12 + if (C > 14) C = 14 + if (C < 12) C = 12 Start = .true. plus_mw = fluid%plus_mw call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & difference_old, characterization) plus_mw_old = plus_mw - plus_mw = 0.9*plus_mw + plus_mw = 0.9 * plus_mw Start = .false. call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & difference, characterization) - - do while (abs(difference)>0.00001) + + ! Iterative adjustment to minimize difference in molecular weight calculations. + do while (abs(difference) > 0.00001) aux = plus_mw - plus_mw = plus_mw - difference*(plus_mw-plus_mw_old)/(difference-difference_old) + plus_mw = plus_mw - difference*(plus_mw-plus_mw_old) / & + (difference - difference_old) plus_mw_old = aux difference_old = difference call difference_mw_plus(fluid, mw_source, "plus_mw" , start, & diff --git a/src/types.f90 b/src/types.f90 index 242f369..ecccc7d 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -1,96 +1,140 @@ module dtypes + !! This module contains data type definitions and subroutines used for + !! manipulatingfluid characterization. The main data types are FluidData + !! and FluidDataOut, which are used to store and handle information about + !! defined components, cut fractions, and other parameters of fluid. + use constants, only: pr implicit none type :: FluidData - integer :: def_comp_nc !! number of defined components being considered in the oil - integer :: scn_nc !! number of single cuts being considered in the oil - integer :: scn_nc_ps !! CN from which all SCN fractions will be lumped into the specified number of pseudos - integer :: numbers_ps !! number of pseudos in which the scn fractions grouped. - logical :: density_setup - integer :: number_plus_density - character(len=:), allocatable :: filename - integer, allocatable :: scn (:) !! set of singles cuts being considered in the oil - character(len=15), allocatable :: def_components (:) !! set of defined components being considered in the oil + !! Structure to save data from input file + integer :: def_comp_nc + !! Number of defined components being considered in the oil + integer :: scn_nc + !! Number of single cuts being considered in the oil + integer :: scn_nc_ps + !! CN from which all SCN fractions will be lumped into the specified + !! number of pseudos + integer :: numbers_ps + !! number of pseudos in which the scn fractions grouped. + character(len=:), allocatable :: filename !! Name of input file + integer, allocatable :: scn (:) + !! set of singles cuts being considered in the oil + character(len=15), allocatable :: def_components (:) + !! set of defined components being considered in the oil character(len=15) :: scn_plus !! name of residual fraction - real(pr), allocatable :: def_comp_z(:) !! set of corresponding mole fractions of defined components - real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts - real(pr), allocatable :: def_comp_mw(:) !! set of corresponding molecular weights of defined components - real(pr), allocatable :: scn_mw(:) !! set of corresponding molecular weights of scn cuts - real(pr), allocatable:: product_z_mw_def_comp(:) !! product between composition and molecular weight of defind components - real(pr), allocatable:: product_z_mw_scn(:) !! product between composition and molecular weight of scn fractions - real(pr) :: sum_z_mw_i !! sum of the product between composition and molecular weight of the fluid's compounds + real(pr), allocatable :: def_comp_z(:) + !! set of corresponding mole fractions of defined components + real(pr), allocatable :: scn_z(:) + !! set of corresponding mole fractions of scn cuts + real(pr), allocatable :: def_comp_mw(:) + !! set of corresponding molecular weights of defined components + real(pr), allocatable :: scn_mw(:) + !! set of corresponding molecular weights of scn cuts + real(pr), allocatable:: product_z_mw_def_comp(:) + !! product between composition and molecular weight of defind components + real(pr), allocatable:: product_z_mw_scn(:) + !! product between composition and molecular weight of scn fractions + real(pr) :: sum_z_mw_i + !! sum of the product between composition and molecular weight of the + !! fluid's compounds real(pr), allocatable :: w(:) !! mass fractions of the fluid's compounds real(pr) :: plus_z !! composition of residual fraction - real(pr) :: plus6_z_exp !! composition of residual fraction - real(pr) :: plus7_z_exp !! composition of residual fraction - real(pr) :: plus12_z_exp !! composition of residual fraction - real(pr) :: plus30_z_exp !! composition of residual fraction + real(pr) :: plus6_z_exp !! composition of residual fraction, for C6+ + real(pr) :: plus7_z_exp !! composition of residual fraction, for C7+ + real(pr) :: plus12_z_exp !! composition of residual fraction, for C12+ + real(pr) :: plus30_z_exp !! composition of residual fraction, for C30+ real(pr) :: plus_mw !! molecular weight of residual fraction - real(pr) :: plus6_mw_exp !! molecular weight of residual fraction - real(pr) :: plus7_mw_exp !! molecular weight of residual fraction - real(pr) :: plus12_mw_exp !! molecular weight of residual fraction - real(pr) :: plus30_mw_exp !! molecular weight of residual fraction - real(pr) :: product_z_mw_plus !! product between composition and molecular weight of residual fraction - real(pr), allocatable :: def_comp_w(:) !! mass fractions of the defined compounds - real(pr), allocatable :: scn_w(:) !! !! mass fractions of the scn-s compounds + real(pr) :: plus6_mw_exp !! molecular weight of residual fraction, for C6+ + real(pr) :: plus7_mw_exp !! molecular weight of residual fraction, for C7+ + real(pr) :: plus12_mw_exp !! molecular weight of residual fraction, for C12+ + real(pr) :: plus30_mw_exp !! molecular weight of residual fraction, for C30+ + real(pr) :: product_z_mw_plus + !! product between composition and molecular weight of residual fraction + real(pr), allocatable :: def_comp_w(:) + !! mass fractions of the defined compounds + real(pr), allocatable :: scn_w(:) + !! mass fractions of the scn-s compounds real(pr) :: plus_w !! mass fractions of the plus fraction - real(pr), allocatable :: scn_density(:) !! set of corresponding densities of scn cuts + real(pr), allocatable :: scn_density(:) + !! set of corresponding densities of scn cuts real(pr) :: plus_density !! experimental density of the plus fraction real(pr) :: plus6_density_exp + !! experimental density of the plus fraction, for C6+ real(pr) :: plus7_density_exp + !! experimental density of the plus fraction, for C7+ real(pr) :: plus12_density_exp + !! experimental density of the plus fraction, for C12+ real(pr) :: plus30_density_exp - + !! experimental density of the plus fraction, for C30+ + end type FluidData type :: FluidDataOut - ! incluir aqui las variables que quiero que salgan como salida - type(FluidData) :: input_data - real(pr), allocatable :: scn_z(:) !! set of corresponding mole fractions of scn cuts calculated - real(pr), allocatable :: log_scn_z(:) !! set of corresponding mole fractions of scn cuts calculated - real(pr) :: plus_mw !! molecular weight of residual fraction - real(pr) :: plus_z !! composition of residual fraction - real(pr) :: C !! C constants which is used in equation \[M = 84-characterization%C(i-6)\] - real(pr), allocatable :: scn_mw(:) !! set of corresponding molecular weights of scn cuts + !! Structure to save output data + type(FluidData) :: input_data !! Structure to save data from input file + real(pr), allocatable :: scn_z(:) + !! set of corresponding mole fractions of scn cuts calculated + real(pr), allocatable :: log_scn_z(:) + !! set of corresponding logarithm of mole fractions of scn cuts calculated + real(pr) :: plus_mw !! Molecular weight of residual fraction compute + real(pr) :: plus_z !! Composition of residual fraction compute + real(pr) :: C + !! Compute C constants which is used in equation \[M = 84-characterization%C(i-6)\] + real(pr), allocatable :: scn_mw(:) + !! set of corresponding molecular weights of scn cuts compute real(pr), allocatable :: scn_zm (:) + !! Product between composition and molecular weight compute real(pr) :: plus_zm + !! Product between composition and molecular weight for residual fraction compute real(pr) :: a_blr !! A constant for best linear regression line. real(pr) :: b_blr !! B constant for best linear regression line. real(pr) :: a_60 !! A constant for Cmax60 line. real(pr) :: b_60 !! B constant for Cmax60 line. real(pr) :: a_lim !! A constant for limit feasible line. real(pr) :: b_lim !! B constant for limit feasible line. - real(pr) :: a !! output real variable. Slope of the best feasible regression line. - real(pr) :: b !! output real variable. Intercept of the best feasible regression line. - real(pr) :: r2 !! output real variable. Square correlation coefficient. - integer :: n_init ! minimum carbon number obtained from the best linear regression - integer :: c_max_blr, c_max_lim, c_max_60 + real(pr) :: a !! Slope of the best feasible regression line. + real(pr) :: b !! Intercept of the best feasible regression line. + real(pr) :: r2 !! Square correlation coefficient. + integer :: n_init !! Minimum carbon number obtained from the best linear regression + integer :: c_max_blr !! Maximum carbon number for best linear regression distribution + integer :: c_max_lim !! Maximum carbon number for limit line distribution + integer :: c_max_60 !! Maximum carbon number for C60 line distribution integer :: nc_plus !! number of elements in the distribution of C20+ - integer :: c_max !! output CN at which plus_z is reached, as the summation of single z(i) from the best linear distribution (blr) - real(pr), allocatable :: carbon_number_plus(:) - real(pr), allocatable :: plus_z_i(:) + integer :: c_max + !! output CN at which plus_z is reached, as the summation of single z(i) + !! from the best linear distribution (blr) + real(pr), allocatable :: carbon_number_plus(:) !! Carbon numbers of plus fractions + real(pr), allocatable :: plus_z_i(:) !! Composition of the residual fraction real(pr), allocatable :: product_z_mw_plus_i(:) - real(pr) :: a_d !! ad constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. - real(pr) :: b_d !! bd constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. - real(pr) :: volume_cal - real(pr) :: volume_exp - integer :: last_C - integer :: i_last - integer :: last - integer :: scn_nc_new - real(pr), allocatable :: scn_i(:) - real(pr), allocatable :: plus6_density(:) - real(pr) :: plus_w - real(pr) :: plus_density - real(pr), allocatable :: mol_fraction(:) - real(pr), allocatable :: lumped_z(:) - real(pr), allocatable :: lumped_mw(:) - real(pr), allocatable :: lumped_densities(:) + !! Product of composition and molecular weight of the residual fraction + real(pr) :: a_d + !! ad constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. + real(pr) :: b_d + !! bd constant which is used in equation \[rho_i = ad*exp(-i/10) +bd\] for density. + real(pr) :: volume_cal !! Calculated volume based on the current density function + real(pr) :: volume_exp !! Experimental volume of the fluid based on input data + integer :: last_C !! Last defined SCN component, typically 19 + integer :: i_last !! Number of plus components, typically 124 + integer :: last !! Total number of elements after lumping + integer :: scn_nc_new !! Define new SCN count after lumping adjustment + real(pr), allocatable :: scn_i(:) !! Adjusted number of single cut fractions + real(pr), allocatable :: plus6_density(:) ! Density values for the 6+ fraction + real(pr) :: plus_w !! molecular weight of residual fraction + real(pr) :: plus_density !! density of residual fraction + real(pr), allocatable :: mol_fraction(:) !! Adjusted molar fraction of the fluid + real(pr), allocatable :: lumped_z(:) !! Adjusted molar fractions after lumping + real(pr), allocatable :: lumped_mw(:) !! Adjusted molecular weights after lumping + real(pr), allocatable :: lumped_densities(:) !! Adjusted densities after lumping real(pr), allocatable :: critical_temperature(:) + !! Critical temperature compute real(pr), allocatable :: critical_pressure(:) + !! Critical pressure compute real(pr), allocatable :: acentric_factor(:) + !! Acentric factor compute real(pr), allocatable :: m_funtion (:) + !! M funtion compute to calculate acentric factor contains private @@ -101,11 +145,12 @@ module dtypes contains subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) + !! This subroutine writes the fluid characterization results to an output file. use ftools__io, only: str use defined_critical_parameters implicit none - class(FluidDataOut), intent(in) :: characterization + class(FluidDataOut), intent(in) :: characterization !! Structure to save results integer :: i, i_prev integer, intent(in) :: unit integer, intent(out) :: iostat @@ -116,6 +161,7 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) character(len=*), parameter :: num_fmt_1 = "(A3,1x,5(E15.5,2x),/)" character(len=*), parameter :: num_fmt_2 = "(A3,1x,7(E15.5,2x),/)" character(len=*), parameter :: str_fmt_2 = "(130('-'),/),/)" + character(len=*), parameter :: str_fmt_3 = "(A70,/)" associate(& def_nc => characterization%input_data%def_comp_nc, & @@ -132,32 +178,8 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) m_funtion => characterization%m_funtion & ) - - - write(unit, fmt=str_fmt_2, iostat = iostat) - !write(*,"(A64,/)") "----------------------------------------------------------------" - !write(*,*) "Best feasible regresion parameters" - !write(unit, fmt=str_fmt_1, iostat = iostat) "" - !write(unit, fmt=num_fmt_2, iostat = iostat) "Init_BFR:", characterization%n_init - !write(*,*) "A:", characterization%a," ","B:", characterization%b - !write(*,*) "C:", characterization%C - !write(*,*) "MW+:", characterization%plus_mw - !write(*,*) "Cmax:", characterization%c_max - !write(*,*) "" - !write(*,*) "Density funtion parameters parameters" - !write(*,*) "" - !write(*,*) "ad:", characterization%a_d," ","bd:", characterization%b_d - !write(*,*) "----------------------------------------------------------------" - !write(*,*) "" - !write(*,*) "compositional result" - !write(*,*) "" - - - - - - - + write(unit, fmt=str_fmt_3, iostat = iostat) "COMPOSITIONAL RESULTS" + write(unit, fmt=str_fmt_3, iostat = iostat) "" write(unit, fmt=str_fmt_1, iostat = iostat) "Comp", "Z", "Mw", "Tc", & "Pc", "Omega", "Rho", "M" ! defined components @@ -180,6 +202,30 @@ subroutine write_result(characterization,unit,iotype,v_list,iostat,iomsg) end associate end subroutine write_result + subroutine general_result(characterization) + !! This subroutine writes the fluid characterization results to an output file. + + implicit none + class(FluidDataOut), intent(in) :: characterization !! Structure to save results + + write(*,*) "-----------------------------------------------------------& + --------------------------------------------------------------" + write(*,*) "Best feasible regresion parameters" + write(*,*) "Init_BFR:", characterization%n_init + write(*,*) "A:", characterization%a," ","B:", characterization%b + write(*,*) "C:", characterization%C + write(*,*) "MW+:", characterization%plus_mw + write(*,*) "Cmax:", characterization%c_max + write(*,*) "Plus_z:", characterization%plus_z + write(*,*) "" + write(*,*) "Density funtion parameters parameters" + write(*,*) "" + write(*,*) "Ad:", characterization%a_d," ","Bd:", characterization%b_d + write(*,*) "-----------------------------------------------------------& + --------------------------------------------------------------" + write(*,*) "" + + end subroutine general_result end module dtypes