{Fpe2dTh.PDE - FPE PDESolns - RGN's advice f: FIRSTPARTS } title 'Fpe2DReTh: 2D FPE Relax. Thesis' {D.Bakewell 21:00 Tu/18/12/2001} SELECT errlim=8e-5 regrid=on firstparts COORDINATES cartesian (x, y) VARIABLES p {solute probability or concentration} DEFINITIONS a=0.1 b=200.1 {50.1} {rectangle horiz and vert dimnsions, respctivly} d=10 d12=d/2 d15=d/5 d10=d/10 dA=0.81 d2=2*d Vo=1.2 xc=d/2 { x center for this boundary is +ve elctrde right edge} Fcm= {0.7185} {0.5925} 0.2756 {0.0269} {0.38 0.64 0.74} {Cl-Msoti Fctr} r=.108 {r=.0465} Kdiff=0.2424/r { = 2.245 = KbT*(10^12)/zeta } Kfpe=(1.637e5)*(r^2)*(Vo^2)*Fcm/(d^3) {nrmlizd valu - see Book III, p.40} Kss=(5.373e4)*(r^3)*(Vo^2)*Fcm/(d^2) xh=(pi*x/d2)+(pi/4) yh=pi*y/d2 { h = 'hat' } ThT=(ARCTAN(sin(xh)/sinh(yh))-ARCTAN(cos(xh)/sinh(yh))) ThL=(LN(cosh(yh)+cos(xh))-LN(cosh(yh)-cos(xh))+ LN(cosh(yh)+sin(xh))-LN(cosh(yh)-sin(xh))) Dpos=cosh(2*yh)+cos(2*xh) Dneg=cosh(2*yh)-cos(2*xh) XiS=sinh(yh)*(cos(xh)/Dneg+sin(xh)/Dpos) {Xi sum} XiD=cosh(yh)*(cos(xh)/Dpos-sin(xh)/Dneg) {Xi diff} GamX=2.*ThT*XiS+ThL*XiD GamY=2.*ThT*XiD-ThL*XiS Vx=Kfpe*0 Vy=Kfpe*0 Areap=(2*dA)*1 {Areap=pi*(dA^2)/2} pA=integral(p,2)/Areap {- steady state p(x,y) from analytical solution and flux monitoring -} ThSq=4.*ThT^2+ThL^2 pss=EXP(Kss*ThSq) Cint=integral(pss) {dble intgrl of p} pInit=pss/Cint pIA=integral(pInit,2)/Areap {pInit=initial p trnsfred in} ps=1/(d)*1/(b-a) psA=integral(ps,2)/Areap pAtau=(psA-pIA)*(1.-exp(-1.))+pIA DpARel=(psA-pIA)/pIA Jx= -(Kdiff*dx(p)-p*Vx) {for monitoring Flux in X direction} Jy= -(Kdiff*dy(p)-p*Vy) {for monitoring Flux in Y direction} J= vector(p*Vx-Kdiff*dx(p),p*Vy-Kdiff*dy(p)) {Transfer in the output p file frm Fpe2dTh.dat'} transfer ("c:\windows\desktop\p0.6V200h0.108r0.72F120s.dat",pini) INITIAL VALUES p=pini EQUATIONS Kdiff*div[grad(p)]-dx[p*Vx]-dy[Vy*p]=dt(p) { 2D FPE equation } BOUNDARIES Region 1 natural(p) = 0 {walking outer boundary} start (xc-d12,a) line to (xc+d12,a) to (xc+d12,b) to (xc-d12,b) to finish {start (xc-d,a) line to (xc+d,a) to (xc+d,b) to (xc-d,b) to finish} Region 2 start 'Aint1' (xc,a) line to (xc-dA,a) to (xc-dA,a+1) to (xc+dA,a+1) to (xc+dA,a) to finish {arc (center=xc,a) to (xc+dA,a) line to finish} FEATURE start (xc-.1,a) line to (xc-.09,a) to (xc-.08,a) to (xc-.07,a) to (xc-.06,a) to (xc-.05,a) line to (xc-.04,a) to (xc-.03,a) to (xc-.02,a) to (xc-.01,a) to (xc,a) to (xc+.01,a) to (xc+.02,a) to (xc+.03,a) to (xc+.04,a) to (xc+.05,a) to (xc+.06,a) to (xc+.07,a) to (xc+.08,a) to (xc+.09,a) to (xc+.1,a) TIME 0 to 120 by 0.00000001 { time range and initial timestep (secs)} MONITORS for cycle=1 contour (p) contour (LOG10(ABS(p))) {elevation (p) from (xc,a) to (xc,b) as "X-section"} elevation (LOG10(ABS(p))) from (xc,a) to (xc,b) as "X-section" elevation (LOG10(ps)) from (xc,a) to (xc,b) as "X-section" contour(GamX) contour(GamY) ! contour(Vx) contour(Vy) {contour(LOG10(k*E2/Cdiff)) } contour (ps ) contour(LOG10(ps)) contour ((p-ps)/ps) !elevation (log10(GamX)) from (xc,a) to (xc,b) as "log10(dxE2) X-section" !elevation (Gamx) from (xc-d12,a) to (xc+d12,a) as "dxE2 X-section" !elevation (log10(-GamY)) from (xc,a) to (xc,b) as "log10(-dyE2) X-section" !elevation (GamY) from (xc-d12,a) to (xc+d12,a) as "dyE2 X-section" HISTORIES history (p) at (xc,a) history (pA) report(r) report(pIA) report(psA) report(pAtau) report(DpARel) PLOTS { plot prints a tag-delimited data list to the file "filename.TXT":} {choose for label: 0.72F 0.59F 0.276F } history(pA) export format "#t#b#" file= "c:\windows\desktop\Re1.2V200h0.108r0.72F120s.txt" elevation (LOG10(ABS(p))) from (xc,a) to (xc,b) as "X-section" export format "#x#b#y#b#" file="c:\windows\desktop\0.6V0.72Fpr_elev.txt" end