program membrane_diffusion implicit none integer :: model, i, n_points, iostat_val double precision :: thickness, area, T_C, pFeed_bar, pPerm_bar, yAf, yAp, yBf, yBp double precision :: DA, SA, DB, SB, rpore, porosity, tau, MWA, MWB, feedFlow double precision :: T, pF, pP, R, permA, permB, permeanceA, permeanceB, JA, JB, JmassA double precision :: dPA, dPB, DKA, DKB, DKeffA, DKeffB, PsolA, PsolB, PeffA, PeffB double precision :: selectivity, flowA, flowB, stageCut, purityA double precision :: pSweep, dPAs, dPBs, JAs, JBs, flowAs, flowBs, thetaS, purityS character(len=80) :: model_name read(*,*,iostat=iostat_val) model if (iostat_val /= 0) then write(*,*) 'ERROR: Invalid model input.' stop end if read(*,*,iostat=iostat_val) thickness read(*,*,iostat=iostat_val) area read(*,*,iostat=iostat_val) T_C read(*,*,iostat=iostat_val) pFeed_bar read(*,*,iostat=iostat_val) pPerm_bar read(*,*,iostat=iostat_val) yAf read(*,*,iostat=iostat_val) yAp read(*,*,iostat=iostat_val) yBf read(*,*,iostat=iostat_val) yBp read(*,*,iostat=iostat_val) DA read(*,*,iostat=iostat_val) SA read(*,*,iostat=iostat_val) DB read(*,*,iostat=iostat_val) SB read(*,*,iostat=iostat_val) rpore read(*,*,iostat=iostat_val) porosity read(*,*,iostat=iostat_val) tau read(*,*,iostat=iostat_val) MWA read(*,*,iostat=iostat_val) MWB read(*,*,iostat=iostat_val) feedFlow if (iostat_val /= 0) then write(*,*) 'ERROR: Failed to read all membrane diffusion inputs.' stop end if if (thickness <= 0.0d0 .or. area <= 0.0d0 .or. pFeed_bar < 0.0d0 .or. pPerm_bar < 0.0d0) then write(*,*) 'ERROR: Thickness, area, and pressures must be valid positive values.' stop end if if (MWA <= 0.0d0 .or. MWB <= 0.0d0) then write(*,*) 'ERROR: Molecular weights must be positive.' stop end if if (T_C <= -273.15d0) then write(*,*) 'ERROR: Temperature must be above absolute zero.' stop end if R = 8.314462618d0 T = T_C + 273.15d0 pF = pFeed_bar * 1.0d5 pP = pPerm_bar * 1.0d5 dPA = yAf*pF - yAp*pP dPB = yBf*pF - yBp*pP PsolA = DA*SA PsolB = DB*SB DKA = (2.0d0/3.0d0)*rpore*sqrt(8.0d0*R*T/(3.141592653589793d0*(MWA/1000.0d0))) DKB = (2.0d0/3.0d0)*rpore*sqrt(8.0d0*R*T/(3.141592653589793d0*(MWB/1000.0d0))) if (tau > 0.0d0) then DKeffA = porosity/tau*DKA DKeffB = porosity/tau*DKB else DKeffA = 0.0d0 DKeffB = 0.0d0 end if select case(model) case(1) model_name = 'Solution-diffusion dense membrane' PeffA = PsolA PeffB = PsolB permeanceA = PeffA/thickness permeanceB = PeffB/thickness JA = permeanceA*dPA JB = permeanceB*dPB case(2) model_name = 'Knudsen porous membrane' PeffA = DKeffA/(R*T) PeffB = DKeffB/(R*T) permeanceA = PeffA/thickness permeanceB = PeffB/thickness JA = permeanceA*dPA JB = permeanceB*dPB case(3) model_name = 'Mixed solution-diffusion plus Knudsen resistance' PeffA = series_perm(PsolA, DKeffA/(R*T)) PeffB = series_perm(PsolB, DKeffB/(R*T)) permeanceA = PeffA/thickness permeanceB = PeffB/thickness JA = permeanceA*dPA JB = permeanceB*dPB case default write(*,*) 'ERROR: Invalid model. Use 1 solution, 2 Knudsen, 3 mixed.' stop end select if (permeanceB > 0.0d0) then selectivity = permeanceA/permeanceB else selectivity = 0.0d0 end if flowA = JA*area flowB = JB*area JmassA = JA*MWA/1000.0d0 if (feedFlow > 0.0d0) then stageCut = (flowA + flowB)/feedFlow else stageCut = 0.0d0 end if if ((flowA + flowB) > 0.0d0) then purityA = flowA/(flowA + flowB) else purityA = 0.0d0 end if write(*,'(A)') '============================================================' write(*,'(A)') ' MEMBRANE DIFFUSION ENGINE' write(*,'(A)') '============================================================' write(*,*) write(*,'(A,A)') ' Model = ', trim(model_name) write(*,'(A,ES12.4,A)') ' Thickness = ', thickness, ' m' write(*,'(A,ES12.4,A)') ' Membrane Area = ', area, ' m2' write(*,'(A,F12.4,A)') ' Temperature = ', T_C, ' deg-C' write(*,'(A,ES12.4,A)') ' Feed Pressure = ', pFeed_bar, ' bar' write(*,'(A,ES12.4,A)') ' Permeate Pressure = ', pPerm_bar, ' bar' write(*,*) write(*,'(A)') '--- TRANSPORT PROPERTIES ------------------------------------' write(*,'(A,ES12.4,A)') ' Solution D A = ', DA, ' m2/s' write(*,'(A,ES12.4,A)') ' Solution D B = ', DB, ' m2/s' write(*,'(A,ES12.4,A)') ' Solubility A = ', SA, ' mol/m3.Pa' write(*,'(A,ES12.4,A)') ' Solubility B = ', SB, ' mol/m3.Pa' write(*,'(A,ES12.4,A)') ' Knudsen D A = ', DKA, ' m2/s' write(*,'(A,ES12.4,A)') ' Knudsen D B = ', DKB, ' m2/s' write(*,'(A,ES12.4)') ' Porosity/Tortuosity = ', merge(porosity/tau, 0.0d0, tau>0.0d0) write(*,*) write(*,'(A)') '--- MEMBRANE RESULTS ----------------------------------------' write(*,'(A,ES12.4,A)') ' Effective Permeability A = ', PeffA, ' mol.m/m2.s.Pa' write(*,'(A,ES12.4,A)') ' Effective Permeability B = ', PeffB, ' mol.m/m2.s.Pa' write(*,'(A,ES12.4,A)') ' Permeance A = ', permeanceA, ' mol/m2.s.Pa' write(*,'(A,ES12.4,A)') ' Permeance B = ', permeanceB, ' mol/m2.s.Pa' write(*,'(A,ES12.4)') ' Ideal Selectivity = ', selectivity write(*,'(A,ES12.4,A)') ' Flux A = ', JA, ' mol/m2.s' write(*,'(A,ES12.4,A)') ' Flux B = ', JB, ' mol/m2.s' write(*,'(A,ES12.4,A)') ' Mass Flux A = ', JmassA, ' kg/m2.s' write(*,'(A,ES12.4,A)') ' Permeate Flow A = ', flowA, ' mol/s' write(*,'(A,ES12.4,A)') ' Permeate Flow B = ', flowB, ' mol/s' write(*,'(A,ES12.4)') ' Stage Cut Estimate = ', stageCut write(*,'(A,ES12.4)') ' Permeate Purity A = ', purityA write(*,*) write(*,'(A)') '--- PRESSURE SWEEP PROFILE ----------------------------------' write(*,'(A)') ' pFeed[bar] J_A[mol/m2.s] J_B[mol/m2.s] selectivity stage_cut purity_A' write(*,'(A)') ' ---------------------------------------------------------------------------' n_points = 45 do i=0,n_points-1 pSweep = max(pPerm_bar + 0.05d0, 0.5d0*pFeed_bar + 1.5d0*pFeed_bar*dble(i)/dble(n_points-1)) dPAs = yAf*pSweep*1.0d5 - yAp*pP dPBs = yBf*pSweep*1.0d5 - yBp*pP JAs = permeanceA*dPAs JBs = permeanceB*dPBs flowAs = JAs*area flowBs = JBs*area if (feedFlow > 0.0d0) then thetaS = (flowAs + flowBs)/feedFlow else thetaS = 0.0d0 end if if ((flowAs + flowBs) > 0.0d0) then purityS = flowAs/(flowAs + flowBs) else purityS = 0.0d0 end if write(*,'(F10.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') & pSweep, JAs, JBs, selectivity, thetaS, purityS end do write(*,*) write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------' write(*,'(A)') ' Solution-diffusion: P = D*S; permeance = P/thickness; J = permeance*Delta p.' write(*,'(A)') ' Knudsen: D_K = (2/3) r_p sqrt(8RT/(pi MW)); effective D_K = eps/tau D_K.' write(*,'(A)') ' Mixed model combines solution and Knudsen permeabilities as series resistances.' contains double precision function series_perm(P1, P2) implicit none double precision, intent(in) :: P1, P2 if (P1 <= 0.0d0) then series_perm = P2 else if (P2 <= 0.0d0) then series_perm = P1 else series_perm = 1.0d0/(1.0d0/P1 + 1.0d0/P2) end if end function series_perm end program membrane_diffusion