ADVECTION AND DIFFUSION MODEL
MODELLING ADVECTION and DIFFUSION WITH FACSIMILE
DESCRIPTION
This example FACSIMILE program models a fluid flowing along a pipe 5m long at a velocity of 1 m/s. The fluid is assumed to be well mixed radially (i.e. concentrations are assumed to be homogeneous across the pipe) and 'plug flow' is assumed (velocity is constant and parallel to the pipe). The pipe is modelled by dividing it into 20 sections in this case (the number of sections can be varied depending on the accuracy and/or computation time required). Components A and C are introduced at one end of the tube at a constant rate and with an initial concentration of 10 and 1 (the units are irrelevant as long as they are consistent with the rest of the code). A and C decompose at a certain rate to produce B and D which react together at a specified rate to produce E. The effect of diffusion of the components is taken into account.
The program models the concentrations along the length of the pipe as they vary with time. It produces two output files; pipe20.out, which gives the concentrations of the components in the first section (0) and the end section (19) as they vary with time, and pipe20.snp which records the concentrations along the length of the pipe at various times.
THE DETAILS
A detailed description of the theory and practice of this type of model can be found in the FACSIMILE User Guide, Chapter 15. Chapter 10 of the same guide discusses the use of arrays in FACSIMILE and should also be reviewed when studying this code.
If you would like a copy of the facsimile input file for this application, please contact us.
1. SETTING UP THE PROBLEM
* #NCELL is number of cells ;
INTEGER #NCELL 20 #WK 200 #NTSTEP 200 ;
INTEGER #NSNAPS 12 ;
* L=Length of pipe (m), U=velocity (m/s) ;
PARAMETER H L 5 U 1 Tmp;
This section of the code defines the pipe length and velocity and sets up the number of sections to be used in the model and the required workspace (see manual for details). It also declares some parameters used in the calculation and the number of equal spaced samples of the concentrations along the pipe to be recorded.
DEFINITIONS PARAMETER A0 10 C0 1 ;
This section defines parameters for the initial concentrations of A and C and assigns their initial values.
PARAMETER K1f 1.0 K1r 0.25 K2f 0.25 K2r 0.08 K3f 100.0 ;
This section defines parameters and values for the reaction rate constants.
* ++ Diffusion coefficients ++ ;
PARAMETER DA 0.1 DB 0.1 DC 0.2 DD 0.1 DE 0.2
This section defines parameters and values for the diffusion constants for the components
* ++ Transport coefficients ++ ; PARAMETER TFA TBA TFB TBB TFC TBC TFD TBD TFE TBE
This section defines parameters for the transport coefficients
* ++ Output Parameters ++ ;
PARAMETER A_0 A_19 B_0 B_19 C_0 C_19 D_0 D_19 E_0 E_19 ;
This section defines parameters to be used for outputting the end concentration values.
PARAMETER <#WK> WORK ;
PARAMETER <#NCELL> DIST ;
PARAMETER <#NTSTEP> OUTT1 ;
PARAMETER <12> OUTT2 ;
PARAMETER TINCR ;
This section defines array parameters for use in the calculation and for storing the computed output times.
VARIABLE <#NCELL> A B C D E ;
This section defines variable arrays for use in the integration. This array stores the concentrations of each component in each section of the model.
2. THE INSTANT BLOCK
This block is executed at the start of the simulation and computes the various parameters needed.
COMPILE INSTANT;
H=L/FLOAT(#NCELL);
TFA = U/H + DA/H/H ;
TBA = DA/H/H ;
TFB = U/H + DB/H/H ;
TBB = DB/H/H ;
TFC = U/H + DC/H/H ;
TBC = DC/H/H ;
TFD = U/H + DD/H/H ;
TBD = DD/H/H ;
TFE = U/H + DE/H/H ;
TBE = DE/H/H ;
The code above computes the forward and backward transport coefficients for each of the components the velocity and the diffusion coefficients.
* Compute time step for output ;
TINCR = 5*L/U/FLOAT(#NTSTEP) ;
* fill OUTT1 array with output times ;
DO 1 FOR #1=0(1)#NTSTEP-1 ;
OUTT1<#1> = TINCR*FLOAT(#1) ;
LABEL 1 ;
* Fill distance array with midpoints of cells ;
DIST<0> = H/2.0 ;
DO 2 FOR #1=1(1)#NCELL-1 ;
DIST<#1> = DIST<0> + H*FLOAT(#1) ;
LABEL 2 ;
DO 3 FOR #1=0(1)#NSNAPS-1;
OUTT2<#1> = 2.27*FLOAT(#1) ;
LABEL 3 ;
*FOR #1=#NSNAPS-1;
*OUTT2<#1> = 2.27*FLOAT(#1) *0.995;
*WRITE 2=6,"Original value ",OUTT2<#1>;
**;
The code above computes the output times for the two output streams and stores them in the arrays OUTT1 and OUTT2 (this purely for convenience for this model, any output times can be specified and they do not need to be equispaced). It also stores the mid-point of each section in the array DIST for plotting purposes.
3. THE EQUATIONS BLOCK
The equations block is slightly more complicated for this problem because it deals with arrayed variables.
COMPILE EQUATIONS;
ARRAY <#NCELL> WORK;
% K1f % K1r : A = B;
K2f % K2r : C = D;
K3f : B + D = E;
ARRAY END;
The code above sets up the work array required for the integration process and defines the chemistry. A decomposes to B and vice versa with a forward rate of K1f and a backward rate of K1r. A similar reaction is defined between B and D. B reacts with D to give E at a rate K3f. This chemistry occurs in each section of the model.
* TRANSPORT statements ;
TRANSPORT <#NCELL> A TFA TBA ;
TRANSPORT <#NCELL> B TFB TBB ;
TRANSPORT <#NCELL> C TFC TBC ;
TRANSPORT <#NCELL> D TFD TBD ;
TRANSPORT <#NCELL> E TFE TBE ;
The code above takes care of the transport of the components between the sections of the model.
* source terms ;
%A0*U/H: = A<0> ;
C0*U/H: = C<0> ;
The code above models the introduction of A and C at the start of the pipe.
* sink terms ;
FOR #1 = #NCELL -1 ;
%U/H: A<#1> = ;
U/H: B<#1> = ;
U/H: C<#1> = ;
U/H: D<#1> = ;
U/H: E<#1> %= ;
**;
The code above models the loss of components at the end of the pipe.
4. GENERATING OUTPUT
COMPILE OUT1;
A_0 = A<0>;
B_0 = B<0>;
C_0 = C<0>;
D_0 = D<0>;
E_0 = E<0>;
FOR #1=#NCELL-1;
A_19 = A<#1>;
B_19 = B<#1>;
C_19 = C<#1>;
D_19 = D<#1>;
E_19 = E<#1>;
PSTREAM 1;
**;
The code above stores the begin and end concentrations for each component before calling the first print stream.
SETPSTREAM 1 8;
TIME A_0 A_19 B_0 B_19 C_0 C_19 D_0 D_19 E_0 E_19 ;
**;
The code above defines what will be printed in the first print stream.
COMPILE SNAPOUT;
XTIME=TIME;
WRITE 1=9, "Snapshot at time ", TIME %;
DO 1 FOR #2 = 0(1)#NCELL-1;
WRITE 1=9, ((E15,3)), DIST<#2>, A<#2>, B<#2>, C<#2>, D<#2>, E<#2>, XTIME %;
LABEL 1;
**;
The code above writes out the concentrations in each section of the pipe when called.
WHENEVER TIME = OUTT1 CALL OUT1 ;
TIME = OUTT2 % CALL SNAPOUT;
**;
The code above causes the output to be written as defined in the OUT1 and SNAPOUT blocks at times defined in the OUTT1 and OUTT2 arrays.
COMPILE INSTANT;
WRITE 1=9, "PRINT STREAM NO. 99" %;
WRITE 1=9, "Distance A B C D E XTIME" %;
**;
This code writes out the header for the pipe20.snp file at the start of the simulation.
5. SIMULATION
BEGIN;
This command starts the simulation.
STOP;
This command terminates the FACSIMILE engine when the simulation is complete.
The FACSIMILE for Windows Flow and Diffusion Wizard makes it easier to generate this type of model, and the snapshot viewer allows you to view the snapshot and output a screen shot.
