c case2_chamber.for c calculates stress and deflection in a chamber consisting of two c flat endplates in the form of circular rings, supported apart by c two cylinders. c based on case2_step.for -- which was for a single endplate c Roark & Young, table 24, case 2, with a step implicit real*8 (a-h), (o-z) c endplate index: 1 = front, 2 = rear common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta dimension amatf(4,4), cvecf(4), vbegf(4), oswire(3,99), area(99), . twire(3), Ewire(3), dwire(3), EAwire(3), fwire(3), . loadfac(3,99), ncelltot(99), nwire(3), dhole(3), . dzf(2,99), amatr(4,4), cvecr(4), vbegr(4) c bi = inner radius of region i data bi / 0.235, ! inner radius of plate . 0.2530, 0.2662, 0.2781, 0.2899, 0.3016, ! superlayer 1 . 0.3124, 0.3243, 0.3361, 0.3479, 0.3596, ! superlayer 2 . 0.3646, 0.3764, 0.3882, 0.4000, 0.4118, ! superlayer 3 . 0.4168, 0.4286, 0.4404, 0.4522, 0.4639, ! superlayer 4 . 0.4747, 0.4867, 0.4984, 0.5101, 0.5219, ! superlayer 5 . 0.5269, 0.5387, 0.5505, 0.5623, 0.5741, ! superlayer 6 . 0.5791, 0.5909, 0.6027, 0.6145, 0.6262, ! superlayer 7 . 0.6370, 0.6489, 0.6606, 0.6724, 0.6842, ! superlayer 8 . 0.6892, 0.7010, 0.7128, 0.7246, 0.7364, ! superlayer 9 . 0.7414, 0.7532, 0.7650, 0.7768, 0.7900, ! superlayer 10 . 0.8090, 47*0.0 / ! outer radius of plate c ihole = 1 -- region has holes data ihole / 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, . 0, 4*1, 0, 4*1, 0, 4*1, 0, 4*1, 49*0 / c layer = drift chamber layer number data layer / 0, 1, 2, 3, 4, 0, 5, 6, 7, 8, . 0, 9, 10, 11, 12, 0, 13, 14, 15, 16, . 0, 17, 18, 19, 20, 0, 21, 22, 23, 24, . 0, 25, 26, 27, 28, 0, 29, 30, 31, 32, . 0, 33, 34, 35, 36, 0, 37, 38, 39, 40, . 49*0 / c desired wire tensions in grams (sense, clear, field) data twire / 34., 86., 182. / c data twire / 34., 86., 86. / ! if use 80 micron wires for field wires c modulus of wire (Pa) data Ewire / 3.55e11, 2*7.0e10 / c diameter of wires (m) data dwire / 0.00002, 0.00008, 0.00012 / c diameter of holes (m) data dhole / 0.0045, 2*0.0025 / c cells per layer data ncelltot / 0, 4*96, 0, 4*112, 0, 4*128, 0, 4*144, 0, 4*176, . 0, 4*192, 0, 4*208, 0, 4*224, 0, 4*240, 0, 4*256, . 49*0 / c number of wires per cell (sense, clear, field) data loadfac / . 3*0, 1,2,3, 1,0,2, 1,0,2, 1,2,2, ! 1 extra field wire . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, ! in first and last layer . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,2, . 3*0, 1,2,2, 1,0,2, 1,0,2, 1,2,3, . 147*0/ pi = 4.*atan(1.) c set up wire parameters do i = 1,3 fwire(i) = twire(i)*9.8/1000. ! tension in N EAwire(i) = Ewire(i)*pi*dwire(i)**2/4. nwire(i) = 0 enddo ! i zlen = 2.764 ! wire length (m) c setup for plate divided into 51 regions n = 51 rstep = 0.47 ! radius of step in front plate c rstep = 0.61 t1r = 0.024 ! thickness of rear endplate t1f = 0.024 ! thickness for r < rstep, front plate t2f = 0.012 ! thickness for r > rstep, front plate c t1f = 0.014 c t2f = 0.014 c delta = possible extra length of outer cylinder delta = 0.0031 c delta = 0. Eal = 7.0e10 ! Young's modulus for 6061aluminum c Eal = 7.7e10 ! Young's modulus for 8090 aluminum anu = 0.33 ! Poissons' ratio c values of holefac and stresscon chosen to match Hodgson's results holefac = 1./1.5 ! modulus reduction factor for holes c holefac = 1. stressfac = 4.1 ! stress concentration factor around holes c stressfac = 1. aload = 0. do i = 1,51 rmid(i) = 0.5*(bi(i+1) + bi(i)) area(i) = pi*(bi(i+1)**2 - bi(i)**2) E(i) = Eal anui(i) = anu t(1,i) = t1f ! front endplate if (bi(i) .gt. rstep) t(1,i) = t2f t(2,i) = t1r ! rear endplate qi(i) = 0. if (ihole(i). eq. 1) then c holes(i) = holefac stresscon(i) = stressfac awire = 0. do j = 1,3 loadfac(j,i) = loadfac(j,i)*ncelltot(i) ! convert to wires per layer ! calculate area of wires in this layer awire = awire + loadfac(j,i)*pi*dhole(j)**2/4. nwire(j) = nwire(j) + loadfac(j,i) oswire(j,i) = twire(j) ! oswire stores wire tension in grams temp = oswire(j,i)*loadfac(j,i)/1000. ! wire tension in N aload = aload + temp ! load in kg qi(i) = qi(i) + temp*9.8/area(i) ! pressure in Pa c activate the following line for no load conditions c qi(i) = 0. c write(1,9898) j, i, loadfac(j,i), oswire(j,i), temp 9898 format (' j, i',2i3,' loadfac, oswire, temp ',i5,2e11.3) enddo ! j ! set hole factor to fraction of area not in holes c holes(i) = 1. - awire/area(i) ! set hole factor to holefac holes(i) = holefac else holes(i) = 1. stresscon(i) = 1. endif enddo ! i write (1,1) zlen, n 1 format (' Stress analysis of a chamber ',f5.3,' m long, with ', . i3,' regions'/) c spring constants on inner and outer cylinders c Al + Be inner cylinder zAl = 1.340 zBe = 1.424 tAl = 0.005 tBe = 0.0011 EBe = 3.e11 scin = zAl/tAl/EAl + zBe/tBe/EBe c scin = 0. c carbon fiber outer cylinder tcf = 0.003 Ecf = 5.e10 scout = zlen/tcf/Ecf c scout = 0. write (1,91) scin, scout, delta 91 format (' spring const. per unit length, inner =, 'e11.3, . ' outer = ',e11.3/' outer cylinder longer by ',f7.4,' m'/) write (1,2) aload, nwire 2 format (' wire load = ',f5.0,' kg for ', i5,' sense, ',i5, . ' clearing, and ',i5,' field wires'// . ' Parameters for the regions:'// . ' region, b(cm), P(kPa), E(Gpa), nu, holefac, stress tf(mm),' . ' Df, tr, Dr') do i = 1,n do j = 1,2 Di(j,i) = E(i)*holes(i)*t(j,i)**3/12./(1. - anui(i)**2) ! plate constant enddo ! j write (1,11) i, bi(i)*100., qi(i)/1000., . E(i)*1.e-9, anui(i), holes(i), stresscon(i), . t(1,i)*1000., Di(1,i), t(2,i)*1000., Di(2,i) 11 format (i4,f10.1,f8.2,f8.2,f6.2,f7.3,f8.2,f7.1,f7.0,f7.1,f7.0) c write on unit 2 for TeX table write (2,112) i, bi(i)*100., (loadfac(j,i),j=1,3), . qi(i)/1000., holes(i), stresscon(i), t(1,i)*1000., Di(1,i) 112 format (i3,' & ',f4.1,' & ',i3,' & ',i3,' & ',i3,' & ', . f4.1,' & ',f4.2,' & ',f3.1,' & ',f3.0,' & ',f7.0,' \t') 10 format ( . ' region',t30,i3/ . ' inner radius',t30,1pg11.3,' m'/ . ' thickness, b < r < c',t30,g11.3,' m'/ . ' pressure',t30,g11.3,' Pascals'/ . ' Youngs modulus',t30,1pg11.3,' Pascals'/ . ' Poissons ratio',t30,g11.3/ . ' Modulus loss to holes',t30,g11.3/ . ' Stress conc. at holes',t30,g11.3/ . ' D = plate constant',t30,g11.3/) enddo ! i write (1,12) bi(n+1) 12 format ( . ' outer radius',t30,g11.3,' m'/) c build the transformation from inner to outer regions c front endplate call build(amatf,cvecf,1) c rear endplate call build(amatr,cvecr,2) c solve for values at inner radius, case 2c call case2c(vbegf,vbegr,amatf,cvecf,amatr,cvecr) c print final values, and store final displacements in dzf call output(1,vbegf,0,1,dzf) call output(1,vbegr,0,2,dzf) stop 9999 continue stop end subroutine build(amat,cvec,ifr) c build transformation between inner and outer radii c vec_j(1) -- y at inner edge of region j c vec_j(2) -- theta c vec_j(3) -- radial bending moment, M_r c vec_j(4) -- shear force density, Q c form of transformation of vec_j from the inner edge of region j to c vec_i at the outer edge of region j: c vec_i(k) = sum_l amat_j(k,l)*vec_j(l) + cvec_j(k), c where amat_j and cvec_j are obtained from routine smatrix c we also want to transform from the inner edge of region 1 to the c outer edge of region j: c vec_i(k) = sum_l amat_i(k,l)*vec_j(l) + cvec_i(k) c therefore c cvec_i(k) = sum_l amat_j(k,l)*cvec_i-1(l) + cvec_j(k), c amat_i(k,l) = sum_m amat_j(k,m)*amat_i-1(m,l). c I store cvec_i-1 in cvec, and amat_i-1 in amat c ifr = 1 -- front, = 2 -- rear implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta dimension amat(4,4), cvec(4), amatj(4,4), cvecj(4), amati(4,4), . cveci(4) do k = 1,4 cvec(k) = 0. do l = 1,4 amat(k,l) = 0. if (k .eq. l) amat(k,l) = 1. enddo ! l enddo ! k do j = 1,n ! get coeffs at outer edge of region call smatrix(j,bi(j+1),amatj,cvecj,ifr) c write (1,93) j, bi(j+1), cvecj, ((amatj(k,l),l=1,4),k=1,4) 93 format (' region ',i3,' outer radius ',f8.4/ . ' cvec',t10,4e11.3/' amat ',(t10,4e11.3)) do k = 1,4 cveci(k) = cvecj(k) do l = 1,4 cveci(k) = cveci(k) + amatj(k,l)*cvec(l) amati(k,l) = 0. do m = 1,4 amati(k,l) = amati(k,l) + amatj(k,m)*amat(m,l) enddo ! m enddo ! l enddo ! k do k = 1,4 cvec(k) = cveci(k) do l = 1,4 amat(k,l) = amati(k,l) enddo ! l enddo ! k enddo ! j return end subroutine case2c(vbegf,vbegr,amatf,cvecf,amatr,cvecr) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta dimension amatf(4,4), cvecf(4), vbegf(4), vfinf(4) dimension amatr(4,4), cvecr(4), vbegr(4), vfinr(4) dimension c(3), bmat(3,3), binv(3,3) c apply B.C.'s at inner and outer radii -- depending on case c case 2c, simply supported at both a and b with load transferred c between the two endplates by cylinders c deflection defined to be = 0 at inner radius of front plate vbegf(1) = 0. c simply supported edges => c bend moment = 0 at inner and outer radii of both plates vbegf(3) = 0. vbegr(3) = 0. vfinf(3) = 0. vfinr(3) = 0. c first solve for vbegf(2), vbegf(4) and vbegr(2) c load the constant vector c--- The commented out lines in c and bmat reduce the problem to that of a c--- single endplate (when restored) c(1) = cvecr(1) + cvecf(1) - scout*cvecf(4) - delta c c(1) = cvecf(1) - scout*cvecf(4) c(2) = -cvecf(3) c(3) = -cvecr(3) c c(3) = 0. c load the 3 by 3 matrix bmat(1,1) = -amatf(1,2) + scout*amatf(4,2) bmat(1,2) = -amatf(1,4) + scout*amatf(4,4) + scin*amatr(1,1) - . amatr(1,4) c bmat(1,2) = -amatf(1,4) + scout*amatf(4,4) bmat(1,3) = -amatr(1,2) c bmat(1,3) = 0. bmat(2,1) = amatf(3,2) bmat(2,2) = amatf(3,4) bmat(2,3) = 0. bmat(3,1) = 0. bmat(3,2) = -scin*amatr(3,1) + amatr(3,4) c bmat(3,2) = 0. bmat(3,3) = amatr(3,2) c bmat(3,3) = 1. write (1,92) bmat(1,1), bmat(1,2), bmat(1,3), c(1), . bmat(2,1), bmat(2,2), bmat(2,3), c(2), . bmat(3,1), bmat(3,2), bmat(3,3), c(3) 92 format ( .' M11 ',e11.3,' M12 ',e11.3,' M13 ',e11.3' C1 ',e11.3/ .' M21 ',e11.3,' M22 ',e11.3,' M23 ',e11.3' C2 ',e11.3/ .' M31 ',e11.3,' M32 ',e11.3,' M33 ',e11.3' C3 ',e11.3/) c invert the 3 x 3 matrix call minv3(bmat,binv) vbegf(2) = binv(1,1)*c(1) + binv(1,2)*c(2) + binv(1,3)*c(3) vbegf(4) = binv(2,1)*c(1) + binv(2,2)*c(2) + binv(2,3)*c(3) vbegr(2) = binv(3,1)*c(1) + binv(3,2)*c(2) + binv(3,3)*c(3) c complete the solution for the inner parameters vbegr(1) = -scin*vbegf(4) vbegr(4) = vbegf(4) c as a check, calculate the remaining outer paramters do i = 1,4 vfinf(i) = cvecf(i) vfinr(i) = cvecr(i) do j = 1,4 vfinf(i) = vfinf(i) + amatf(i,j)*vbegf(j) vfinr(i) = vfinr(i) + amatr(i,j)*vbegr(j) enddo ! j enddo ! i write (1,10) vbegf, vbegr, vfinf, vfinr 10 format (' vbeg(front) = ',4e11.3/' vbeg(rear) = ',4e11.3/ . ' vfin(front) = ',4e11.3/' vfin(rear) = ',4e11.3/) return end subroutine minv3(a,b) c b = inverse of 3 x 3 matrix a implicit real*8 (a-h), (o-z) dimension a(3,3), b(3,3) det = a(1,1)*a(2,2)*a(3,3) + . a(1,2)*a(2,3)*a(3,1) + . a(1,3)*a(2,1)*a(3,2) - . a(1,1)*a(2,3)*a(3,2) - . a(1,2)*a(2,1)*a(3,3) - . a(1,3)*a(2,2)*a(3,1) if (det .eq. 0.) then write (1,10) 10 format (' singular 3 x 3 matrix') stop endif ! det b(1,1) = (a(2,2)*a(3,3) - a(3,2)*a(2,3))/det b(1,2) =-(a(1,2)*a(3,3) - a(3,2)*a(1,3))/det b(1,3) = (a(1,2)*a(2,3) - a(2,2)*a(1,3))/det b(2,1) =-(a(2,1)*a(3,3) - a(3,1)*a(2,3))/det b(2,2) = (a(1,1)*a(3,3) - a(3,1)*a(1,3))/det b(2,3) =-(a(1,1)*a(2,3) - a(2,1)*a(1,3))/det b(3,1) = (a(2,1)*a(3,2) - a(3,1)*a(2,2))/det b(3,2) =-(a(1,1)*a(3,2) - a(3,1)*a(1,2))/det b(3,3) = (a(1,1)*a(2,2) - a(2,1)*a(1,2))/det return end subroutine case2d(vbeg,amat,cvec) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta dimension amat(4,4), cvec(4), vbeg(4), vfin(4) c apply B.C.'s at inner and outer radii -- depending on case c case 2d, fixed at inner radius, simply supported at outer radius c height = 0 at inner and outer radii vbeg(1) = 0. vfin(1) = 0. c slope = 0 at inner radius vbeg(2) = 0. c bend moment = 0 at outer radius vfin(3) = 0. c solve for vbeg(3) and veg(4) using c vfin(1) = 0 = amat(1,3)*vbeg(3) + amat(1,4)*vbeg(4) + cvec(1) c vfin(3) = 0 = amat(3,3)*vbeg(3) + amat(3,4)*vbeg(4) + cvec(3) denom = amat(1,3)*amat(3,4) - amat(3,3)*amat(1,4) vbeg(3) = -(cvec(1)*amat(3,4) - cvec(3)*amat(1,4))/denom vbeg(4) = -(amat(1,3)*cvec(3) - amat(3,3)*cvec(1))/denom return end subroutine output(icase,vbeg,iout,ifr,dz) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta dimension vbeg(4), vbegi(4,99), vfin(4), dz(2,99) character cases(2) / 'c', 'd' / pi = 4.*atan(1.) if (iout .eq. 0 .and. ifr. eq. 1) write(1,2) 2 format (' Results for Front Endplate') if (iout .eq. 0 .and. ifr. eq. 2) write(1,3) 3 format (' Results for Rear Endplate') if (iout .eq. 0) write (1,1) cases(icase) 1 format (/' case 2',a1// . ' Values at inner edges of regions' // . ' region, y(mm), theta(rad), M_r, Load(kg)') c find values at beginning of each region do j = 1,4 vbegi(j,1) = vbeg(j) enddo ! j do i = 1,n+1 if (i. gt. 1) . call values(vbegi(1,i),vbegi(1,i-1),bi(i),i-1,ifr) c write (1,91) (vbegi(j,i-1),j=1,4), (vbegi(j,i),j=1,4) 91 format (' vbeg ',4e11.3/' vfin ',4e11.3) alb = 2.*pi*bi(i)*vbegi(4,i)/9.8 if (iout .eq. 0) write (1,11) i, vbegi(1,i)*1000., vbegi(2,i), . vbegi(3,i), alb 11 format (i4,f10.3,f9.4,f9.1,f8.0) enddo ! i c print values at middles of regions that have wires if (iout .eq. 0) write (1,29) 29 format (/' Values at middle of regions that have wires'/) if (iout. eq. 0) write (1,21) 21 format(' layer, r(cm), y(mm), theta(rad), M_r(N-m/m), M_t,' . ' Q(N), sigma_r(MPa), sigma_t') do i = 1,n r = rmid(i) call values(vfin,vbegi(1,i),r,i,ifr) c store displacement; change sign!!!! dz(ifr,i) = -vfin(1) amt = vfin(2)*Di(ifr,i)*(1. - anui(i)**2)/r + . anui(i)*vfin(3) sigmar = 6.*vfin(3)/t(ifr,i)**2*1.e-6*stresscon(i) ! include stress sigmat = 6.*amt/t(ifr,i)**2*1.e-6*stresscon(i) ! concentration factor if (ihole(i). eq. 1 .and. iout .eq. 0) . write (1,30) layer(i), r*100., . vfin(1)*1000., vfin(2), . vfin(3), amt, vfin(4), sigmar, sigmat c write on unit 2 for TeX table write (2,31) i, r*100., . vfin(1)*1000., vfin(2)*1000., . vfin(3), amt, vfin(4)*2.*pi*bi(i)/9.8, sigmar, sigmat 31 format (i3,' & ',f5.1,' & ',f6.3,' & ',f6.1,' & ',f7.1,' & ', . f7.1,' & ',f7.0,' & ',f6.1,' & ',f6.1,' \t') enddo ! i return 9999 continue c print values for 100 radii do i = 1,101 r = bi(1) + 0.01*(i- 1)*(bi(n+1) - bi(1)) ! find which region we are in j = 0 100 j = j + 1 if (r .gt. bi(j+1)) go to 100 if (j .gt. n) j = n call values(vfin,vbegi(1,j),r,j,ifr) amt = vfin(2)*Di(ifr,j)*(1. - anui(j)**2)/r + anui(j)*vfin(3) sigmar = 6.*vfin(3)/t(ifr,j)**2*1.e-6*stresscon(j) ! include stress sigmat = 6.*amt/t(ifr,j)**2*1.e-6*stresscon(j) ! concentration factor write (1,30) j, r*100., vfin(1)*1000., vfin(2), vfin(3), . amt, vfin(4), sigmar, sigmat 30 format (1x,i3, f7.2,',',f8.3,',',f9.4,',',f9.1,',', . f9.1,',',f9.0,',',f9.2,',',f9.2) enddo ! i return end subroutine values(vfin,vbeg,r,i,ifr) implicit real*8 (a-h), (o-z) dimension vfin(4), vbeg(4), amat(4,4), cvec(4) call smatrix(i,r,amat,cvec,ifr) do j = 1,4 vfin(j) = cvec(j) do k = 1,4 vfin(j) = vfin(j) + amat(j,k)*vbeg(k) enddo ! k enddo ! j return end subroutine smatrix(i,r,amat,cvec,ifr) implicit real*8 (a-h), (o-z) common/data/ bi(99), Di(2,99), qi(99), anui(99), E(99), t(2,99), . holes(99), stresscon(99), n, rmid(99), ihole(99), . layer(99), scin, scout, delta common/nu/anu dimension amat(4,4), cvec(4) anu = anui(i) b = bi(i) D = Di(ifr,i) q = qi(i) amat(1,1) = 1. amat(1,2) = r*f1(r,b) amat(1,3) = r**2/D*f2(r,b) amat(1,4) = r**3/D*f3(r,b) amat(2,1) = 0. amat(2,2) = f4(r,b) amat(2,3) = r/D*f5(r,b) amat(2,4) = r**2/D*f6(r,b) amat(3,1) = 0. amat(3,2) = D/r*f7(r,b) amat(3,3) = f8(r,b) amat(3,4) = r*f9(r,b) amat(4,1) = 0. amat(4,2) = 0. amat(4,3) = 0. amat(4,4) = b/r cvec(1) = -q*r**4/D*f11(r,b) cvec(2) = -q*r**3/D*f14(r,b) cvec(3) = -q*r**2*f17(r,b) cvec(4) = -q*(r**2 - b**2)/2./r return end real*8 function f1(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f1 = 0.5*(1. + anu)*b/r*log(r/b) + . 0.25*(1. - anu)*(r/b - b/r) return end real*8 function f2(r,b) implicit real*8 (a-h), (o-z) f2 = 0.25*(1. - (b/r)**2*(1. + 2.*log(r/b))) return end real*8 function f3(r,b) implicit real*8 (a-h), (o-z) f3 = 0.25*b/r*((1. + (b/r)**2)*log(r/b) + (b/r)**2 - 1.) return end real*8 function f4(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f4 = 0.5*((1. + anu)*b/r + (1. - anu)*r/b) return end real*8 function f5(r,b) implicit real*8 (a-h), (o-z) f5 = 0.5*(1. - (b/r)**2) return end real*8 function f6(r,b) implicit real*8 (a-h), (o-z) f6 = 0.25*b/r*((b/r)**2 - 1. + 2.*log(r/b)) return end real*8 function f7(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f7 = 0.5*(1. - anu**2)*(r/b - b/r) return end real*8 function f8(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f8 = 0.5*(1. + anu + (1. - anu)*(b/r)**2) return end real*8 function f9(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f9 = b/r*(0.5*(1. + anu)*log(r/b) + 0.25*(1. - anu)* . (1. - (b/r)**2)) return end real*8 function f11(r,b) implicit real*8 (a-h), (o-z) f11 = (1. + 4.*(b/r)**2 - 5.*(b/r)**4 - 4.*(b/r)**2* . (2. + (b/r)**2)*log(r/b))/64. return end real*8 function f14(r,b) implicit real*8 (a-h), (o-z) f14 = (1. - (b/r)**4 - 4.*(b/r)**2*log(r/b))/16. return end real*8 function f17(r,b) implicit real*8 (a-h), (o-z) common/nu/anu f17 = 0.25*(1. - 0.25*(1. - anu)*(1. - (b/r)**4) - . (b/r)**2*(1. + (1. + anu)*log(r/b))) return end