SUBROUTINE INITAL C... C... Three Component Diffusion C... C... Consided the following system consisting of two bulbs connected C... by a tube C... C... C... ........ ........ C... . .+---------------ZL---------------+. . C... . .................................... . C... . VO .................................... VL . C... . . tube . . C... ........ ........ C... left bulb right bulb C... C... We make the following assumptions: C... C... (1) The pressure and temperature in the system are uniform. Also, C... the ideal gas law can be used to calculate the total molar C... density. C... C... (2) The diffusion in the tube is essentially at steady state, C... i.e., the concentration profiles in the tube adjust in- C... stantaneously to changes in the concentrations in the C... bulbs. C... C... (3) The concentration profiles in the tube are linear. C... C... (4) The fluxes of the three components in the tubes are evaluated C... at the midpoint. C... C... (5) The concentrations in the bulbs are uniform (i.e., dependent C... on time only, and not position). C... C... The following three component system is considered: C... C... i = 1; hydrogen C... C... i = 2; nitrogen C... C... i = 3; carbon dioxide C... C... The objective then is to compute the concentrations in the bulbs C... as a function of time. C... C... The following equations are used to model the system (capital C... letters are used for consistency with the Fortran programming) C... C... Component balances on the left bulb: C... C... dXIO/dt = -NI*A/(CT*VO), I = 1,2 (1) C... C... where C... C... XIO mole fraction of component i (= I) in the left bulb (gm C... mols I/gm mol gas) C... C... t time (s) C... C... NI flux of component I in the tube (gm mols of I/s-cm**2) C... C... A tube cross sectional area (cm**2) C... C... CT total molar density (gm mols gas/cm**3) C... C... VO volume of the left bulb (cm**3) C... C... Similarly a component balance on the right bulb gives (with "O" in C... eq. (1) replaced by "L") C... C... dXIL/dt = NL*A/(CT*VL), I = 1,2 (2) C... C... The third component concentration and flux can be computed from C... C... X1O + X2O + X3O = 1 (3) C... C... X1L + X2L + X3L = 1 (4) C... C... N1 + N2 + N3 = 0 (5) C... C... The fluxes are related to the gradients by the Stefan-Maxwell C... equation [Bird, et al (1960), pp 569-572] C... C... NC XI*NJ - XJ*NI C... DEL(XI) = SUM ------------- (6) C... J=1 CT*DIJ C... I NE J C... where C... C... DEL "del" differential operator C... C... NC number of components (= 3) C... C... DIJ binary diffusion coefficient for components I and J C... C... For I = 1, eq. (6) is C... C... 3 XI*NJ - XJ*NI C... DEL(X1) = SUM ------------- C... J=2 CT*DIJ C... J NE 1 C... C... = (1/(CT*D12))*(X1*N2 - X2*N1) + (1/(CT*D13))*(X1*N3 - X3*N1) (7) C... C... Eq. (7) can be written in standard algebraic form as C... C... A11*N1 + A12*N2 = B1 (8) C... C... where C... C... A11 = (1/(CT*D12))*(-X2) + (1/(CT*D13))*(-X1-X3) C... C... A12 = (1/(CT*D12))*( X1) + (1/(CT*D13))*(-X1) C... C... B1 = DEL(X1) C... C... For I = 2, eq. (6) is C... C... 3 XI*NJ - XJ*NI C... DEL(X2) = SUM ------------- C... J=1 CT*DIJ C... J NE 2 C... C... = (1/(CT*D21))*(X2*N1 - X1*N2) + (1/(CT*D23))*(X2*N3 - X3*N2) (9) C... C... Eq. (9) can be written in standard algebraic form as C... C... A21*N1 + A22*N2 = B2 (10) C... C... where C... C... A21 = (1/(CT*D21))*( X2) + (1/(CT*D23))*(-X2) C... C... A22 = (1/(CT*D21))*(-X1) + (1/(CT*D23))*(-X2-X3) C... C... B2 = DEL(X2) C... C... For the H2-N2-CO2 SYSTEM, the following diffusivities (cm**2/s) C... are reported by Taylor, et al (1993), pp 105-110, 131-133 C... C... D12 = 8.33D-01 C... C... D13 = 6.8D-01 C... C... D23 = 1.68D-01 C... C... Eqs. (1) to (11) are programmed in the following coding to compute C... the three component concentrations and fluxes. C... C... Bird, R. B., (1960), Transport Phenomena, Wiley, New York C... C... Taylor, R., and R. Krishna, (1993), Multicomponent Mass Transfer, C... Wiley, New York C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NOSTP, NORUN + /Y/ X1O, X2O, X1L, X2L + /F/ X1OT, X2OT, X1LT, X2LT + /C/ A, ZL, VO, VL, + TC, P, R, CT, + D12, D13, D23, X1, X2, + X3O, X3L, DELX1, DELX2, + N1, N2, N3 + /I/ IP DOUBLE PRECISION N1, N2, N3 C... C... MODEL PARAMETERS C... C... TUBE LENGTH (CM) ZL=10.0D0 C... C... TUBE CROSS SECTIONAL AREA (CM**2) A=0.04D0 C... C... TERMINAL VOLUMES (CM**3) VO=80.0D0 VL=80.0D0 C... C... TEMPERATURE (C) TC=35.2D0 C... C... PRESSURE (PA) P=101300.0D0 C... C... GAS CONSTANT (PA-CM**3/K-GM MOL) R=101300.0D0*22400.0D0/273.16D0 C... C... MOLAR DENSITY (GM MOL/CM**3) CT=P/(R*(TC+273.16D0)) C... C... BINARY DIFFUSIVITIES (CM**2/HR) D12=8.33D-01*3600.0D0 D13= 6.8D-01*3600.0D0 D23=1.68D-01*3600.0D0 C... C... INITIAL CONDITIONS X1O=0.0D0 X2O=0.5D0 X1L=0.5D0 X2L=0.5D0 C... C... INITIAL DERIVATIVES CALL DERV IP=0 RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/T/ T, NOSTP, NORUN + /Y/ X1O, X2O, X1L, X2L + /F/ X1OT, X2OT, X1LT, X2LT + /C/ A, ZL, VO, VL, + TC, P, R, CT, + D12, D13, D23, X1, X2, + X3O, X3L, DELX1, DELX2, + N1, N2, N3 + /I/ IP DOUBLE PRECISION N1, N2, N3 C... C... GRADIENTS DELX1=(X1L-X1O)/ZL DELX2=(X2L-X2O)/ZL C... C... MIDPOINT CONCENTRATIONS X1=(X1L+X1O)/2.0D0 X2=(X2L+X2O)/2.0D0 C... C... DEFINE LINEAR ALGEBRAIC EQUATIONS FOR N1 AND N2 X3=1.0D0-X1-X2 A11=1.0D0/(CT*D12)*(-X2)+(1.0D0/(CT*D13))*(-X1-X3) A12=1.0D0/(CT*D12)*( X1)+(1.0D0/(CT*D13))*(-X1) B1=DELX1 A21=1.0D0/(CT*D12)*( X2)+(1.0D0/(CT*D23))*(-X2) A22=1.0D0/(CT*D12)*(-X1)+(1.0D0/(CT*D23))*(-X2-X3) B2=DELX2 C... C... SOLUTION OF LINEAR ALGEBRAIC EQUATIONS BY KRAMER'S RULE C... C... DISCRIMINANT DIS=A11*A22-A12*A21 C... C... CHECK FOR NEAR-SINGULAR SYSTEM IF(DABS(DIS).LT.1.0D-20)THEN WRITE(*,1)DIS,T 1 FORMAT(' Discriminant = ',D12.3,' at t = ',F6.2) STOP END IF C... C... FLUXES N1=(A22*B1-A12*B2)/DIS N2=(A11*B2-A21*B1)/DIS C... C... ODES X1OT=-N1*A/(CT*VO) X2OT=-N2*A/(CT*VO) X1LT= N1*A/(CT*VL) X2LT= N2*A/(CT*VL) RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NT=101) COMMON/T/ T, NOSTP, NORUN + /Y/ X1O, X2O, X1L, X2L + /F/ X1OT, X2OT, X1LT, X2LT + /C/ A, ZL, VO, VL, + TC, P, R, CT, + D12, D13, D23, X1, X2, + X3O, X3L, DELX1, DELX2, + N1, N2, N3 + /I/ IP DOUBLE PRECISION N1, N2, N3 C... C... MONITOR THE OUTPUT IP=IP+1 WRITE(*,*)IP,T C... C... HEADING FOR THE SOLUTION IF(IP.EQ.1)WRITE(NO,1) 1 FORMAT(9X,'t',7X,'X1O',7X,'X2O',7X,'X3O',/, + 17X,'X1L',7X,'X2L',7X,'X3L',/, + 17X,' N1',7X,' N2',7X,' N3',/, + 3X,'dX1O/dt',3X,'dX2O/dt',3X,'dX1L/dt',3X,'dX2L/dt') C... C... WRITE THE SOLUTION EVERY 10TH CALL TO PRINT X3O=1.0D0-X1O-X2O X3L=1.0D0-X1L-X2L IF(((IP-1)/10*10).EQ.(IP-1))THEN WRITE(NO,2)T,X1O,X2O,X3O,X1L,X2L,X3L, + N1,N2,N3,X1OT,X2OT,X1LT,X2LT 2 FORMAT(F10.2,3F10.5,/, + 10X,3F10.5,/, + 10X,3D10.2,/, + 4D10.2,//) END IF C... C... PLOT THE SOLUTION WITH MATLAB CALL MPLOT RETURN END SUBROUTINE MPLOT IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER(NT=101) COMMON/T/ T, NOSTP, NORUN + /Y/ X1O, X2O, X1L, X2L + /F/ X1OT, X2OT, X1LT, X2LT + /C/ A, ZL, VO, VL, + TC, P, R, CT, + D12, D13, D23, X1, X2, + X3O, X3L, DELX1, DELX2, + N1, N2, N3 + /I/ IP DOUBLE PRECISION N1, N2, N3 C... C... OPEN FILES FOR MATLAB PLOTTING IF(IP.EQ.1)THEN OPEN(1,FILE='ma3c.out') OPEN(2,FILE='ma3f.out') END IF C... C... WRITE THE SOLUTION FOR PLOTTING C... C... CONCENTRATIONS X3O=1.0D0-X1O-X2O X3L=1.0D0-X1L-X2L WRITE(1,1)T, X1O, X2O, X3O, X1L, X2L, X3L 1 FORMAT(7F9.5) C... C... FLUXES N3=-N1-N2 WRITE(2,2)T, N1, N2, N3 2 FORMAT(F8.4,3F13.7) RETURN END