* DSMC0R.FOR * PROGRAM DSMC0R * *--test of rotational relaxation in a uniform gas. * *--SI units are used throughout * * The number of molecules=100000,no. of sub-cells,gas species=1 PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) *--All variables listed below are common variables *--MNM is the maximum number of molecules *--MNC is the maximum number of sub *--MNSC is the maximum number of sub-cells *--MNSP is the maximum number of molecular species *--MNSG is the number of species groups for collision sampling * * *--variables as defined in DSMC0.FOR * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) *--All variables listed below are common variables *--COL(M,N) is the number of collisions between species N-M molecules *--NCOL is the total number of collisions *--MOVT the total number of molecular moves *--SELT the total number of pair selections *--SEPT the sum of collision pair separations *--CS(N,M,L) sampled information on species L in cell M *----N=1 number sum *----N=2,3,4 sum of u,v,w *----N=5,6,7 sum of u*u,v*v,w*w * * *--variables as defined in DSMC0.FOR * DOUBLE PRECISION CSR(MNC,MNSP) * * CSR(M,L) the sum of the rotational energy of species L in cell M * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) *--All variables listed below are common variables *--NM is the number of molecules *--PP(M) is the x coordinate molecule M *--PV(1 to 3,M) u,v,w velocity components of molecule M *--IPL(M) sub-cell number for molecule M *--IPS(M) species code number *--IR(M) cross-reference array (molecule numbers in order of sub-cells) * COMMON /MOLSR / PR(MNM) * *--PR(M) is the rotational energy of molecule M * COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) * *--All variables listed below are common variables *--CC(M) is the cell volume *--CCG(N,M,L,K) is for collisions between species groups L-K in cell M *----N=1 is the maximum value of (relative speed)*(coll. cross-section) *----N=2 is the remainder when the selection number is rounded *--CG(N,M) is the geometry related information on cell M *----N=1 the minimum x coordinate *----N=2 the maximum x coordinate *----N=3 the cell width *--IC(N,M,L) information on the molecules of species group L in cell M *----N=1 (start address -1) of the molecule numbers in the array IR *----N=2 the number of molecules in the cell *--ISC(M) the cell in which the sub-cell lies *--ISCG(N,M,L) is the information on species group L in sub-cell M *----N=1 (start address -1) of the molecule numbers in the array IR *----N=2 the number of molecules in the sub-cell *--IG(2,M) information on group L molecules *----N=1 (start address -1) of the molecule numbers in the array IR *----N=2 the number of molecules in the cell * COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) *--All variables listed below are common variables *--SP(N,M) information on species M *----N=1 the reference cross-section (diameter in the data) *----N=2 the reference temperature *----N=3 the viscosity-temperature power law *----N=4 the reciprocal of the VSS scattering parameter *----N=5 the molecular mass *--SPM(N,M,L) information on the interaction between L-M molecules *----N=1 the reference cross-section (diameter in the data) *----N=2 the reference temperature *----N=3 the viscosity-temperature power law *----N=4 the reciprocal of the VSS scattering parameter *----N=5 the reduced mass *----N=6 the Gamma function of (5/2 - viscosity-temperature power law) *--ISP(M) the colision sampling group in which species M lies * COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) *--All variables listed below are common variables *--SPR(N,M,L) information on rotational relaxation properties of a *----species M molecule in a collision with a species L molecule *----N=1 the const. in the temperature polynomial for the collision numb *----N=2 the coefficient of T in this polynomial *----N=3 the coefficient of T**2 in this polynomial *--ISPR(N,M) integer information on rotational properties of species M *----N=1 the number of rotational degrees of freedom *----N=2 0, 1 for constant, polynomial for relaxation collision number *----N=3 0, 1 for a common, collision partner species dependent rate *--CT(M) the macroscopic temperature in cell M * COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD *--All variables listed below are common variables *--TIME time *--NPR the number of output/restart file update cycles *--NSMP the total number of samples *--FND the stream number density *--FTMP the stream temperature *--FSP(M) the fraction of species M in the stream *--ISPD relates to the setting of data for colls. between unlike mols. *----set to 0 if data is set automatically to the mean values *----set to 1 if the values are set explicitly in the data * COMMON /SAMPR / CSR * COMMON /SAMPD / CSDV(MNSP,400),CSDE(MNSP,400) *--All variables listed below are common variables *--CSDV samples the velocity distribution in 400 divisions *--CSDE samples the rotational energy distribution in 400 divisions * COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT *--All variables listed below are common variables *--FNUM is the number of real molecules represented by a simulated mol. *--DTM is the time step *--NIS is the number of time steps between samples *--NSP is the number of samples between restart and output file updates *--NPS is the estimated number of samples to steady flow *--NPT is the number of file updates to STOP * COMMON /GEOM / CW,NSC,XF,XR *--All variables listed below are common variables *--CW is the cell width *--NSC is the number of sub-cells per cell *--XF is the minimum x coordinate *--XR is the maximum x coordinate * COMMON /CONST / PI,SPI,BOLTZ *--All variables listed below are common variables *--PI is pi and SPI is the square root of pi *--BOLTZ is the Boltzmann constant * *--Basec on value of user-input to NQL call INIT0R routine for new calculation WRITE (*,*) ' INPUT 0,1 FOR CONTINUING,NEW CALCULATION:- ' READ (*,*) NQL *--Basec on value of user-input to NQLS call SAMPI0R WRITE (*,*) ' INPUT 0,1 FOR CONTINUING,NEW SAMPLE:- ' READ (*,*) NQLS * IF (NQL.EQ.1) THEN * CALL INIT0R * *--Pick input from .RES file if NQL=0-RES file is a restart file which has i/p values from a previous run ELSE * WRITE (*,*) ' READ THE RESTART FILE' OPEN (4,FILE='DSMC0R.RES',STATUS='OLD',FORM='UNFORMATTED') READ (4) BOLTZ,CC,CCG,CG,COL,CS,CSDV,CSDE,CSR,CT,CW,DTM,FNUM, & FTMP,IC,IPL,IPS,IR,ISC,ISCG,ISP,ISPR,MOVT,NCOL,NIS,NM, & NPS,NSC,NSMP,NPR,NPT,NSP,PI,PP,PR,PV,SELT,SEPT,SP,SPI, & SPM,SPR,TIME,TIMI,XF,XR CLOSE (4) * END IF * IF (NQLS.EQ.1) CALL SAMPI0R * *--NPR is incremented after each output/restart file update 100 NPR=NPR+1 * *--Call SAMPI0R untill no. of samples to steady flow is achieved IF (NPR.LE.NPS) CALL SAMPI0R * *--FOR 1 TO no. of samples between restart and output file updates *--JJ, III are local variables DO 200 JJJ=1,NSP DO 150 III=1,NIS TIME=TIME+DTM * WRITE (*,99001) III,JJJ,NIS,NSP,IFIX(NCOL) 99001 FORMAT (' DSMC0R:- Move ',2I5,' of ',2I5,I14,' Collisions') * CALL MOVE0 * CALL INDEXM * CALL COLLMR * 150 CONTINUE * CALL SAMPLE0R * 200 CONTINUE * WRITE (*,*) ' WRITING RESTART AND OUTPUT FILES',NPR,' OF ',NPT OPEN (4,FILE='DSMC0R.RES',FORM='UNFORMATTED') WRITE (4) BOLTZ,CC,CCG,CG,COL,CS,CSDV,CSDE,CSR,CT,CW,DTM,FNUM, & FTMP,IC,IPL,IPS,IR,ISC,ISCG,ISP,ISPR,MOVT,NCOL,NIS,NM, & NPS,NSC,NSMP,NPR,NPT,NSP,PI,PP,PR,PV,SELT,SEPT,SP,SPI, & SPM,SPR,TIME,TIMI,XF,XR CLOSE (4) * CALL OUT0R * IF (NPR.LT.NPT) GO TO 100 STOP END *--END of Main Body *--Start of INITOR subroutine * INIT0R.FOR * SUBROUTINE INIT0R * *--Maximum number of molecules set to 10000, Only 1 (MNSP) molecular species PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) *--All variables listed below are common variables DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) DOUBLE PRECISION CSR(MNC,MNSP) *--All variables listed below are common variables COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /MOLSR / PR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /SAMPR / CSR COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT COMMON /GEOM / CW,NSC,XF,XR COMMON /CONST / PI,SPI,BOLTZ * *--set constants- Pi, square root of Pi & Boltzman constant *--All variables listed below are common variables PI=3.141592654 SPI=SQRT(PI) BOLTZ=1.3806E-23 * CALL DATA0R **--Local Variables used in the subroutine=>M,N,L,K,VMP,MM,REM **--Here the physical constants are set. Like number **--density, temperature, molecular mass, diameter **--and all the time stepping parameters. Also seeds the particles in different cells and sub cells * *--set additional data on the gas * IF (MNSP.EQ.1) ISPD=0 *--Run the loop for each molecular species (Here MNSP=1) *--ISPD relates to the setting of data for colls. between unlike mols. Set ISPD=0 DO 100 N=1,MNSP DO 50 M=1,MNSP *--ISPR(3,N)=0 is for a common species dependent rate *--And m=n means same species IF ((ISPR(3,N).EQ.0).AND.(M.NE.N)) THEN SPR(1,N,M)=SPR(1,N,N) *--SPR(1,N,M) is const. in the temperature polynomial for the collision numb SPR(2,N,M)=SPR(2,N,N) *--SPR(2,N,M) the coefficient of T in this polynomial SPR(3,N,M)=SPR(3,N,N) *--SPR(2,N,M) the coefficient of T**2 in this polynomial END IF *--ISPD relates to the setting of data for colls. between unlike mols. Was set to ISPD=0 at begining of loop earlier IF ((ISPD.EQ.0).OR.(N.EQ.M)) THEN SPM(1,N,M)=0.25*PI*(SP(1,N)+SP(1,M))**2 *--SPM(1,N,M) is total collision C.S of the hard-sphere *--the collision cross section is assumed to be given by eqn (1.35) SPM(2,N,M)=0.5*(SP(2,N)+SP(2,M)) *--SPM(2,N,M) is ref. temp.-mean of species M & N SPM(3,N,M)=0.5*(SP(3,N)+SP(3,M)) *--SPM(3,N,M) is viscosity-temperature power law-mean of species M & N SPM(4,N,M)=0.5*(SP(4,N)+SP(4,M)) *--SPM(4,N,M) reciprocal of the VSS scattering parameter-mean of species M & N *--mean values are used for ISPD=0 ELSE SPM(1,N,M)=PI*SPM(1,N,M)**2 *--the cross-collision diameter is converted to the cross-section END IF SPM(5,N,M)=(SP(5,N)/(SP(5,N)+SP(5,M)))*SP(5,M) *--the reduced mass= { molecular mass of N * molecular mass of M } /(molecular mass of N+molecular mass of M) *--the reduced mass is defined in eqn (2.7) SPM(6,N,M)=GAM(2.5-SPM(3,N,M)) *--GAM is gamma function of (5/2-viscosity-temperature power law mean of N & M) *--SPM(6,N,M) is the Gamma function of viscosity-temperature power law 50 CONTINUE 100 CONTINUE * *--initialise variables, set NM (is the number of molecules)=0, *--set NPR(number of output/restart file update cycles)=0, set NCOL (total number of collisions)=0, *--set MOVT (total number of molecular moves)=0 *--SELT(total number of pair selections)=0,SEPT (sum of collision pair separations)=0 TIME=0. NM=0 NPR=0 NCOL=0 MOVT=0. SELT=0. SEPT=0. * *--Set all values of COL(M,N) =0 DO 200 M=1,MNSP DO 150 N=1,MNSP COL(M,N)=0. 150 CONTINUE 200 CONTINUE * *--Set minimum x-cordiate of cell 1 to XF (which is i/p parameter=minimum x coordinate) CG(1,1)=XF *--Cell width= (maximum x coordinate - minimum x coordinate)/maximum number of sub CW=(XR-XF)/MNC DO 300 M=1,MNC CT(M)=FTMP *--the macroscopic temperature is set to the free stream temperature IF (M.GT.1) CG(1,M)=CG(2,M-1) *--For M>1, set the starting x cordinate as end x cordinate of previous cell CG(2,M)=CG(1,M)+CW *--Set the cell end cordinate as begin x + cell width CG(3,M)=CW *--Setting cell width CC(M)=CW *--Setting cell volume=cell width DO 250 L=1,MNSG DO 220 K=1,MNSG *--Set remainder of CCG(1,M,L,K) (when the selection number is rounded) as a random number CCG(2,M,L,K)=RF(0) *--CCG(1,M,L,K) maximum value of (relative speed)*(coll. cross-section) for collision between L & K in cell M *--CCG(1,M,L,K)=(reference cross section of interation between L &K * 300 * (Free stream temperature/300)) CCG(1,M,L,K)=SPM(1,1,1)*300.*SQRT(FTMP/300.) 220 CONTINUE 250 CONTINUE *--the maximum value of the (rel. speed)*(cross-section) is set to a *--reasonable, but low, initial value and will be increased as necessary 300 CONTINUE * *--set sub-cells * *--For N=1 to MNC=the maximum number of sub DO 400 N=1,MNC *--For each of values of N for M=1 to NSC(number of sub-cells per cell) DO 350 M=1,NSC L=(N-1)*NSC+M ISC(L)=N *--ISC(L)=cell in which the sub-cell lies is set as N 350 CONTINUE 400 CONTINUE * *--generate initial gas with translational temperature FTMP, *--but with zero rotational temperature. * DO 500 L=1,MNSP *--For counter L = 1 to MNSP (maximum number of molecular species) *--Set reminder REM=0 REM=0 *--VMP=1/Beta VMP=SQRT(2.*BOLTZ*FTMP/SP(5,L)) *--Same as probable speed=Root of (2kT/m) m=molecular mass of species L *--VMP is the most probable speed in species L, see eqns (4.1) and (4.7) DO 450 N=1,MNC *--A = { FND() * CG(cell width) * FSP(L)-->fraction of species L in the stream / (FNUM-->no. of simulated to real molecules + REM(set to 0 earlier)) } A=FND*CG(3,N)*FSP(L)/FNUM+REM *--A is the number of simulated molecules of species L in cell N to *--simulate the required concentrations at a total number density of FND (stream number density) *--N < maximum number of sub IF (N.LT.MNC) THEN MM=A REM=(A-MM) *--Set reminder to (A-MM) *--the remainder REM is carried forward to the next cell ELSE *--Round off A to nearest whole number MM=NINT(A) END IF DO 420 M=1,MM *--For all values of no. of molecules NM < Max. no. of molecules MNM IF (NM.LE.MNM) THEN *--round-off error could have taken NM to MNM+1 *--Increment NM untill it reaches mx. no of molecules NM=NM+1 *--set Species Code No. as L (which is counter set from 1 to no of molecular species MSNP) IPS(NM)=L *--Position: x cordinate of molecules NM is set as = (min. x cordinate of N) + Random no. * (Difference of min. & max. x cordinates of N) PP(NM)=CG(1,N)+RF(0)*(CG(2,N)-CG(1,N)) *--Set sub-cell number for NM = {x cordinate of NM - min. x cordinate of N} * {No. of sub-cells per cell-0.001} / {Cell width + 1 + No. of sub-cells per cell * N-1 } IPL(NM)=(PP(NM)-CG(1,N))*(NSC-.001)/CG(3,N)+1+NSC*(N-1) *--species, position, and sub-cell number have been set *--u,v,w components which is equivalent to K=1,2 & 3 DO 405 K=1,3 *--RVELC generates random velocity components based on probable Velocity VMP & returns cosine component in A & sine components of molecule NM as u,v,w => PV(1 to 3,NM) *--based on inverse probability Maxwellian. The velocity is always three-dimensional CALL RVELC(PV(K,NM),A,VMP) 405 CONTINUE *--velocity components have been set *--set the rotational energy *--the initial rotational energy is zero for the relaxation test *--If (rotational degree of freedom for L > 0) then set rotational energy of molecule NM =0 IF (ISPR(1,L).GT.0) PR(NM)=0. END IF 420 CONTINUE 450 CONTINUE 500 CONTINUE WRITE (*,99001) NM 99001 FORMAT (' ',I6,' MOLECULES') * RETURN END * SAMPI0R.FOR * SUBROUTINE SAMPI0R * *--initialises all the sampling variables * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) DOUBLE PRECISION CSR(MNC,MNSP) * COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /SAMPR / CSR COMMON /SAMPD / CSDV(MNSP,400),CSDE(MNSP,400) COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT * NSMP=0 TIMI=TIME DO 100 L=1,MNSP DO 50 N=1,MNC CS(1,N,L)=1.E-6 DO 20 M=2,7 CS(M,N,L)=0. 20 CONTINUE CSR(N,L)=0. 50 CONTINUE IF (MNC.EQ.1) THEN DO 60 K=1,400 CSDV(L,K)=0. CSDE(L,K)=0. 60 CONTINUE END IF 100 CONTINUE RETURN END * COLLMR.FOR * *--This subroutine determines which molecules should collide. *--We assume binary collision & calculate the normalized probability *--for collision based on NTC (No Time Counter) method based on *--maximum collision cross section of the cell & time step *--The probabilty of collision is then compared to a random no. & based *--on acceptance-rejection method we determine if 2 molecules should collide SUBROUTINE COLLMR * *--calculates collisions appropriate to DTM in a gas mixture * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) DOUBLE PRECISION CSR(MNC,MNSP) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /MOLSR / PR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /SAMPR / CSR COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT COMMON /CONST / PI,SPI,BOLTZ COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N * *--VRC(3) are the pre-collision components of the relative velocity * DO 100 N=1,MNC *--For each molecule *--consider collisions in cell N *--For each of species groups for collision sampling DO 50 NN=1,MNSG DO 20 MM=1,MNSG SN=0. *--For each of molecular species DO 10 K=1,MNSP *--If colision sampling group of K *--matches with given loop counter MM *--then to SN which a variable to calculate number *--sum, add the number sumfor sampled information of *--species K in cell N IF (ISP(K).EQ.MM) SN=SN+CS(1,N,K) 10 CONTINUE *--If number sum>1 IF (SN.GT.1.) THEN *--set AVN i.e. the average number of group MM *--molecules in the cell = number sum normalized *--over the total number of samples NSMP AVN=SN/FLOAT(NSMP) ELSE *--set AVN = the number of molecules in the cell for *--species MM in cell N AVN=IC(2,N,MM) END IF *--AVN is the average number of group MM molecules in the cell *--No. of representative collisions based on NTC (No Time Counter) ASEL where CC(N) is cell volume of cell N *--i.e. Eqn. from Page 25 of classs notes ASEL=0.5*IC(2,N,NN)*AVN*FNUM*CCG(1,N,NN,MM)*DTM/CC(N) & +CCG(2,N,NN,MM) *--ASEL is the number of pairs to be selected, see eqn (11.5) NSEL=ASEL *--CCG is a variable to store the rounding value CCG(2,N,NN,MM)=ASEL-NSEL *--For more than one zero pairs (else no collisions) IF (NSEL.GT.0) THEN *--Check for no. of molecules & checking if it's not for same molecules (NN=MM) *--i.e. if insufficient molecules present-hence no collisions *--add NSEL to reminder rounding value IF (((NN.NE.MM).AND.(IC(2,N,NN).LT.1.OR.IC(2,N,MM).LT.1)) & .OR.((NN.EQ.MM).AND.(IC(2,N,NN).LT.2))) THEN CCG(2,N,NN,MM)=CCG(2,N,NN,MM)+NSEL *--if there are insufficient molecules to calculate collisions, *--the number NSEL is added to the remainer CCG(2,N,NN,MM) ELSE *--CCG=maximum value of (relative speed)*(coll. cross-section) *--Update the largest rate for comparison with CVR calculated in SELECT CVM=CCG(1,N,NN,MM) *--Add no. of pairs calculated above to total number of pair selections SELT=SELT+NSEL *--For total such calculated collisions call SELECT *--to select a potential collision pair DO 12 ISEL=1,NSEL * *--To calculate CVR based on relative vel. & collison C.S CALL SELECT * *--CVR is caluclated in SELECT subroutine *--If CVR is larger than the largest rate make the value i.e. CVM = CVR IF (CVR.GT.CVM) CVM=CVR *--Acceptance-Rejection for comparison with *--probability of collision (P25 of lecture notes) *--the probability of collision is *--compared with Rand no. for acceptance/non-acceptance IF (RF(0).LT.CVR/CCG(1,N,NN,MM)) THEN *--Acceptance loop *--the collision is accepted with the probability of eqn (11.6) *--If collision is accepted *--increment total number of collisions by 1 NCOL=NCOL+1 *--To sum of collision pair separations *--add absolute value of diff. bet. x cordinates of 2 cells L & M SEPT=SEPT+ABS(PP(L)-PP(M)) *--number of collisions between species LS-MS molecules COL(LS,MS)=COL(LS,MS)+1.D00 COL(MS,LS)=COL(MS,LS)+1.D00 * *--call INELRS which is serial *-- application of the Larsen-Borgnakke method *--Larsen-Borganakke distribution defines the division *--of energy between translational & internal modes between molecules *--For this we assume collisions as inelastic *--If number of rotational degrees of freedom of LS or MS > 0 *--i.e. bypass rotational redistribution if *--both molecules are monatomic IF (ISPR(1,LS).GT.0.OR.ISPR(1,MS).GT.0) CALL INELRS *--note the three options for the selection of rotational energy *----INELR the hierarchical application of the Larsen-Borgnakke method *----INELRS the serial application of the Larsen-Borgnakke method *----INELRA the serial application of the alternative L-B method * *--To find post-collision velocity vectors using VHS & VSS model CALL ELASTIC * END IF 12 CONTINUE *--Min. x cordinate is updated from variable CVM CCG(1,N,NN,MM)=CVM END IF END IF 20 CONTINUE 50 CONTINUE 100 CONTINUE RETURN END * INELR.FOR * SUBROUTINE INELR * *--adjustment of rotational energy in a collision *--The output is the fraction ERM which is fraction to total energy ECC *--which is redistributed for each molecule *--Local variables are L,M,LS,MS,MM,NN,N,ECI,ECF,ECC,XIB,IRT,KS,JS,XIA PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) *--All variables listed below are common COMMON /MOLSR / PR(MNM) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N * DIMENSION IR(2) *--IR is the indicator for the rotational redistribution ETI=0.5*SPM(5,LS,MS)*VRR *--ETI is the initial translational energy ECI=0. *--ECI is the initial energy in the active rotational modes ECF=0. *--ECF is the final energy in these modes ECC=ETI *--ECC is the energy to be divided *--or the available energy to be redistributed by L-R method XIB=2.5-SPM(3,LS,MS) *--XIB is th number of modes in the redistribution *--i.e. translational & rotational IRT=0 *--IRT is 0,1 if no,any redistribution is made DO 100 NSP=1,2 *--consider the molecules in turn *--NSP=no. of samples between restart and output file updates IF (NSP.EQ.1) THEN K=L KS=LS JS=MS ELSE K=M KS=MS JS=LS END IF IR(NSP)=0 *--If no. of degrees of freedom>0 i.e. molecule is monoatomic IF (ISPR(1,KS).GT.0) THEN *--If relaxation collision number is assumed to be a constant IF (ISPR(2,KS).EQ.0) THEN *--relaxation collision number= Inverse of fraction of inelastic collisions ATK=1./SPR(1,KS,JS) ELSE *--relaxation collision number with polynomial term ATK=1./(SPR(1,KS,JS)+SPR(2,KS,JS)*CT(N)+SPR(3,KS,JS)*CT(N) & **2) END IF *--ATK is the probability that rotation is redistributed to molecule L *--Applying Acceptance/Rejection method to see if ATK-probability *--of rotation redistribution can be accepted or not IF (ATK.GT.RF(0)) THEN *--Acceptance loop-Accept ATK *--IRT set to 1 if redistribution is made IRT=1 *--IR is the indicator for the rotational redistribution IR(NSP)=1 *--Add rotational energy of molecule K to energy to be divided ECC ECC=ECC+PR(K) *--Set initial energy in the active rotational modes ECI=ECI+PR(K) *--Set number of modes in the redistribution=1/2 of number of rotational degrees of freedom XIB=XIB+0.5*ISPR(1,KS) END IF END IF 100 CONTINUE *--apply the general Larsen-Borgnakke distribution function *--IRT was set to 1 if redistribution is made IF (IRT.EQ.1) THEN *--For samples between restart and output file updates 1 &2 DO 150 NSP=1,2 IF (IR(NSP).EQ.1) THEN IF (NSP.EQ.1) THEN K=L KS=LS ELSE K=M KS=MS END IF *--Update number of modes in the redistribution XIB=XIB-0.5*ISPR(1,KS) *--the current molecule is removed from the total modes *--If number of rotational degrees of freedom > 2 for given molecule IF (ISPR(1,KS).EQ.2) THEN *--ERM= Ea/(Ea+Eb) which is ratio of distributed energies for different modes ERM=1.-RF(0)**(1./XIB) ELSE *--If number of rotational degrees of freedom < 2 for given molecule XIA=0.5*ISPR(1,KS) *--This subroutine selects ratio of probability *--of a value of Ea to maximum probability CALL LBS(XIA-1.,XIB-1.,ERM) END IF *--set rotational energy of molecule based on ratio returned by LBS routine *--i.e. ERM fraction of total energy ECC PR(K)=ERM*ECC *--Reduce the redistributed rotational energy from total energy ECC=ECC-PR(K) *--the available energy is reduced accordingly *--Update the final energy for mode ECF=ECF+PR(K) END IF 150 CONTINUE *--ETF is the post-collision translational energy *--ETF= initial translational energy + initial energy in the *--active rotational modes - final energy in the mode ETF=ETI+ECI-ECF *--adjust VR and, for the VSS model, VRC for the change in energy A=SQRT(2.*ETF/SPM(5,LS,MS)) *--If reciprocal of the VSS scattering parameter < 1e-3 (VSS model) IF (ABS(SPM(4,LS,MS)-1.).LT.1.E-3) THEN *--Update relative velocity VR=A ELSE DO 160 K=1,3 *--Update the 3 (u,v,w) pre-collision components of the relative velocity VRC(K)=VRC(K)*A/VR 160 CONTINUE VR=A END IF END IF RETURN END * INELRS.FOR SUBROUTINE INELRS *--This method transfers only a fraction of the calculated change in rotational *--energy at each collision==> *--Code same except for above mentioned deviation *--adjustment of rotational energy in a collision *--alternative version with serial application of L-B method * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) *--All variables listed below are common COMMON /MOLSR / PR(MNM) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N *--Local variables are L,M,LS,MS,MM,NN,N,ECC,XIB,IRT,KS,JS,XIA,ETI ETI=0.5*SPM(5,LS,MS)*VRR *--ETI is the initial translational energy XIB=2.5-SPM(3,LS,MS) *--XIB is the number of translational modes in the redistribution *--i.e. translational & rotational IRT=0 *--IRT remains zero if there is no redistribution DO 100 NSP=1,2 *--consider the molecules in turn for serial redistribution *--NSP=no. of samples between restart and output file updates IF (NSP.EQ.1) THEN K=L KS=LS JS=MS ELSE K=M KS=MS JS=LS END IF *--If no. of degrees of freedom>0 i.e. molecule is monoatomic IF (ISPR(1,KS).GT.0) THEN *--If relaxation collision number is assumed to be a constant IF (ISPR(2,KS).EQ.0) THEN *--relaxation collision number= Inverse of fraction of inelastic collisions ATK=1./SPR(1,KS,JS) ELSE *--relaxation collision number with polynomial term ATK=1./(SPR(1,KS,JS)+SPR(2,KS,JS)*CT(N)+SPR(3,KS,JS)*CT(N) & **2) END IF *--Applying Acceptance/Rejection method to see if ATK-probability *--of rotation redistribution can be accepted or not *--ATK is the probability that rotation is redistributed to molecule L IF (ATK.GT.RF(0)) THEN *--Acceptance loop-Accept ATK *--IRT set to 1 if redistribution is made IRT=1 *--Add rotational energy of molecule K to energy to be divided ECC ECC=ETI+PR(K) *--Set initial energy in the active rotational modes ECI=PR(K) *--apply the general Larsen-Borgnakke distribution function IF (ISPR(1,KS).EQ.2) THEN ERM=1.-RF(0)**(1./XIB) ELSE *--Update number of modes in the redistribution XIA=0.5*ISPR(1,KS) CALL LBS(XIA-1.,XIB-1.,ERM) END IF PR(K)=ERM*ECC ETI=ETI+ECI-PR(K) END IF END IF 100 CONTINUE IF (IRT.EQ.1) THEN *--adjust VR and, for the VSS model, VRC for the change in energy A=SQRT(2.*ETI/SPM(5,LS,MS)) IF (ABS(SPM(4,LS,MS)-1.).LT.1.E-3) THEN VR=A ELSE DO 120 K=1,3 VRC(K)=VRC(K)*A/VR 120 CONTINUE VR=A END IF END IF RETURN END * INELRA.FOR * SUBROUTINE INELRA * *--adjustment of rotational energy in a collision *--alternative version with serial application of L-B method *--and with fractional adjustment (Larsen-Borgnakke second option) *----this option is for demonstration only and it is not recommended * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * COMMON /MOLSR / PR(MNM) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N * ETI=0.5*SPM(5,LS,MS)*VRR *--ETI is the initial translational energy XIB=2.5-SPM(3,LS,MS) *--XIB is the number of translational modes in the redistribution IRT=0 *--IRT remains zero if there is no redistribution DO 100 NSP=1,2 *--consider the molecules in turn for serial redistribution IF (NSP.EQ.1) THEN K=L KS=LS JS=MS ELSE K=M KS=MS JS=LS END IF IF (ISPR(1,KS).GT.0) THEN IF (ISPR(2,KS).EQ.0) THEN ATK=1./SPR(1,KS,JS) ELSE ATK=1./(SPR(1,KS,JS)+SPR(2,KS,JS)*CT(N)+SPR(3,KS,JS)*CT(N) & **2) END IF *--ATK is the fraction of rotation that is redistributed to molecule L IRT=1 ECC=ETI+PR(K) ECI=PR(K) *--apply the general Larsen-Borgnakke distribution function IF (ISPR(1,KS).EQ.2) THEN ERM=1.-RF(0)**(1./XIB) ELSE XIA=0.5*ISPR(1,KS) CALL LBS(XIA-1.,XIB-1.,ERM) END IF PR(K)=PR(K)+(ERM*ECC-PR(K))*ATK ETI=ETI+ECI-PR(K) END IF 100 CONTINUE IF (IRT.EQ.1) THEN *--adjust VR and, for the VSS model, VRC for the change in energy A=SQRT(2.*ETI/SPM(5,LS,MS)) IF (ABS(SPM(4,LS,MS)-1.).LT.1.E-3) THEN VR=A ELSE DO 120 K=1,3 VRC(K)=VRC(K)*A/VR 120 CONTINUE VR=A END IF END IF RETURN END * LBS.FOR * SUBROUTINE LBS(XMA,XMB,ERM) *--selects a Larsen-Borgnakke energy ratio using eqn (11.9) *--uses acceptance-rejection to check if the ratio of probability of particular value *--of Ea to max. probility is to be accepted *--Local variable used is XMA,XMB,ERM *--Set random no. to ERM 100 ERM=RF(0) *--XMA & XMB are sum of average degrees of freedom for 2 molecules IF (XMA.LT.1.E-6.OR.XMB.LT.1.E-6) THEN *--If energy levels are less than 1e-6 return 0 IF (XMA.LT.1.E-6.AND.XMB.LT.1.E-6) RETURN *--If one the energy levels are less pass the value which is > 1e-6 or significant IF (XMA.LT.1.E-6) P=(1.-ERM)**XMB IF (XMB.LT.1.E-6) P=(1.-ERM)**XMA ELSE *--Calculate the probability ratio P/Pmax P=(((XMA+XMB)*ERM/XMA)**XMA)*(((XMA+XMB)*(1.-ERM)/XMB)**XMB) END IF *--Acceptance/Rejection method to see if the ratio of probability of particular value *--of Ea to max. probility is to be accepted *--Don't accept & repeat steps if value < Rand. no. IF (P.LT.RF(0)) GO TO 100 RETURN END * SAMPLE0R.FOR * SUBROUTINE SAMPLE0R *--Local variable used is J,L,N,M,I,K,NN,LL *--Here time ensembled average in each cell of quantities such as *--the number density, velocities, rotational energies is done *--samples the molecules in the flow. * *--All variables listed below are global PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) DOUBLE PRECISION CSR(MNC,MNSP) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /MOLSR / PR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /SAMPR / CSR COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT * *--Increment total number of samples performed NSMP=NSMP+1 *--Perform sampling for each species DO 100 NN=1,MNSG *--For each cell DO 50 N=1,MNC *--Set L=No. of molecules in the cell N for species NN L=IC(2,N,NN) *--If there is at least 1 molecules in the cell IF (L.GT.0) THEN *--Run for each molecule in cell DO 10 J=1,L *--Accumlated no. of simulated molecules K=IC(1,N,NN)+J *--Pointer into cross reference M=IR(K) *--I=species code number for molecule I=IPS(M) *--Update number sum CS(1,N,I)=CS(1,N,I)+1 *--For 3 cordinates DO 5 LL=1,3 *--sum of u,v,w or accumlated V CS(LL+1,N,I)=CS(LL+1,N,I)+PV(LL,M) *--sum of u*u,v*v,w*w or accumlated V^2 CS(LL+4,N,I)=CS(LL+4,N,I)+PV(LL,M)**2 5 CONTINUE *--Accumlated rotational energy of molecules *--CSR(N,I) the sum of the rotational energy of species I in cell N CSR(N,I)=CSR(N,I)+PR(M) 10 CONTINUE END IF 50 CONTINUE 100 CONTINUE RETURN END * OUT0R.FOR * SUBROUTINE OUT0R * *--output a progressive set of results to file DSMC0R.OUT. * *--Following Output is written *-- Time flow sampled from-to *-- Total number of collisions *-- Total number of samples *-- Total number of molecules & molecular moves *-- Mean collision separation *-- Flow-field propertiesą X co-ordinate, density, translational temp, rotational temp, overall temp, stream velocity components U, V & W *-- Properties of separate species-relaxation information, rotational & translational temp *-- molecular distribution function PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) DOUBLE PRECISION CSR(MNC,MNSP) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /MOLSR / PR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /SAMPR / CSR COMMON /SAMPD / CSDV(MNSP,400),CSDE(MNSP,400) COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT COMMON /CONST / PI,SPI,BOLTZ * DOUBLE PRECISION VEL(3),SMU(3),SVEL(3,MNC),SN,SM,SMCC,SRDF,SRE,TT, & TROT,DBOLTZ DBOLTZ=BOLTZ * *--File open before writing o/p OPEN (4,FILE='DSMC0R.OUT',FORM='FORMATTED') OPEN (3,FILE='RELAX.OUT',FORM='FORMATTED',ACCESS='DIRECT',RECL=80) * *--Time flow sampled from-to WRITE (4,*) ' FLOW SAMPLED FROM TIME ',TIMI,' TO TIME ',TIME *--Total number of collisions WRITE (4,*) ' COLLISIONS:-' WRITE (4,99001) ((IFIX(COL(M,L)),M=1,MNSP),L=1,MNSP) 99001 FORMAT (5I12) *--Total no. of samples WRITE (4,*) ' TOTAL NUMBER OF SAMPLES ',NSMP WRITE (4,*) NM,' MOLECULES' WRITE (4,*) MOVT,' TOTAL MOLECULAR MOVES' *--Total molecular moves WRITE (4,*) INT(SELT),' SELECTIONS ',INT(NCOL), & ' COLLISION EVENTS, RATIO ',REAL(NCOL/SELT) IF (NCOL.GT.0) WRITE (4,*) ' MEAN COLLISION SEPARATION ', & REAL(SEPT/NCOL) WRITE (4,*) 'SAMPLES' WRITE (4,*) ' CELL N SP 1 N SP 2 ETC ' DO 100 N=1,MNC WRITE (4,99002) N,(IFIX(CS(1,N,L)),L=1,MNSP) 100 CONTINUE 99002 FORMAT (' ',I6,5I9) * WRITE (4,*) ' FLOWFIELD PROPERTIES' WRITE (4,*) &' CELL X COORD DENSITY TR TEMP ROT TEMP OV TEMP U V & W ' *--first the mixture properties-translational temp, rotational temp, overall temp, stream velocity components U, V & W DO 300 N=1,MNC A=FNUM/(CG(3,N)*NSMP) SN=0. SM=0. DO 150 K=1,3 SMU(K)=0. 150 CONTINUE SMCC=0. SRE=0. SRDF=0. DO 200 L=1,MNSP SN=SN+CS(1,N,L) *--SN is the number sum SM=SM+SP(5,L)*CS(1,N,L) *--SM is the sum of molecular masses DO 160 K=1,3 SMU(K)=SMU(K)+SP(5,L)*CS(K+1,N,L) *--SMU(1 to 3) are the sum of mu, mv, mw 160 CONTINUE SMCC=SMCC+(CS(5,N,L)+CS(6,N,L)+CS(7,N,L))*SP(5,L) *--SMCC is the sum of m(u**2+v**2+w**2) SRE=SRE+CSR(N,L) *--SRE is the sum of rotational energy SRDF=SRDF+ISPR(1,L)*CS(1,N,L) *--SRDF is the sum of the rotational degrees of freedom 200 CONTINUE DENN=SN*A *--DENN is the number density, see eqn (1.34) DEN=DENN*SM/SN *--DEN is the density, see eqn (1.42) DO 250 K=1,3 VEL(K)=SMU(K)/SM SVEL(K,N)=VEL(K) 250 CONTINUE *--VEL and SVEL are the stream velocity components, see eqn (1.43) UU=VEL(1)**2+VEL(2)**2+VEL(3)**2 TT=(SMCC-SM*UU)/(3.D00*DBOLTZ*SN) *--TT is the translational temperature, see eqn (1.51) TROT=(2.D00/DBOLTZ)*SRE/SRDF *--TROT is the rotational temperature, see eqn (11.11) TEMP=(3.D00*TT+(SRDF/SN)*TROT)/(3.+SRDF/SN) *--TEMP is the overall temperature, see eqn (11.12) CT(N)=TEMP XC=0.5*(CG(1,N)+CG(2,N)) *--XC is the x coordinate of the midpoint of the cell WRITE (4,99003) N,XC,DEN,TT,TROT,TEMP,VEL(1),VEL(2),VEL(3) 99003 FORMAT (' ',I5,F10.4,1P,E12.4,0P,6F10.4) 300 CONTINUE * WRITE (4,*) DO 500 L=1,MNSP *--now the properties of the separate species WRITE (4,*) ' SPECIES ',L WRITE (4,*) &' CELL X COORD N DENS DENSITY TR TEMP ROT TEMP U DIF & VEL V DIF VEL W DIF VEL ' DO 400 N=1,MNC A=FNUM/(CG(3,N)*NSMP) DENN=CS(1,N,L)*A *--DENN is the partial number density DEN=SP(5,L)*DENN *--DEN is the partial density, see eqn (1.13) DO 320 K=1,3 VEL(K)=CS(K+1,N,L)/CS(1,N,L) *--VEL defines the average velocity of the species L molecules 320 CONTINUE UU=VEL(1)**2+VEL(2)**2+VEL(3)**2 TT=(SP(5,L)/(3.D00*DBOLTZ)) & *((CS(5,N,L)+CS(6,N,L)+CS(7,N,L))/CS(1,N,L)-UU) *--TT is the translational temperature, see eqn (1.29) IF (ISPR(1,L).GT.0) THEN TROT=2.D00*CSR(N,L)/(ISPR(1,L)*DBOLTZ*CS(1,N,L)) ELSE TROT=0. END IF *--TROT is the rotational temperature, see eqn (11.10) DO 340 K=1,3 VEL(K)=VEL(K)-SVEL(K,N) *--VEL now defines the diffusion velocity of species L, see eqn (1.45) 340 CONTINUE XC=0.5*(CG(1,N)+CG(2,N)) WRITE (4,99004) N,XC,DENN,DEN,TT,TROT,VEL(1),VEL(2),VEL(3) 99004 FORMAT (' ',I5,F9.4,1P,2E12.4,0P,5F10.4) *--output the relaxation information IF (NPR.LE.NPS) THEN CRATE=0. DO 350 M=1,MNSP CRATE=CRATE+COL(L,M)*NSMP/CS(1,N,L) 350 CONTINUE WRITE (3,99005,REC=MNSP*(NPR-1)+L) L,CRATE,TT,TROT 99005 FORMAT (' REC ',I3,3E14.5) END IF 400 CONTINUE 500 CONTINUE * IF (MNC.EQ.1) THEN *--output the molecular distribution function WRITE (4,*) WRITE (4,*) & ' V/VM f(V/VM) Er/kT f(Er/kT)' DO 550 K=1,MNSP DKT=0.025*BOLTZ*CT(1) DVMP=0.01*SQRT(2.*BOLTZ*CT(1)/SP(5,K)) *--DKT and DVMP are the steps in kT and most probable speed DO 520 N=1,NM KK=IPS(N) IF (KK.EQ.K) THEN M=SQRT(PV(1,N)**2+PV(2,N)**2+PV(3,N)**2)/DVMP+1 IF (M.LE.400) CSDV(K,M)=CSDV(K,M)+1. M=PR(N)/DKT+1 IF (M.LE.400) CSDE(K,M)=CSDE(K,M)+1. END IF 520 CONTINUE DO 540 M=1,400 WRITE (4,99006) M*.01,CSDV(K,M)*NSP/(CS(1,1,K)*.01),M*0.025, & CSDE(K,M)*NSP/(CS(1,1,K)*0.025) 540 CONTINUE 99006 FORMAT (' ',4E13.5) 550 CONTINUE END IF * CLOSE (3) CLOSE (4) * *--check the total energy TOTE=0. DO 600 N=1,NM L=IPS(N) TOTE=TOTE+0.5*SP(5,L)*(PV(1,N)**2+PV(2,N)**2+PV(3,N)**2) TOTE=TOTE+PR(N) 600 CONTINUE WRITE (*,*) ' TOTAL ENERGY = ',TOTE * RETURN END * MOVE0.FOR * SUBROUTINE MOVE0 *--Move molecules to new position for every time step DTM=>r(t + DTM) = r(t) + {v(t) * DTM} where v(t) is velocity *--Compute new cell number IPL(N) & x cordinate position PP(N) for each of molecules NM *--the NM molecules are moved over the time interval DTM *--Also the coloumn data to be used in indexing is set up here PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT COMMON /GEOM / CW,NSC,XF,XR * *--For each of the molecules i.e. 1 to no. of molecules NM DO 100 N=1,NM *--Increment MOVT=total number of molecular moves for each molecule MOVT=MOVT+1 *--Pass sub-cell number for molecule N as MSC=new sub-cell number MSC=IPL(N) *--Set initial cell no. MC = cell in which the sub-cell lies ISC(N) *--MC is the initial cell number MC=ISC(MSC) *--XI is initial position= PP(N) which is x cordinate of molecule N XI=PP(N) *--DX is delta/distance moved during timestep DTM = velocity component u * time-step DX=PV(1,N)*DTM *--X is the new position= Initial position XI + dist. moved in DTM i.e. DX *--molecule N at XI is moved by DX to X X=XI+DX *--If new position X < XF minimum x cordinate/boundary then set X = ( 2 * (XF-X) ) i.e. for specular reflection at x=XF IF (X.LT.XF) THEN *--specular reflection from the min x boundary at x=XF (eqn (11.7)) X=2.*XF-X PV(1,N)=-PV(1,N) END IF *--If new position X > XR max x cordinate/boundary then set X = ( 2 * (XR-X) ) i.e. for specular reflection at x=XR IF (X.GT.XR) THEN *--specular reflection from the maximum x boundary at x=XR (eqn (11.7)) X=2.*XR-X PV(1,N)=-PV(1,N) END IF *--If molecular position is out of the cell cordinates viz. min-CG(1,MC) & max-CG(2,MC) IF (X.LT.CG(1,MC).OR.X.GT.CG(2,MC)) THEN *--the molecule has moved from the initial cell *--Set new cell number as diff. of x position & min. x cordinate of the cell / Cell width + 0.99999 MC=(X-XF)/CW+0.99999 *--Set MC=1 if value << 1 to avoid error in next step IF (MC.EQ.0) MC=1 *--MC is the new cell number (note avoidance of round-off error) END IF **--Set the new cell number = (current x position - min. x cordinate / Cell width) * (NSC-No. of sub-cells per cell - 0.001)+1+(NSC*new cell no.-1) MSC=((X-CG(1,MC))/CG(3,MC))*(NSC-.001)+1+NSC*(MC-1) *--MSC is the new sub-cell number *--Set new sub-cell number for molecule N as MSC IPL(N)=MSC *--Set new x-cordinate position for molecule N as X PP(N)=X 100 CONTINUE RETURN END * INDEXM.FOR * *--After having moved the molecules according to their velocities a reindexing *--is performed. i.e. in all sub-cells the number of particles are counted, *--and each molecule are given a specific address.This is basically creating the cross reference table. *--The NM molecule numbers are arranged in order of the molecule groups *--and, within the groups, in order of the cells and, within the cells, *--in order of the sub-cells SUBROUTINE INDEXM PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) * *--Foe each molecules of the species groups for collision sampling DO 200 MM=1,MNSG *-- (start address -1) of the molecule numbers in the array IR-> is reset to zero IG(2,MM)=0 *--For each cell (1 to MNC) DO 50 NN=1,MNC *--The number of molecules in the cell NN for species MM is set to zero initially IC(2,NN,MM)=0 50 CONTINUE *--For each sub-cell (1 to MNSC) DO 100 NN=1,MNSC *--number of molecules in sub-cell NN for species MM is set to zero initially ISCG(2,NN,MM)=0 100 CONTINUE 200 CONTINUE *--For each of the molecules (NM-No. of molecules) DO 300 N=1,NM *--species code number for N is passed to LS. LS is just a local counter LS=IPS(N) *--colision sampling group in which species LS lies is passed to local variable MG MG=ISP(LS) *--(start address -1) of the molecule numbers in the array IR-> is incremented by 1 IG(2,MG)=IG(2,MG)+1 *--sub-cell number for molecule N passed to MSC MSC=IPL(N) *--number of molecules in the sub-cell N is incremented by 1 ISCG(2,MSC,MG)=ISCG(2,MSC,MG)+1 *--cell in which the sub-cell lies is passed to MC MC=ISC(MSC) *--the number of molecules in the cell MC is incremented by 1 IC(2,MC,MG)=IC(2,MC,MG)+1 300 CONTINUE *--number in molecule groups in the cells and sub-cells have been counted M=0 *--For each of the number of species groups for collision sampling DO 400 L=1,MNSG *--(start address -1) of the molecule numbers in the array IR is set as M IG(1,L)=M *--the (start address -1) has been set for the groups M=M+IG(2,L) 400 CONTINUE DO 600 L=1,MNSG M=IG(1,L) DO 450 N=1,MNC *--the (start address -1) has been set as MNC IC(1,N,L)=M *--Increment M by number of molecules in the cell M=M+IC(2,N,L) 450 CONTINUE *--the (start address -1) has been set for the cells M=IG(1,L) *--For each sub-cell (1 to MNSC) DO 500 N=1,MNSC *--(start address -1) of the molecule numbers in the array IR set as M ISCG(1,N,L)=M *--Increment M by number of molecules in the sub-cell M=M+ISCG(2,N,L) *--Set number of molecules in the sub-cell=0 ISCG(2,N,L)=0 500 CONTINUE 600 CONTINUE *--the (start address -1) has been set for the sub-cells *--For each of the molecules (NM-No. of molecules) DO 700 N=1,NM *--Set species code number LS=IPS(N) *--set colision sampling group in which species LS lies MG=ISP(LS) *--new sub-cell number is set MSC=IPL(N) *--number of molecules in the sub-cell MSC is incremented by 1 ISCG(2,MSC,MG)=ISCG(2,MSC,MG)+1 *--sum up (start address -1) of the molecule numbers in the *--array & no. of molecules in sub-cell K=ISCG(1,MSC,MG)+ISCG(2,MSC,MG) *--Set cross reference array IR(K)=N *--the molecule number N has been set in the cross-reference array 700 CONTINUE RETURN END * SELECT.FOR * SUBROUTINE SELECT *--selects a potential collision pair and calculates the product of the *--collision cross-section and relative speed to be used in the COLLM subroutine =>CVR * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /CELLS / CC(MNC),CG(3,MNC),IC(2,MNC,MNSG),ISC(MNSC), & CCG(2,MNC,MNSG,MNSG),ISCG(2,MNSC,MNSG),IG(2,MNSG) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /CONST / PI,SPI,BOLTZ COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N * *--Set local K= { Rand no. * number of molecules in the cell N for species NN - 0.0001) + { (start address -1) of the molecule numbers in the array IR } + 1 K=INT(RF(0)*(IC(2,N,NN)-0.0001))+IC(1,N,NN)+1 *--Set local variable L = cross-reference array value for K or molecule numbers in order of sub-cells *--i.e. first molecule L has been chosen at random from group NN in cell L=IR(K) *--Set MSC=sub-cell number for molecule L 100 MSC=IPL(L) IF ((NN.EQ.MM.AND.ISCG(2,MSC,MM).EQ.1).OR. & (NN.NE.MM.AND.ISCG(2,MSC,MM).EQ.0)) THEN *--if MSC has no type MM molecule find the nearest sub-cell with one *--INC, NSG, NST are local variables to calculate sub-cell MSC NST=1 NSG=1 150 INC=NSG*NST NSG=-NSG NST=NST+1 MSC=MSC+INC IF (MSC.LT.1.OR.MSC.GT.MNSC) GO TO 150 IF (ISC(MSC).NE.N.OR.ISCG(2,MSC,MM).LT.1) GO TO 150 END IF *--choose the 2nd molecule M at random from the group MM *--molecules that are in the sub-cell MSC K=INT(RF(0)*(ISCG(2,MSC,MM)-0.0001))+ISCG(1,MSC,MM)+1 *--Set local variable M = cross-reference array value for K or molecule numbers in order of sub-cells *--i.e. first molecule M has been chosen at random from group MM in cell M=IR(K) *--choose a new second molecule if the first is again chosen IF (L.EQ.M) GO TO 100 DO 200 K=1,3 *--loop for 1-3 i.e. u,v & w. Set VRC are the pre-collision components of the relative velocity *--as diff. of velocity components of molecule L & M VRC(K)=PV(K,L)-PV(K,M) 200 CONTINUE *--VRC(1 to 3) are the components of the relative velocity VRR=VRC(1)**2+VRC(2)**2+VRC(3)**2 *--VR is the relative speed VR=SQRT(VRR) *--pass species code number of L & M LS=IPS(L) MS=IPS(M) *--this is the diamter at relative speed VR *--m=reduced mass (m1*m2/(m1+m2)), w=viscosity tempe. power law *-- i.e. rel. vel. * dia. * ( (2kT/m*VRR^2) }^(w-0.5) / {Gamma function (5/2 - w) } CVR=VR*SPM(1,LS,MS)*((2.*BOLTZ*SPM(2,LS,MS)/(SPM(5,LS,MS)*VRR)) & **(SPM(3,LS,MS)-0.5))/SPM(6,LS,MS) *--the collision cross-section is based on eqn (4.63) RETURN END * ELASTIC.FOR * SUBROUTINE ELASTIC * *--Local variables used are VRCP,VCCM,RML,RMM,OC,SC,A,B,C,D *--generate the post-collision velocity components. *--i.e. compute the direction of the particles after the collision, *--also known as the angle of deflection *--This subroutine has different implementation based on *--whether we use Hard Sphere Model (angles are computed *--with random picks for the hard sphere case) *--or alternatively Variable Hard Sphere (VHS) *--and Variable Soft Sphere (VSS) *--In the below implementation for DSMC0R we use VHS & VSS PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) *--Variables listed below are global COMMON /MOLS / NM,PP(MNM),PV(3,MNM),IPL(MNM),IPS(MNM),IR(MNM) COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /CONST / PI,SPI,BOLTZ COMMON /ELAST / VRC(3),VRR,VR,L,M,LS,MS,CVR,MM,NN,N * DIMENSION VRCP(3),VCCM(3) *--VRCP(3) are the post-collision components of the relative velocity *--VCCM(3) are the components of the centre of mass velocity * *--Ratio of reduced mass of MS / molecular mass of MS =>m1 *--Ratio of reduced mass of MS / molecular mass of LS =>m2 RML=SPM(5,LS,MS)/SP(5,MS) RMM=SPM(5,LS,MS)/SP(5,LS) DO 100 K=1,3 *--Compute three (K=1,2,3) components of the centre of mass velocity *--based on mass & momentum conservation i.e. m1v1=m2v2 & m1v1^2=m2v2^2 *--VCCM=m1v1/m1 + m2v2/m2 VCCM(K)=RML*PV(K,L)+RMM*PV(K,M) *--VCCM defines the components of the centre of mass velocity (eqn 2.1) 100 CONTINUE *--When scatering parameter alpha close to unity IF (ABS(SPM(4,LS,MS)-1.).LT.1.E-3) THEN *--use the VHS logic *--i.e. i.e. cos X = (2R)^(1/alpha) - 1 & alpha=1 where X= deflection angle *--alpha=scatering parameter *--B is the cosine of a random elevation angle B=2.*RF(0)-1. A=SQRT(1.-B*B) VRCP(1)=B*VR *--C is a random azimuth angle (Phi in notes) C=2.*PI*RF(0) *--cos & sine post-collision components of the relative velocity VRCP(2)=A*COS(C)*VR VRCP(3)=A*SIN(C)*VR ELSE *--use the VSS logic *--i.e. i.e. cos X = (2R)^(1/alpha) - 1 & alpha=1 where X= deflection angle *--alpha=scatering parameter *--B is the cosine of a random elevation angle *--B is the cosine of the deflection angle for the VSS model (eqn (11.8) B=2.*(RF(0)**SPM(4,LS,MS))-1. A=SQRT(1.-B*B) *--C is a random azimuth angle (Phi in notes) C=2.*PI*RF(0) *--Sine & Cosine of random azimuth angle OC=COS(C) SC=SIN(C) D=SQRT(VRC(2)**2+VRC(3)**2) *--u,v & w velocity vector equations *--Given in P38 of notes (for hard-sphere model). Here we use VHS & VSS IF (D.GT.1.E-6) THEN *--When root of sqaure of v & w are significant *--u relative = cos X * u + sin X * sin E (D)^0.5 VRCP(1)=B*VRC(1)+A*SC*D *--v relative velocity VRCP(2)=B*VRC(2)+A*(VR*VRC(3)*OC-VRC(1)*VRC(2)*SC)/D *--w relative velocity VRCP(3)=B*VRC(3)-A*(VR*VRC(2)*OC+VRC(1)*VRC(3)*SC)/D ELSE *--When root of sqaure of v & w are significant very low ignore them *--u,v & w velocity vector equations ignoring D VRCP(1)=B*VRC(1) VRCP(2)=A*OC*VRC(1) VRCP(3)=A*SC*VRC(1) END IF *--the post-collision rel. velocity components are based on eqn (2.22) END IF *--VRCP(1 to 3) are the components of the post-collision relative vel. *--by adding components of the centre of mass velocity to *--product of post-collision components of the relative velocity *--and mass of molecules *--For 3 components DO 200 K=1,3 *--Calculate velocities for L & M PV(K,L)=VCCM(K)+VRCP(K)*RMM PV(K,M)=VCCM(K)-VRCP(K)*RML 200 CONTINUE RETURN END * RVELC.FOR * SUBROUTINE RVELC(U,V,VMP) *--RVELC generates random velocity components based on variable velocity VMP *--& passes the sine value back in U & cosine value in V * *--generates two random velocity components U an V in an equilibrium *--gas with most probable speed VMP (based on eqns (C10) and (C12)) * *--A= Square root of -log(Rand)=>based on inverse probability Maxwellian A=SQRT(-LOG(RF(0))) *--B=(2*Pi*Random Number) is Random Azimuth Angle B=6.283185308*RF(0) *--U & V are Sine & Cos components of random velocity components U=A*SIN(B)*VMP V=A*COS(B)*VMP RETURN END * GAM.FOR * FUNCTION GAM(X) *--Calculates the Gamma of X which is (x-1)! when x is a positive integer * *--calculates the Gamma function of X. * A=1. Y=X IF (Y.LT.1.) THEN A=A/Y ELSE 50 Y=Y-1 IF (Y.GE.1.) THEN A=A*Y GO TO 50 END IF END IF GAM=A*(1.-0.5748646*Y+0.9512363*Y**2-0.6998588*Y**3+ & 0.4245549*Y**4-0.1010678*Y**5) RETURN END *--Random Number generation routine * * RF.FOR * FUNCTION RF(IDUM) *--generates a uniformly distributed random fraction between 0 and 1 *----IDUM will generally be 0, but negative values may be used to *------re-initialize the seed SAVE MA,INEXT,INEXTP PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1.E-9) DIMENSION MA(55) DATA IFF/0/ IF (IDUM.LT.0.OR.IFF.EQ.0) THEN IFF=1 MJ=MSEED-IABS(IDUM) MJ=MOD(MJ,MBIG) MA(55)=MJ MK=1 DO 50 I=1,54 II=MOD(21*I,55) MA(II)=MK MK=MJ-MK IF (MK.LT.MZ) MK=MK+MBIG MJ=MA(II) 50 CONTINUE DO 100 K=1,4 DO 60 I=1,55 MA(I)=MA(I)-MA(1+MOD(I+30,55)) IF (MA(I).LT.MZ) MA(I)=MA(I)+MBIG 60 CONTINUE 100 CONTINUE INEXT=0 INEXTP=31 END IF 200 INEXT=INEXT+1 IF (INEXT.EQ.56) INEXT=1 INEXTP=INEXTP+1 IF (INEXTP.EQ.56) INEXTP=1 MJ=MA(INEXT)-MA(INEXTP) IF (MJ.LT.MZ) MJ=MJ+MBIG MA(INEXT)=MJ RF=MJ*FAC IF (RF.GT.1.E-8.AND.RF.LT.0.99999999) RETURN GO TO 200 END *--Start of DATA0R subroutine * DATA0R.FOR * SUBROUTINE DATA0R * *--defines the data for a particular run of DSMC0R.FOR. * PARAMETER (MNM=100000,MNC=1,MNSC=1,MNSP=1,MNSG=1) * DOUBLE PRECISION COL(MNSP,MNSP),MOVT,NCOL,SELT,SEPT,CS(7,MNC,MNSP) * COMMON /GAS / SP(5,MNSP),SPM(6,MNSP,MNSP),ISP(MNSP) COMMON /GASR / SPR(3,MNSP,MNSP),ISPR(3,MNSP),CT(MNC) COMMON /SAMP / COL,NCOL,MOVT,SELT,SEPT,CS,TIME,NPR,NSMP,FND,FTMP, & TIMI,FSP(MNSP),ISPD COMMON /COMP / FNUM,DTM,NIS,NSP,NPS,NPT COMMON /GEOM / CW,NSC,XF,XR * *--set data (must be consistent with PARAMETER variables) * FND=1.E20 *--FND is the number densty FTMP=500. *--FTMP is the temperature FSP(1)=1. *--FSP(N) is the number fraction of species N FNUM=1.0E15 *--FNUM is the number of real molecules represented by a simulated mol. DTM=2.E-5 *--DTM is the time step NSC=1 *--NSC is the number of sub-cells in each cell XF=0. XR=1. *--the simulated region is from x=XF to x=XR--i.e. minimum & maximum X cordinate SP(1,1)=3.5E-10 SP(2,1)=273. SP(3,1)=0.75 SP(4,1)=1. SP(5,1)=5.E-26 *--SP(1,N) is the molecular diameter of species N *--SP(2,N) is the reference temperature *--SP(3,N) is the viscosity-temperatire index *--SP(4,N) is the reciprocal of the VSS scattering parameter *--SP(5,N) is the molecular mass of species N ISP(1)=1 *--ISP(N) is the group for species N ISPR(1,1)=2. SPR(1,1,1)=5. ISPR(2,1)=0 *--ISPR(1,N) is the number of degrees of freedom of species N *--SPR(1,N,M) is the constant in the polynomial for the rotational *--relaxation collision number of species N with species M *--ISPR(2,N) is 0,1 for constant, polynomial for collision number NIS=1 *--NIS is the number of time steps between samples NSP=1 *--NSP is the number of samples between restart and output file updates NPS=100 *--NPS is the number of updates to reach assumed steady flow NPT=1000 *--NPT is the number of file updates to STOP * RETURN END