{fpe_1D_Co.PDE   FPE 1-D COllection model PDE Solutions Inc.}
title 'FPE_1D_Co' {lastest update D.Bakewell 19:30 Su/26/11/2000}
SELECT
errlim=8e-5  regrid=on  firstparts  {RGN's advice f: FIRSTPARTS}
COORDINATES
cartesian (x, y)
VARIABLES
 p     {solute probability or concentration}
DEFINITIONS
 a=0 b=20 {rectangle horiz & vert dimensns (microns)} 
 c=0.5       {rectangle 1/2 width} 
 ktz=298*1.38/(6.00*pi*9.00*1.08) { = 2.245 = KbT*(10^12)/zeta = Cdiff}
 pa0=1./(b-a)
 { ------- hyperbolic force prfile Book I p.    --- }
 Kratio=10         	{DEP/thermal energy ratio=(k/zeta)/(KbT/zeta)  }
 k_z=Kratio*ktz        {k/zeta = Force coeff/friction coeff}
 !Qy=-k_z/y  Qx=0    {solute velocity in Y & X directions}
 { ------ exponential force profile Book II p. 92 - 93; 132 - 5 ---- }
  Cv=15.38 {normalized valu - see Book II pp. }
 d=10  d2=2*d  Vo=1 k1=Cv*(Vo^2) k2=pi/d 
 Qy=-k1*exp(-k2*y) Qx=0
  pLint=BINTEGRAL(p, 'Lint')/(2*0.01)
 { ------ constant force e.g. gravity ----------- }
 !Vbgz=1.4112e-3 {gravity - bouyancy Force/friction Book II p 40-1}  
 !Qy=-Vbgz Qx=0
 { --------- steady state p(x,y) from analytical solution --------- }
 { ---  steady state: hyperbolic Book I p.    --- }
  al2=1.-Kratio {2*alpha=(k/zeta)/(KbT/zeta)}
 !ps = (al2*y^(al2-1.) )/ (b^al2-a^al2)
 !pas= (al2*a^(al2-1.) )/ (b^al2-a^al2) 
 !pTau=(pas-pa0)*(1.-exp(-1.))+pa0
 { ---  steady state: exponential Book I p.    ; Book II p.    }
 beta=k1/(k2*ktz) 
 pss=exp(beta*exp(-k2*y))
 Cexp=EXPINT(beta*exp(-k2*b)) - EXPINT(beta*exp(-k2*a))
 ps= -k2*pss/Cexp 
 pas=-k2*exp(beta*exp(-k2*a))/Cexp
 Cint=integral(pss)/(2*c)
 !ps=-k2*pss/Cint {report dif btw Cint & -Cexp/k2 } 
 { --- steady state Constant force e.g. gravity --- }
 ! Cexp=exp(-b*Vbgz/ktz)-exp(-a*Vbgz/ktz)
 ! ps=(-Vbgz/ktz)*exp(-y*Vbgz/ktz)/Cexp
 ! pas=(-Vbgz/ktz)*exp(-a*Vbgz/ktz)/Cexp
 { -- parameters common to all steady state -- } 
  pTau=(pas-pa0)*(1.-exp(-1.))+pa0
  IQx=integral(Qx) IQy=integral(Qy)
{ ------------ flux monitoring ------------- }
Fluy1= ktz*dy(p)-p*Qy   {for monitoring Flux in Y direction}
Flux1= ktz*dx(p)-p*Qx   {for monitoring Flux in X direction} 
 J= vector(p*Qx-ktz*dx(p),p*Qy-ktz*dy(p)) 
INITIAL VALUES
p=1./(b-a) 
EQUATIONS 
ktz*div[grad(p)]-dx[p*Qx]-dy[Qy*p]=dt(p)      { 2D FPE equation }
BOUNDARIES
Region 1
natural(p) = 0
start (-c,a) line to (c,a) to (c,b) to (-c,b) to finish {walk outer boundary }
Region 2
start  'Lint' (-0.01,a) line to (0.01,a)
FEATURE
start (-c,a) line to (-c+0.1,a) to (-c+0.2,a) to (-c+0.3,a) to (-c+0.4,a) to (-c+0.5,a)
start (0,a) line to (0.1,a) to (0.2,a) to (0.3,a) to (0.4,a) to (0.5,a)
start(-0.01,a+0.02) line to (0.01,a+0.02)
TIME 0 to 200000 by 0.00000001     { establish time range and initial timestep }
MONITORS
for cycle=10     contour (p) report('IQx')  report('IQy')  report('IE2')
elevation (p) from  (0,a) to (0,b)  as "X-sectional"
elevation (LOG10(p)) from(0,a) to (0,b)  as "X-sectional"
elevation (LOG10(ps)) from (0,a) to (0,b)  as "X-sectional"
contour(Qx) report 'IQx' contour(-LOG10(-Qy)) 
contour (ps) contour((p-ps)/ps) report(-Cexp/k2/Cint-1)
HISTORIES
 !history (p) at (0,a) {report(Kratio)} {report(Vbgz/ktz)} report(a) report(b) report(pas)
 ! report(pa0)  report(pTau)  
 history (p) at (0,a) report(a) report(b) report(k1) report(k2) report(pa0) report(pas)
 report(pTau)
PLOTS     { plot prints a tag-delimited data list to the file "filename.TXT":}
 grid(x,y)
 history(pLint) at (0,a) export format "#t#b#" file="c:\windows\desktop\Co_15.38_pi-d_160.txt"
! for T= 1
! elevation(p) from (0,a) to (0,b) export format "#y#b#" file="c:\windows\desktop\Co_0_grv_20.txt"
END  26300