🔬
Solver Purpose & Physical Scope
Estimate single-stage, crosscurrent, or countercurrent liquid-liquid extraction using distribution coefficient KD and Hunter-Nash style stage balances.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 459 times
📄 Source File:
liquid_liquid_extraction.f90
📁 calcul/MassTransfer /
liquid_liquid_extraction.f90
program liquid_liquid_extraction
implicit none
integer :: mode, N, i, reqStages, iostat_val
double precision :: F, xF, S, yS, KD, targetRecovery, E, solventRatio, residualFrac
double precision :: xR, yE, recovery, totalSolute, raffSolute, extSolute, extractYavg
read(*,*,iostat=iostat_val) mode
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid mode input.'
stop
end if
read(*,*,iostat=iostat_val) F
read(*,*,iostat=iostat_val) xF
read(*,*,iostat=iostat_val) S
read(*,*,iostat=iostat_val) yS
read(*,*,iostat=iostat_val) KD
read(*,*,iostat=iostat_val) N
read(*,*,iostat=iostat_val) targetRecovery
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all extraction inputs.'
stop
end if
if (F <= 0.0d0 .or. S < 0.0d0 .or. KD <= 0.0d0) then
write(*,*) 'ERROR: Feed flow, solvent flow and KD must be valid.'
stop
end if
if (xF < 0.0d0 .or. xF >= 1.0d0 .or. yS < 0.0d0 .or. yS >= 1.0d0) then
write(*,*) 'ERROR: Compositions must be between 0 and less than 1.'
stop
end if
if (N < 1) N = 1
if (targetRecovery < 0.0d0 .or. targetRecovery >= 100.0d0) then
write(*,*) 'ERROR: Target recovery must be between 0 and less than 100 percent.'
stop
end if
totalSolute = F*xF + S*yS
solventRatio = S/F
E = KD*S/F
if (mode == 1) then
xR = (F*xF + S*yS)/(F + KD*S)
yE = KD*xR
raffSolute = F*xR
extSolute = totalSolute - raffSolute
recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
else if (mode == 2) then
residualFrac = (1.0d0/(1.0d0 + E/dble(N)))**N
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
yE = KD*xR
raffSolute = F*xR
extSolute = max(0.0d0, F*xF - raffSolute)
recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
else
residualFrac = kremser_residual(E,N)
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
yE = KD*xR
raffSolute = F*xR
extSolute = max(0.0d0, F*xF - raffSolute)
recovery = 100.0d0*extSolute/max(1.0d-30, F*xF)
end if
if (S > 0.0d0) then
extractYavg = (S*yS + extSolute)/S
else
extractYavg = 0.0d0
end if
reqStages = 1
do i=1,200
if (mode == 1) then
xR = (F*xF + S*yS)/(F + KD*S)
else if (mode == 2) then
residualFrac = (1.0d0/(1.0d0 + E/dble(i)))**i
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
else
residualFrac = kremser_residual(E,i)
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
end if
recovery = 100.0d0*(F*xF - F*xR)/max(1.0d-30,F*xF)
if (recovery >= targetRecovery) then
reqStages = i
exit
end if
reqStages = i
end do
! Recompute specified case after required-stage loop
if (mode == 1) then
xR = (F*xF + S*yS)/(F + KD*S)
yE = KD*xR
else if (mode == 2) then
residualFrac = (1.0d0/(1.0d0 + E/dble(N)))**N
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
yE = KD*xR
else
residualFrac = kremser_residual(E,N)
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
yE = KD*xR
end if
raffSolute = F*xR
extSolute = max(0.0d0, F*xF - raffSolute)
recovery = 100.0d0*extSolute/max(1.0d-30,F*xF)
if (S > 0.0d0) extractYavg = (S*yS + extSolute)/S
write(*,'(A)') '============================================================'
write(*,'(A)') ' LIQUID-LIQUID EXTRACTION ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,I8)') ' Mode = ', mode
write(*,'(A,ES12.4)') ' Feed Flow (F) = ', F
write(*,'(A,ES12.4)') ' Feed Composition x_F = ', xF
write(*,'(A,ES12.4)') ' Solvent Flow (S) = ', S
write(*,'(A,ES12.4)') ' Solvent Composition y_S = ', yS
write(*,'(A,ES12.4)') ' Distribution Coefficient (K_D) = ', KD
write(*,'(A,I8)') ' Specified Stages = ', N
write(*,*)
write(*,'(A)') '--- DESIGN RESULTS ------------------------------------------'
write(*,'(A,I8)') ' Number of Stages = ', reqStages
write(*,'(A,ES12.4)') ' Solvent-to-Feed Ratio = ', solventRatio
write(*,'(A,ES12.4)') ' Extraction Factor = ', E
write(*,'(A,ES12.4,A)') ' Solute Recovery = ', recovery, ' percent'
write(*,'(A,ES12.4)') ' Final Raffinate x_R = ', xR
write(*,'(A,ES12.4)') ' Final Extract y_E = ', extractYavg
write(*,'(A,ES12.4)') ' Equilibrium y*_E = ', yE
write(*,'(A,ES12.4)') ' Solute Extracted = ', extSolute
write(*,'(A,ES12.4)') ' Raffinate Solute = ', raffSolute
write(*,*)
write(*,'(A)') '--- STAGE PROFILE -------------------------------------------'
write(*,'(A)') ' stage x_R y_E recovery[%] raff_solute ext_solute'
write(*,'(A)') ' -----------------------------------------------------------------------------'
do i=1,max(N,reqStages)
if (mode == 1) then
xR = (F*xF + S*yS)/(F + KD*S)
else if (mode == 2) then
residualFrac = (1.0d0/(1.0d0 + E/dble(i)))**i
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
else
residualFrac = kremser_residual(E,i)
xR = xF*residualFrac + yS/KD*(1.0d0-residualFrac)
end if
yE = KD*xR
raffSolute = F*xR
extSolute = max(0.0d0,F*xF-raffSolute)
recovery = 100.0d0*extSolute/max(1.0d-30,F*xF)
if (S > 0.0d0) then
extractYavg = (S*yS + extSolute)/S
else
extractYavg = 0.0d0
end if
write(*,'(I8,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') &
i, xR, extractYavg, recovery, raffSolute, extSolute
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Distribution equilibrium: y* = K_D x.'
write(*,'(A)') ' Single stage balance: x1 = (F xF + S yS)/(F + KD S).'
write(*,'(A)') ' Crosscurrent approximation uses equal solvent split across stages.'
write(*,'(A)') ' Countercurrent Kremser/Hunter-Nash style residual fraction = (E-1)/(E^(N+1)-1).'
contains
double precision function kremser_residual(Ein,Nin)
implicit none
double precision, intent(in) :: Ein
integer, intent(in) :: Nin
if (abs(Ein-1.0d0) < 1.0d-10) then
kremser_residual = 1.0d0/dble(Nin+1)
else if (Ein > 0.0d0) then
kremser_residual = (Ein-1.0d0)/(Ein**dble(Nin+1)-1.0d0)
if (kremser_residual < 0.0d0) kremser_residual = abs(kremser_residual)
if (kremser_residual > 1.0d0) kremser_residual = 1.0d0
else
kremser_residual = 1.0d0
end if
end function kremser_residual
end program liquid_liquid_extraction
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 liquid_liquid_extraction.f90 -o liquid_liquid_extraction
2. Execution with input.txt redirection:
liquid_liquid_extraction < input.txt
📄 Sample input.txt File Structure
Sample Data:
2 100.0 0.12 80.0 0.0 2.1 4 90.0
Parameter Description:
Extraction configuration (1=Single-stage, 2=Crosscurrent, 3=Countercurrent) Feed flow F Feed solute fraction xF Solvent flow S Solvent solute fraction yS Distribution coefficient KD Number of stages N Target recovery [%]