option nolet ! transientdiffusionIMPLICITSINGLETIMEsphereFINAL.TRU print "--------------------------------------------------" print " Transient Diffusion (single time=dt) (solve TDMA)" print " Transient conduction in sphere :unsteady " print "--------------------------------------------------" print " peter.vlasschaert@gmail.com " print " 5/07/2015 version 1.0 " print "-----------------------------" print " unsteady state => " print " dC/dt=Di*(d/dr(dC/dr)+2/r*dC/dr) (1) " print " C=concentration,t=time,x=space " print " IC: C=C(t,r)=C(0,r) given concentration (drug) in the gel matrix " print " BC: r=0 left -> dC/dr = 0 (2) " print " r=m*dr right -> C=C_surface=0 (3) " input prompt " example =1 , 2 do input = ":sel print if sel=1 then Di = 3e-7 ! cm^2/s dt = 100 ! s ce = 68.9 ! mg/cm^3 c_surface = 0 ! mg/cm^3 r=0.326 ! cm n=10 else input prompt " time step (sec) = ":dt input prompt " Diffusion (cm^2/s) = ":Di input prompt " initial concentratio in bead (= cylinder ) (mg/cm^3) = ":ce input prompt " surface concentration (mg/cm^3) = ":c_surface input prompt " radius of sphere ( cm ) = ":r input prompt " radius of the bead d for 'n' parts 'integer ' = ":n end if print dr = r/n n_eq=n no=n_eq mm = Di*dt/dr^2 ! big M on the website print " stability mm = ";mm dim a(1),b(1),c(1),d(1) mat redim a(no),b(no),c(no),d(no) ! combination (1) and (2) FD schemes a_1=1+1/(6*mm) a_2= -1 d1 = 1/(6*mm)*Ce d(1)=d1 ! general equations TDMA ( except last one use (3) ),2->n_eq-1 ! => inside call in_tridiag (depent on mm and m) 'see website ' ! combination (1) and (3) FD schemes q2=1-1/(n-1) a_nmin = -q2 a_n = (1/mm +2) dn=(1/mm)*ce x1=a_1 x2=a_2 x3=a_nmin x4=a_n for i=2 to n - 1 call in_tridiag(no,a(),b(),c(),x1,x2,x3,x4,i,mm) next i ! TDMA solve 'no' equations !///////////////////////////////////////////////// call in_d(no,d(),d1,dn) dim beta(1),gamma(1) mat redim beta(no),gamma(no) call beta_gamma(no,a(),b(),c(),d(),beta(),gamma()) !///////////////solve (back substitution) ////////////////// dim x(1) mat redim x(no) call back_solve (no,beta(),gamma(),c(),x()) dim aa(1,1) mat redim aa(no,no) print " d(i)=aa(i,j)*x(j) => show: d(i) " print mat print d print " matrix aa(i,j) => " print call show_tridiag(no,a(),b(),c(),aa(,)) print " column C(j) concentration distribution in the sphere => " print mat print x end sub beta_gamma(no,a(),b(),c(),d(),beta(),gamma()) beta(1) = b(1) gamma(1)= d(1)/beta(1) for j=2 to no beta (j)=b(j)-a(j)*c(j-1)/beta(j-1) gamma(j)=(d(j)-a(j)*gamma(j-1))/beta(j) next j end sub sub back_solve (no,beta(),gamma(),c(),x()) x(no)=gamma(no) ! back substitution : solve n->n-1->....1 for k=no-1 to 1 step -1 x(k)=gamma(k)-c(k)*x(k+1)/beta(k) next k end sub sub show_tridiag(no,a(),b(),c(),aa(,)) for i=1 to no !=> main diagonal aa(i,i)=b(i) !=>above diagonal if ibelow diagonal if i>=2 then aa(i,i-1)=a(i) end if next i mat print aa end sub sub in_tridiag(no,a(),b(),c(),x1,x2,x3,x4,i,mm) !=> coef main diag b(1)=x1 ! correction first line (diagonal matrix) b(i)=1/mm+2 b(no)=x4 ! correction last line (diagonal matrix) !=> coef above diag c(1)=x2 ! second correction for first line (diagonal matrix) c(i)=-(1+1/(i-1)) !=> coef below diag a(i)=-(1-1/(i-1)) a(no)=x3 ! second correction for last line (diagonal matrix) end sub sub in_d(no,d(),d1,dn) d(1)=d1 for i=2 to no-1 d(i)=dn next i d(no)=dn end sub