option nolet print "------------------------------------------" print " galerkin 'finite element method' " print " weight function '= shape function' " print " 1D dT/dx(k*dT/dx)+S = 0 ,S=0 version 1.0 " print " peter.vlasschaert@gmail.com , 28/07/2015 " print " BC : x=0,T=T0 ; x=L,T=TL " print " find : qx=-A*k*dT/dx at x=0,x=L ;A=1 " print "------------------------------------------" n_div = 5 ! 'integer' number finite elements T0 = 30 ! '°C ' temperature at x=0 TL = 100 ! '°C ' temperarature at x=L ! thermal conductivity kk=1.2 ! 'W/(cm*°C)' print "thermal conductivity ( W/cm°c) = ";kk ! info : plane slab L = 10 ! 'cm' dx = 10/n_div ! 'cm' ! stifness coef:for each each element (e)='NiNj' AA= 1 ! 1D coef =1 ! build : stiffness matrix [k(e)],linear dim ke(2,2) ke(1,1)=1 ke(1,2)=-1 ke(2,1)=-1 ke(2,2)=1 ! build element ' matrix ' :sum [k(e)] nt = n_div * 2 dim kh(2,2) mat redim kh(nt,2) for i=1 to nt step 2 kh(i,1)=coef*ke(1,1) kh(i,2)=coef*ke(1,2) kh(i+1,1)=coef*ke(2,1) kh(i+1,2)=coef*ke(2,2) next i mat print kh ! global stiffness matrix :found by summing these equations [ks] no=n_div+1 dim ks(1,1) mat redim ks(no,no) !mat print ks for k=1 to no-1 for i=k to k+1 for j=k to k+1 ks(i,j)=kh(i+k-1,j-(k-1)) next j next i next k ! diagonal : global matrix correction ks(2,2) = kh(2,2)+kh(3,1) ! ************************************** !ks(3,3) = kh(4,2)+kh(5,1) ! 3 -> 4 -> 5 !ks(4,4) = kh(6,2)+kh(7,1) ! 4 -> 6 -> 7 !ks(5,5) = kh(8,2)+kh(9,1) ! 5 -> 8 -> 8 !*************************************** ! see : sub below to do this !*************************************** k=0 for i=3 to no-1 k=k+1 ks(i,i)=kh(i+k,2)+kh(i+k+1,1) next i print " global stiffness matix = " mat print ks ! Next step is to incorporate BC in the stiffness matrix ! 1e) Step delete ( first and last Row from [ks] matrix ) ! T1=T0 T2,T3,T4,T5 T6=TL ! inside temperature distribution : T2,T3,T4,T5 dim kf(1,1) mat redim kf(no,no) mat kf = 0 ! T1,T6 known T1 removes first element from the matrix = 0 ! T6 removes last element from the matrix = 0 ! -1 sign -> other side -> d(1)= T1=T0 (TDMA) ! -1 sign -> other side -> d(4)= T6=TL (TDMA) for i=2 to 5 ! start 2 to 4,T2,T3,T4,T5 for j=1 to 6 kf(i-1,j)=ks(i,j) next j next i ! removes zero rows mat redim kf(no-2,no) ! delete elements kf(1,1) = 0 kf(no-2,no) = 0 print " TDMA system " v=no-2 ! number of equations for TDMA dim a(1),b(1),c(1),d(1) mat redim a(v),b(v),c(v),d(v) mat print kf a(1)=0 b(1)=2 c(1)=-1 d(1)= T0 for i=2 to v-1 a(i) = -1 b(i) = 2 c(i) =-1 d(i) = 0 next i a(v) = -1 b(v) = 2 c(v) = 0 d(v) = TL dim beta(1),gamma(1) mat redim beta(v),gamma(v) call beta_gamma(v,a(),b(),c(),d(),beta(),gamma()) dim x(1) mat redim x(v) call back_solve (v,beta(),gamma(),c(),x()) print "temperature distribution inside the plane slab ( °C) = " mat print x print " total temperature distribution ( °C) = " print dim g(1) mat redim g(no) g(1) = T0 for i=2 to v+1 g(i) = x( i-1 ) next i g(no) = TL for i=1 to no print g(i);" °C", next i print print print " heat transfer rate ( q = -k*dT/dx ) = " print mq = -(g(2)-g(1))/dx print " dT/dx | (x=0) =" ;mq ; " °C/cm" mp = (g(6)-g(5))/dx print " dT/dx | (x=L) =" ;mp ; " °C/cm" print print " qx ( x=0) = "; -kk*mq; " W/cm^2" print " qx ( x=L) = "; -kk*mp; " W/cm^2" 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