🔬
Solver Purpose & Physical Scope
Estimate gas or vapor transport through dense or porous membranes using solution-diffusion, Knudsen diffusion, or mixed resistance models. Outputs permeance, selectivity, and transmembrane flux.
📂 Discipline: Masstransfer
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 255 times
📄 Source File:
membrane_diffusion.f90
📁 calcul/MassTransfer /
membrane_diffusion.f90
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 membrane_diffusion.f90 -o membrane_diffusion
2. Execution with input.txt redirection:
membrane_diffusion < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1.0e-6 1.0 25.0 5.0 1.0 0.15 0.80 0.85 0.20 1.0e-10 3.0e-2 2.0e-11 5.0e-3 50e-9 0.35 3.0 44.01 28.014 1.0
Parameter Description:
Membrane Model (1=Dense Solution-Diffusion, 2=Porous Knudsen, 3=Mixed Resistance) Membrane thickness delta [m] Membrane area A [m2] Temperature [C] Feed pressure [bar] Permeate pressure [bar] Feed mole fraction A Permeate mole fraction A Feed mole fraction B Permeate mole fraction B Diffusivity A [m2/s] Solubility A [mol/m3-Pa] Diffusivity B [m2/s] Solubility B [mol/m3-Pa] Pore radius rp [m] Porosity epsilon Tortuosity tau Molar mass A [g/mol] Molar mass B [g/mol] Feed molar flow rate [mol/s]